A rapid and intelligent method for extracting multiple parameters of the underlying surface of a watershed

Abstract

This invention discloses a rapid and intelligent method for extracting multiple parameters of the underlying surface of a watershed, comprising: acquiring multimodal remote sensing data, meteorological observation data, and measured flow data of the target watershed; performing fusion and classification based on the remote sensing data before rainfall to obtain the underlying surface classification result and classification quality label, and determining the initial hydrological parameters accordingly; detecting changes in water level during rainfall based on real-time remote sensing data and model simulation to obtain inundation state labels; performing hydrophysical mechanism correction on the initial parameters according to the inundation state and meteorological data to obtain a physically corrected parameter field; combining the classification quality label and measured flow data, performing data assimilation correction constrained by physical state on the physically corrected parameter field, and finally outputting a dynamically corrected parameter field. This invention overcomes the shortcomings of traditional static and fixed model parameters, and realizes dynamic updating of hydrological parameters.

Description

Technical Field

[0001] This invention relates to the field of watershed flood control numerical simulation and hydrological information technology, and in particular to a rapid and intelligent method for extracting multiple parameters of the underlying surface of a watershed. Background Technology

[0002] In distributed hydrological and hydrodynamic simulations, the spatial refinement and temporal dynamic updating capability of underlying surface parameters, such as roughness, permeability, and soil moisture content, directly determine the computational accuracy of runoff generation and the timeliness of flood evolution forecasts. High-quality, high-frequency extraction of surface parameters is the core data foundation for the reliable operation of watershed flood control models and has significant engineering and technical importance for ensuring watershed flood control calculations and disaster prevention, mitigation, and early warning capabilities.

[0003] Currently, the mainstream methods for obtaining underlying surface hydrological parameters typically rely on static land cover data combined with empirical lookup tables for fixed assignment, or on post-interpretation using single optical remote sensing images. In actual operational use, single-modality remote sensing data is susceptible to interference from extreme weather conditions, such as persistent cloud cover, leading to large-area data gaps. Furthermore, traditional hydrological models generally use static parameters throughout the simulation period, failing to establish a dynamic response mechanism between the expansion of flood inundation areas during rainfall and sudden changes in surface hydrophysical properties. This results in the model exhibiting increasing systematic biases in simulating peak runoff and runoff timing when dealing with complex rainstorm and flood scenarios.

[0004] Existing parameter extraction and calculation systems face common technical challenges in dealing with high-intensity rainfall and complex geological changes, including insufficient data source robustness and delayed spatiotemporal dynamic response. There is a need to investigate a method that can effectively overcome extreme weather interference and achieve high-frequency dynamic transfer and accurate adaptation of surface parameters during flood evolution, thereby improving the adaptability and forecast accuracy of flood control numerical simulations under complex conditions. Summary of the Invention

[0005] The purpose of this invention is to provide a rapid and intelligent method for extracting multiple parameters of the underlying surface of a watershed, in order to solve at least one of the aforementioned problems in the existing technology.

[0006] Technical solution: A rapid and intelligent method for extracting multiple parameters of the underlying surface of a watershed, comprising:

[0007] Acquire multimodal remote sensing data, meteorological observation data, and measured flow data for the target watershed;

[0008] Based on multimodal remote sensing data before rainfall, feature fusion and classification reasoning are performed to obtain the underlying surface classification results of the target watershed and the corresponding classification quality label;

[0009] Based on the underlying surface classification results and the pre-configured baseline hydrological parameter field, the initial hydrological parameters of the target watershed are determined.

[0010] Based on multimodal remote sensing data acquired during rainfall and pre-acquired model simulation water levels, change detection is performed on the underlying surface classification results to obtain inundation status labels;

[0011] Based on the inundation status label and meteorological observation data, the initial hydrological parameters are corrected using hydrophysical mechanisms to obtain the physical correction parameter field;

[0012] By combining classification quality labels and measured flow data, a data assimilation correction constrained by physical state is performed on the physical correction parameter field, and a dynamic correction parameter field is output.

[0013] Beneficial effects: This invention overcomes the shortcomings of traditional model parameters being statically fixed, and realizes dynamic updating of hydrological parameters. Attached Figure Description

[0014] Figure 1 This is a schematic diagram of the overall process of a rapid and intelligent extraction method for multiple parameters of the underlying surface of a watershed, provided in an embodiment of this application.

[0015] Figure 2 This is a schematic diagram of the process of obtaining the underlying surface classification result of the target watershed by performing feature fusion and classification reasoning based on multimodal remote sensing data before rainfall, as provided in the embodiments of this application.

[0016] Figure 3 This is a schematic diagram of the process of obtaining the corresponding classification quality label based on feature fusion and classification reasoning of multimodal remote sensing data before rainfall, provided in the embodiments of this application.

[0017] Figure 4 This is a schematic diagram of the process of obtaining the inundation status marker by detecting changes in the underlying surface classification results based on multimodal remote sensing data acquired during rainfall and pre-acquired model simulation water levels, as provided in the embodiments of this application. Detailed Implementation

[0018] Example 1: A general flowchart of a rapid and intelligent method for extracting multiple parameters of the underlying surface of a watershed is provided, as follows: Figure 1 As shown, this paper elaborates on the overall architecture of multi-source data-driven dynamic parameter extraction and assimilation, and provides a technical system with spatial fine expression, high-frequency temporal response and physical consistency of parameters, which solves the technical problems of static solidification of underlying surface parameters and lag in real-time update in traditional flood numerical simulation.

[0019] Step 101: Acquire multimodal remote sensing data, meteorological observation data, and measured flow data of the target watershed.

[0020] Specifically, the multimodal remote sensing data includes medium-to-high resolution optical remote sensing images and synthetic aperture radar (SAR) images covering the target watershed. Optical remote sensing images can specifically utilize data from domestically produced high-resolution optical remote sensing satellites, such as the Jilin-1 series, to extract the spectral reflectance and texture features of ground features. SAR images can specifically utilize data from satellites such as the Water Resources-1 or Sentinel-1 to obtain information on surface backscattering intensity and polarization response. Furthermore, the multimodal remote sensing data can be supplemented with lidar point cloud data and UAV orthophotos. Lidar point cloud data is primarily used to generate high-precision digital elevation models to provide topographic benchmarks, while UAV orthophotos are used for local texture stitching and enhancement of predetermined boundary zones. Meteorological observation data includes real-time rainfall point observations from surface rain gauge networks within the watershed and quantitative precipitation estimation surface data from meteorological radar, with acquisition frequencies set to minute-level or hourly-level. Measured flow data originates from real-time hydrological monitoring data at the watershed outlet and major control sections. After acquiring the above data, orthorectification, geometric registration and radiometric calibration are performed uniformly to align the heterogeneous data sources to a unified spatial reference datum and resolution grid.

[0021] In some optional implementations, to ensure the effectiveness of subsequent calculations, the timestamp, spatial coverage, and effective pixel ratio of data acquisition are recorded for each modality during the data acquisition phase, and an availability metadata table is generated to provide objective criteria for diagnosing whether a modality is missing in subsequent phases.

[0022] Step 102: Based on the multimodal remote sensing data before rainfall, feature fusion and classification reasoning are performed to obtain the underlying surface classification results of the target watershed and the corresponding classification quality label.

[0023] In this embodiment, based on the constructed shared semantic space model, feature extraction and cross-modal fusion are performed on the aforementioned standardized and aligned multimodal remote sensing data. Specifically, considering the physical response characteristics of watershed runoff generation and confluence modeling, this embodiment predefines six typical underlying surface types, including water bodies, urban hardened areas, farmland, forest land, grassland, and bare land. The model outputs semantic classification results for these six underlying surface types at the patch scale, summarizing them to form underlying surface classification results covering the entire watershed. Simultaneously, to address potential modal loss scenarios in actual operations, such as optical imagery being obscured by clouds or radar imagery being outdated, the inference process generates corresponding classification quality labels for each spatial patch based on the modal completeness of the features extracted. The classification quality label is represented in data structure as a two-dimensional mask matrix consistent with the watershed grid scale. The values ​​of the elements in the matrix represent whether the corresponding area has undergone full-modal complete inference or single-modal degraded inference. Through the quality label mechanism, the confidence uncertainty in the pre-rainfall static base map extraction stage can be quantified and transferred to subsequent hydrological calibration and data assimilation stages.

[0024] Step 103: Based on the underlying surface classification results and the pre-configured baseline hydrological parameter field, determine the initial hydrological parameters of the target watershed.

[0025] In this step, it is necessary to establish the physical boundaries of relevant terms to eliminate ambiguity. Current soil moisture content is defined as a state variable that dynamically evolves with rainfall and infiltration, while water storage capacity is defined as a static parameter of the model characterizing the soil's maximum water-holding capacity. The pre-configured baseline hydrological parameter field is a priori parameter library derived from historical flood data after offline calibration. Based on the underlying surface classification results, corresponding parameter values ​​are extracted from the baseline hydrological parameter field according to the land use attributes of the land parcels and the basic soil texture information of the watershed. For example, urban hardened areas correspond to lower permeability and higher roughness of artificial impermeable surfaces, while forests and grasslands correspond to medium-high permeability and medium roughness. For the initial state of current soil moisture content, the surface soil moisture value obtained from the most recent synthetic aperture radar image before the start of rainfall is preferentially used, i.e., the previous water content state, for assignment. When remote sensing inversion is unavailable, measured values ​​from soil moisture monitoring stations within the watershed are used as a substitute through spatial interpolation. The previous water content state refers to the surface soil moisture content of the target area at the start of the rainfall event, and is used as a known initial condition input parameter for correction calculation. This yields spatially distributed initial hydrological parameters, which serve as the common starting state for subsequent rainfall-driven updates.

[0026] Step 104: Based on the multimodal remote sensing data acquired during rainfall and the pre-acquired model simulation water level, change detection is performed on the underlying surface classification results to obtain the inundation status label.

[0027] During continuous rainfall and flood evolution, some low-lying non-water bodies may be submerged, causing instantaneous abrupt changes in the underlying surface type. This step utilizes real-time synthetic aperture radar (SAR) imagery arriving during rainfall to identify areas of water expansion by detecting significant negative jumps in pixel backscattering intensity within the map patch area. Due to a several-hour observation blind period between two satellite transits, the simulated water level from the previous calculation step is compared with the digital elevation model (DEM). When the calculated water level exceeds the ground elevation of the map patch, it is predicted that the patch will be submerged. Through alternating verification between remote sensing measurements and model predictions, the classification attributes of the corresponding areas are dynamically updated. The output submersion status label not only includes a Boolean state indicating whether the area is submerged but also records time-varying attributes such as the inundation start time and the count of continuous submersion durations, used to distinguish between short-term transient submersion and long-term deep submersion.

[0028] In some alternative implementations, the flooding prediction at the current time step is strictly based on the model output results completed at the previous time step, and the model calculation at the current time step is performed based on the corrected parameters. This explicit time-progression mechanism avoids the circular dependency between parameter correction and hydrodynamic calculation.

[0029] Step 105: Based on the inundation status label and meteorological observation data, perform hydrophysical mechanism correction on the initial hydrological parameters to obtain the physical correction parameter field.

[0030] This step adjusts the parameters in the evolution process according to physical laws. For areas determined to be in a state of flooding, since the surface is covered by water and the underlying soil is nearly saturated, the saturated hydraulic conductivity is much smaller than the surface runoff. The system constrains its permeability parameter to a minimal constant approaching zero, no longer calculates the vertical infiltration component, and gradually transitions its roughness parameter to the characteristic value of pure water, i.e., the characteristic roughness value of water. For areas that have not been flooded, the cumulative rainfall index calculated from meteorological observation data, combined with the previous water content, drives the exponential infiltration decay calculation of the permeability parameter. Simultaneously, during prolonged periods of heavy rainfall, a controlled degradation reduction is applied to the water storage capacity parameter to reflect the decrease in water storage capacity caused by surface soil compaction or pore filling. Through these coordinated adjustments, a physically corrected parameter field conforming to the current hydrological, meteorological, and physical evolution laws is generated.

[0031] Step 106: Combining the classification quality label and measured flow data, perform data assimilation correction constrained by physical state on the physical correction parameter field, and output the dynamic correction parameter field.

[0032] Beyond simple physics-driven corrections, model outputs may still deviate from actual observations. This step introduces a data assimilation mechanism based on an ensemble Kalman filter framework. Specifically, the classification quality label from step 102 is obtained, and regions with high classification confidence are assigned smaller initial perturbation variances, while regions with low classification confidence and degraded inference are assigned larger initial perturbation variances, giving the data-driven algorithm prior spatial constraints. More importantly, for regions identified as inundated in step 104 and undergoing physical changes, their corresponding parameter perturbation variances are forcibly set to zero. According to the Kalman filter principle, the Kalman gain matrix obtained from parameter components with zero prior variance is also zero, achieving assimilation locking. This assimilation strategy, constrained by physical state, uses measured flow data to inversely correct parameter deviations in conventional regions, while protecting the physical correction rules for inundated regions from blind manipulation by purely data-driven algorithms. After assimilation and update, the optimal parameter set for the current time period, i.e., the dynamically corrected parameter field, is output, directly used to drive the flood control numerical simulation calculations for the next time step.

[0033] Example 2 details the offline construction process of a pre-configured baseline hydrological parameter field, solving technical challenges faced by high-dimensional spatial distributed hydrological models during parameter calibration, such as the curse of dimensionality, equivalence issues, and scale conflicts when jointly optimizing parameters of different dimensions.

[0034] Step 201: Extract the historical flood event database and decompose the hydrological parameters to be calibrated in each region into land type benchmark values ​​determined by the underlying surface type and soil texture, as well as spatial disturbance components characterizing local deviations.

[0035] Specifically, typical historical flood events covering different magnitudes and rainfall types are selected from historical hydrological and meteorological records and compiled into a historical flood event database. This database includes time-by-time areal rainfall processes, evaporation data, and measured flow sequences at control sections. In the patch-level parameterization framework of the distributed hydrological model, each patch has independent hydrological parameters. Directly calibrating all patch parameters one by one would lead to a sharp decrease in the search efficiency of the optimization algorithm in the ultra-high-dimensional parameter space. This embodiment constructs a two-layer dimensionality reduction encoding method for the parameter space to decompose and calculate each hydrological parameter to be calibrated in any region. Taking a single parameter as an example, its decomposition logic is implemented through the following formula:

[0036] θ i =β cs +Δ i ;

[0037] Where, θ i Let β be the composite value of the hydrological parameters to be calibrated in region i. cs The land use baseline value, Δ, is uniquely determined by the combination of underlying surface type c and soil texture type s in this region. i This is to characterize the spatial perturbation component in region i that deviates locally from the benchmark value of the same type.

[0038] In this embodiment, the hydrological parameters to be calibrated specifically include roughness parameters, permeability parameters, and water storage capacity parameters. All three parameters are decomposed into baseline and disturbance dimensions according to decomposition logic. It should be noted that the initial states of each flood event differ significantly. The current soil moisture content state variable before rainfall in each event is independently determined by station or remote sensing data and serves as the known initial condition driving the water balance equation; it is not included in the parameter decomposition and optimization scope of this step.

[0039] Step 202: Analyze the classification quality labels corresponding to the historical underlying surface classification results, and set differentiated perturbation boundaries for regions with different confidence levels to constrain the feasible domain of spatial perturbation components.

[0040] The classification quality identifier corresponding to the historical underlying surface classification result is obtained by performing the same classification and quality assessment process as in Example 3 on the historical multimodal remote sensing archive images within the coverage period of historical flood events, and is stored in the historical flood event database along with the historical flood-driven data.

[0041] Because the confidence levels of underlying surface classification results vary across different regions, the expected systemic bias introduced by the initial assignment also differs. This step utilizes the classification quality labels generated and transmitted during the pre-rainfall classification stage to set varying allowable ranges for the spatial perturbation components of each region. Specifically, for regions labeled with high confidence levels, the classification results are highly reliable, and parameter deviations mainly stem from generalization errors arising from empirical lookup tables; therefore, a narrower perturbation boundary is set. In this embodiment, the upper and lower limits of the perturbation are constrained within ±15% of the baseline value. For regions labeled with low confidence levels, the classification results undergo downgraded inference due to modality loss, posing a risk of type misclassification; therefore, a wider perturbation boundary is set. For example, the upper and lower limits of the perturbation are constrained within ±40% of the baseline value. This grants the optimization algorithm greater adjustment freedom to compensate for potential type misclassification.

[0042] Furthermore, since hydrological parameters such as roughness, permeability, and storage capacity are all strictly positive values ​​in a physical sense, non-negativity constraints must be added when setting differentiated disturbance boundaries. This non-negativity constraint relationship is expressed by the following inequality:

[0043] Δ i >=-p'*β cs ;

[0044] Where, Δ i β represents the spatial perturbation component, p' is the set maximum negative perturbation ratio limit and is strictly less than 1.0. cs This is the benchmark value for land use categories.

[0045] Through the above constraints, even under extreme negative perturbation conditions, the synthesized parameter values ​​are still greater than zero, ensuring the physical feasibility of the parameters at the algorithm level.

[0046] Step 203: Fix the spatial disturbance components, execute distributed hydrological simulation driven by the historical flood event database, and optimize the global objective function for the land use benchmark value.

[0047] Based on the dimensionality reduction encoding in the parameter space, the first layer of global optimization is implemented. Specifically, the spatial perturbation components of all regions are set as zero vectors, at which point the parameters of all regions are determined only by the land use baseline values ​​of the corresponding land use and texture combinations. The distributed hydrological model is driven by rainfall events from a historical flood event database. The simulated flow process output by the model is compared with the measured flow process to calculate the joint objective function for multiple stations. Commonly used low- to medium-dimensional continuous optimization algorithms are employed to solve this objective function; for example, hybrid evolutionary algorithms or particle swarm optimization algorithms can be used to search within a limited number of dozens of land use baseline value dimensions until the objective function converges.

[0048] In some alternative implementations, the joint objective function for multiple sites is constructed by weighting the overall process fit, the relative error of peak flow, and the relative error of flood volume. Its logic is expressed by the following formula:

[0049] J global =∑(α*(1-NSE)+w1*|RE peak |+w2*|RE vol |);

[0050] Among them, J global The objective function value is ∑, which represents the summation over all control sections and all calibrated flood events within the basin. α is the weighting coefficient for process fitting, NSE is the Nash efficiency coefficient, w1 is the error weighting coefficient for peak flow, and RE... peak RE represents the relative error of the peak flow, w2 is the error weighting coefficient of the flood volume, and RE vol This represents the relative error of the flood volume.

[0051] In actual operational configuration, to highlight the need for peak accuracy in flood forecasting, w1 can be set to 0.5, α to 0.3, and w2 to 0.2.

[0052] Step 204: After the land use benchmark value is optimized, the spatial disturbance components of each region are converted into dimensionless relative disturbance rates, and under the constraint of the differentiated disturbance boundary, a spatial smoothing regularization term is introduced to locally fine-tune the spatial disturbance components.

[0053] This step implements a second-level perturbation optimization to capture local biases in the hydrological response caused by differences in vegetation density and micro-topography. Since the roughness, permeability, and water storage capacity to be calibrated are different physical quantities, directly calculating their perturbation values ​​to account for spatial differences would lead to severe dimensional conflicts, causing the regularization term to lose its physical meaning. Therefore, dimensionless normalization is first performed on each perturbation value. This process is achieved through the following formula:

[0054] Δ tilde_ip =Δ ip / β csp ;

[0055] Where, Δ tilde_ip Δ is the dimensionless relative perturbation rate calculated for region i with respect to parameter type p. ip Let β be the spatial perturbation component of region i with parameter type p. csp For region i, the fixed land use benchmark value for parameter type p after optimization is completed.

[0056] After transforming to a dimensionless form, a local optimization objective function is constructed, comprising a hydrological fitting term and a spatial smoothing regularization term. A gradient optimization algorithm with boundary constraints is then used to iteratively optimize within the defined differential perturbation boundaries. The formula for the local optimization objective function is as follows:

[0057] J local =J hydro +∑(λ p *(Δ tilde_ip -Δ tilde_jp ) 2 );

[0058] Among them, J local To fine-tune the objective function value locally, J hydro To provide the hydrological fitting objective function value with the same form as the first layer, ∑ represents the double summation over all parameter types and all spatially adjacent region pairs, and λ p Δ is the predefined regularization intensity coefficient for parameter type p. tilde_ip Let Δ be the dimensionless relative perturbation rate of the central region i. tilde_jp Let be the dimensionless relative perturbation rate of adjacent region j.

[0059] To more clearly illustrate the dimensionless calculation process of the above regularization logic, the following simplified numerical example is provided. Assume that region i and its adjacent region j both belong to the farmland land category, and the optimized permeability benchmark value for both land categories is 20.0. In the current iteration, the permeability perturbation component of region i is 2.0, and the permeability perturbation component of region j is -1.0. The regularization intensity coefficient corresponding to permeability is set to 0.5. First, the dimensionless relative perturbation rate is calculated: for region i, it is 2.0 / 20.0 = 0.10, and for region j, it is -1.0 / 20.0 = -0.05. The penalty contribution of the squared difference between the two perturbations to the objective function is calculated as 0.5 * (0.10 - (-0.05)). 2 =0.5*0.0225=0.01125. This calculation process deviates from the original physical dimensions of permeability, and the magnitude of the penalty term is stable and controllable, effectively ensuring the smoothness of the update of the high-dimensional space perturbation field.

[0060] Step 205: Integrate the optimized land use benchmark values ​​with the locally fine-tuned spatial disturbance components to generate a pre-configured benchmark hydrological parameter field.

[0061] After the first-level baseline optimization and the second-level perturbation optimization are completed, the optimal land use baseline value is merged with the optimal spatial perturbation component constrained by the differentiated perturbation boundary using the inverse summation process of the aforementioned parameter decomposition logic. This process restores the absolute values ​​of various hydrological parameters for each region. These parameter values, along with their corresponding spatial coordinate indices and land use attributes, are then structured and stored in the database. This process ultimately generates a pre-configured baseline hydrological parameter field for online extraction and initial value assignment before rainfall occurs.

[0062] As an alternative or further improvement to this embodiment, after performing the local fine-tuning in step 204, the obtained global spatial disturbance components can be fixed, and the process can return to step 203 to perform a second round of global optimization on the land use benchmark values, iterating alternately between the two layers. When the improvement of the multi-site joint objective function in two adjacent iterations is lower than a preset minimum threshold, the iteration stops and the integration operation in step 205 is performed. The two-layer alternating iteration mechanism can further squeeze the fitting potential of the parameter space and improve the prediction accuracy limit of the benchmark hydrological parameter field.

[0063] Example 3 details the implementation process of the core algorithm for maintaining the integrity of the classification space and quantifying uncertainty when multimodal data is partially missing.

[0064] Step 301: Extract multimodal remote sensing training samples. During the forward propagation of the network, implement modal-level random masking on the input branch of the predetermined modality according to the preset probability distribution determined based on historical meteorological statistics during the flood season.

[0065] Specifically, during the offline training phase of the multimodal fusion classification model, a masking mechanism is introduced for the entire modal channel. The physical basis for setting the masking probability comes from historical meteorological statistics of the target watershed during the flood season. For example, during the flood season, the frequency of optical images becoming unusable due to cloud cover is significantly higher than the frequency of synthetic aperture radar (SAR) images failing. Based on this statistical pattern, the preset masking probability for the optical mode is set to be greater than that for the SAR mode. For example, the preset masking probability for the optical mode can be set to 0.4, and the preset masking probability for the SAR mode can be set to 0.1.

[0066] Step 302: When the modal-level random mask is triggered, all input channel features of the selected modal branch are set to zero to simulate the extreme situation in real hydrological observations where the corresponding modal data is completely unavailable due to cloud cover or time expiration.

[0067] In this step, the masking operation targets the entire modal data branch, rather than a single neuron node in the hidden layer. This design takes into account that in real-world applications, optical images often exhibit large areas of inefficiency, and traditional node-level feature discarding mechanisms are insufficient to simulate this overall modal-level failure. Specifically, the corresponding modal weight parameters in the mask matrix are all set to 0.0, cutting off all signals fed forward by that feature branch.

[0068] Step 303: Mix full-modal available samples and single-modal masked samples in the same training batch, and calculate the classification loss gradient simultaneously. This forces the independent feature encoders that are not set to zero as a whole to extract discriminative features independently under the condition of missing cross-modal complementary information, so as to obtain pre-trained model parameters that have both full-modal fusion reasoning and downgraded classification reasoning capabilities.

[0069] By mixing samples with different mask states in the same data batch, the model can simultaneously optimize the discriminative performance of both full-branch collaborative work and single-branch independent work when calculating the loss function and backpropagating to update the gradient. After training, the obtained multimodal fusion classification model can extract associated semantics using cross-modal attention mechanisms when data is complete, and can also seamlessly switch to degraded mode when a certain modality is missing, avoiding the technical problem of abnormal program exit due to the breakage of the input tensor dimension.

[0070] During training, the classification loss is calculated using the cross-entropy loss function. Those skilled in the art can use conventional model training methods, such as backpropagation combined with Adam or stochastic gradient descent optimizers, to train the model. The hyperparameters during training can be determined using conventional hyperparameter tuning methods.

[0071] Based on feature fusion and classification inference using multimodal remote sensing data prior to rainfall, the classification results of the underlying surface of the target watershed are obtained, such as... Figure 2 As shown, it includes the following steps:

[0072] Step 304: Using a pre-built independent feature encoder, extract the branch features of each modality in the multimodal remote sensing data.

[0073] During the online inference phase, separate independent feature encoders are invoked for optical and synthetic aperture radar modes, respectively. To accommodate the scale differences among various underlying land cover types, a hollow spatial pyramid pooling module is embedded in the backbone network of the independent feature encoders. This module aggregates information from multi-scale receptive fields to enhance the overall recognition capability of large-scale consistent categories such as grassland and woodland. Simultaneously, the network introduces a spatial residual connection structure, directly feeding shallow texture features to the deep semantic output to preserve the edge morphology of water bodies or hardened urban areas. LiDAR point cloud data and UAV orthophotos serve as auxiliary sources in this scheme. The former provides elevation benchmarks, and the latter provides local boundary stitching information. These are not directly used as regular inputs to the backbone classification network to avoid structural inference failures due to a lack of auxiliary observations in local areas.

[0074] Independent feature encoders can use encoder-decoder structures commonly used in the field as the backbone network, such as the DeepLab series architecture based on residual networks or the segmentation network architecture based on U-Net.

[0075] Step 305: Diagnose the spatial availability status of each modality data region by region.

[0076] For the acquired real-time multimodal remote sensing data, the effective observation quality within each grid area is calculated. For optical modes, it is determined whether the proportion of cloud-obscured pixels in the area exceeds a preset interference threshold, such as 50%. If it does, the optical mode in that area is considered missing. For synthetic aperture radar modes, the time difference between the acquisition timestamp of the most recent image and the current calculation time step is calculated. If this time difference is greater than a preset effective revisit period, the radar mode in that area is considered outdated.

[0077] Step 306: Based on the spatial availability, for regions with missing modal data, the feature signal paths of the corresponding missing modalities are masked, and downgraded classification reasoning is performed only based on the branch features of the remaining available modalities.

[0078] When the system diagnoses that a certain mode is unavailable, it triggers an adaptive gating mechanism within the model. This mechanism sets the gating weight parameter of the missing mode branch to 0.0, shutting down the signal transmission path of the corresponding independent feature encoder. The model relies entirely on the branch features extracted from the remaining available modes, and after shared semantic space mapping and classification head calculation, outputs the semantic label for the corresponding region. For example, in a region where only the synthetic aperture radar mode is available, the model extracts discriminative information solely based on the radar backscatter intensity.

[0079] Step 307: For regions with complete modal data, perform full modality fusion classification inference using the branch features of all modalities.

[0080] For regions effectively covered by data from all modalities, the model maintains all gate weight parameters at 1.0. Spectral reflectance and vegetation index features extracted from the optical modality, and polarization response features extracted from the synthetic aperture radar modality, are simultaneously projected into a shared semantic space. Within this space, a multi-head cross-modal attention mechanism is used to calculate the semantic correlation of features in each branch, and the output classification results are fused to achieve the highest classification accuracy. The multi-head cross-modal attention mechanism employs the scaling dot product attention calculation method, which is well-known in the field.

[0081] Step 308: Merge the inference outputs of each region to obtain the underlying surface classification results covering the entire target watershed.

[0082] The downgraded inference result generated in step 306 is geospatially stitched together with the complete inference result generated in step 307 according to the spatial coordinate index of each region. This combined processing method ensures that, under complex meteorological conditions, the watershed underlying surface classification map does not have blank areas in spatial coverage due to missing data, providing a spatially continuous base map for subsequent hydrological parameter extraction.

[0083] For example, during a flood season rainstorm, the upstream area of ​​the target watershed was continuously obscured by dense clouds. The proportion of cloud-obscured pixels in the optical remote sensing imagery of this area reached 85%, far exceeding the 50% interference threshold. The system determined that the optical mode was missing in this area. At this time, the system automatically shut down the signal path of the optical encoder in this area and performed downgraded inference based solely on the backscattering characteristics of the concurrently acquired synthetic aperture radar imagery, outputting the underlying surface classification result for this area and assigning it a second classification quality level. Meanwhile, the optical and radar data for the downstream area of ​​the watershed were both available, and full-modal fusion inference was performed, assigning it a first classification quality level. The inference results from the two areas were spatially stitched together to form a complete watershed classification base map, without any gaps due to the missing upstream data. In the subsequent data assimilation stage, the upstream second-quality-level area gained greater perturbation adjustment space to compensate for any biases that might be introduced by the downgraded inference.

[0084] Furthermore, feature fusion and classification reasoning are performed based on multimodal remote sensing data prior to rainfall to obtain corresponding classification quality labels, such as... Figure 3 As shown, it includes the following steps:

[0085] Step 309: Determine the actual reasoning path experienced by the corresponding region based on the available space status of each region.

[0086] After classification reasoning is completed, the system backtracks the combination of encoders activated at the input end during the calculation of each region. If the diagnostic record shows that both optical and radar encoders are active, it is determined that the region has undergone a full-modal fusion calculation path; if the diagnostic record shows that only a single encoder is active, it is determined that the region has undergone a single-modal degradation calculation path.

[0087] Step 310: For regions where full-modal fusion classification reasoning is performed, assign a first classification quality level representing high confidence.

[0088] For regions that have undergone a full-modal computational path, the probability of land cover misclassification is low because their classification criteria incorporate sufficient information from diverse and heterogeneous sources. The system assigns predetermined numerical identifiers to these regions, such as a value of 1, as a code for the first classification quality level, to declare to the downstream hydrological calculation module that the underlying surface classification results of this grid have high confidence.

[0089] Step 311: For regions where downgraded classification reasoning is performed, assign a second classification quality level representing low confidence.

[0090] For regions that have undergone a downgraded computation path, the classification results have potential uncertainty because single-modal inference inevitably loses some complementary information. The system assigns these regions another predetermined numerical identifier, for example, a value of 2, as a code for the second classification quality level, to warn downstream modules that the classification criteria for this grid may be biased.

[0091] Step 312: Integrate the first and second classification quality levels of each region to form a classification quality label that spatially corresponds to the underlying surface classification results.

[0092] The quality level identifiers of all regions are reorganized into a two-dimensional matrix to generate a quality mask matrix with the same resolution as the underlying surface classification results. This matrix serves as the classification quality identifier and is attached to the classification results as key metadata, which is then passed to the subsequent data assimilation and parameter calibration calculation modules, enabling the cross-domain transfer of data uncertainty from the classification stage to the hydrophysical extrapolation stage.

[0093] In some alternative implementations, to address the problem of missing modal data, in addition to implementing the masking and degradation strategy provided in this embodiment, an alternative modal reconstruction scheme can be introduced. Specifically, when a large area of ​​optical master modal data is missing, the system prioritizes retrieving UAV orthophotos covering the same area, and uses a regression algorithm to perform dimensionality reduction fitting on high-resolution multispectral bands, using them as alternative inputs for the optical branch. The masking operation of the characteristic signal path is only triggered when UAV data is also difficult to obtain. This alternative scheme further reduces the frequency of entering the degradation inference mode.

[0094] Example 4 details how to break the static assumption of the underlying surface by alternating verification of remote sensing measurements and model predictions during the flood evolution process, and how to use rainfall to drive the time-varying correction of physical parameters to prevent singularities, thus solving the technical problems of parameter response lag and mathematical calculation collapse in traditional hydrological models.

[0095] Based on multimodal remote sensing data acquired during rainfall and pre-acquired model simulation water levels, changes in the underlying surface classification results are detected to obtain inundation status labels, such as... Figure 4 As shown, it includes the following steps:

[0096] Step 401: Extract the temporal abrupt change features from the multimodal remote sensing data during the rainfall period.

[0097] Specifically, the surface conditions change drastically during rainfall and flooding. This step analyzes the acquired synthetic aperture radar (SAR) image data. According to the principle of electromagnetic wave scattering, when a region changes from a non-water body state to a water-covered state, the specular reflection from the water surface causes a significant negative jump in the radar backscattering intensity of that region. The difference in backscattering intensity between the current image and the pre-rainfall baseline image at each grid cell is extracted, and this difference sequence is defined as a temporal abrupt change feature.

[0098] Step 402: Compare the temporal abrupt change characteristics with the preset change threshold to identify the first type of inundation area where water body expansion has occurred.

[0099] Based on this, the system sets a trigger judgment mechanism according to empirical calibration. For example, the preset change threshold can be set to 3.0 dB. When the decrease in backscattering intensity extracted from a certain grid area is greater than or equal to 3.0 dB, the system determines that the surface of that area has been covered by water and marks it as a Class I flooded area. This judgment process is based on objective physical quantities from remote sensing observations and does not rely on subjective prior assumptions.

[0100] Step 403: During the observation blind period when no multimodal remote sensing data is acquired, compare the pre-acquired model simulation water level with the pre-configured regional elevation data to predict the second type of inundation area where water body expansion will occur.

[0101] Due to the limitations of satellite orbit revisit cycles, there may be observation blind periods of up to several hours between two consecutive remote sensing image acquisitions. To fill these data gaps in the time series, this embodiment introduces an explicit time-progressive prediction mechanism. Specifically, at the current calculation time step, the simulated water level from the model output in the previous calculation time step is directly retrieved and numerically compared with the regional elevation data pre-generated from the lidar point cloud. If the water level value of a certain region in the previous time step is strictly greater than the ground elevation value of that region, the region is predicted to be flooded and marked as a second-type flooded area. By relying on the known state of the previous time step to predict the state of the current time step, the deadlock cycle dependency between parameter correction and hydrodynamic calculation is eliminated from the algorithm's underlying layer.

[0102] In the first calculation step, since there is no model output from the previous step, the simulated water level can be initialized using steady-state hydraulic calculation results based on initial hydrological conditions or historical constant water level data. Subsequent steps use the actual output value of the hydrodynamic model from the previous step.

[0103] Step 404: Integrate the measured results of the first type of inundation area with the predicted results of the second type of inundation area, dynamically update the classification attributes of the corresponding areas in the underlying surface classification results, and generate inundation status labels.

[0104] The detection results from the two sources mentioned above are spatially stitched and logically merged. When the latest remote sensing data is available, the judgment result of the first type of inundation area is used to forcibly override the model's prediction result, forming a closed loop of alternating verification. The underlying attributes of the identified inundation areas are temporarily rewritten, and their classification attributes are changed to water body type at the current time step. The generated inundation status label contains a Boolean variable representing whether the area is inundated, and simultaneously records the timestamp of the first inundation of the area and the cumulative count of consecutive inundation periods to support subsequent differentiated physical corrections.

[0105] Step 405: Extract real-time rainfall point observations from the surface rain gauge network and quantitative precipitation estimation surface data from the meteorological radar from the meteorological observation data.

[0106] Before performing physical parameter corrections, a high-precision rainfall driving force field must be constructed. While data from ground-based rain gauge networks is highly accurate, it is spatially discrete and sporadic, whereas data from weather radar exhibits good spatial continuity but suffers from systematic estimation biases. The system simultaneously acquires heterogeneous data sources generated by these two independent observation devices, preparing for the fusion of spatial scale and numerical accuracy.

[0107] Step 406: Using real-time rainfall point observations as the true benchmark, calculate the spatial deviation correction factor between the observed data and the quantitative precipitation estimation surface data.

[0108] The system locates each ground rain gauge station in a spatial coordinate system and extracts the estimated precipitation values ​​from the meteorological radar grid at the corresponding coordinate locations. The ratio between the observed ground point values ​​and the estimated radar values ​​is calculated to obtain the local deviation correction coefficients for each station location. Using spatial inverse distance weighted interpolation or Kriging interpolation algorithms, the discrete local coefficients are interpolated and extended to the entire watershed grid, generating a continuously distributed spatial deviation correction factor matrix across the entire watershed.

[0109] Step 407: Use the spatial bias correction factor to perform global grid correction on the quantitative precipitation estimation surface data to generate a real-time fused precipitation field.

[0110] The raw quantitative precipitation estimation surface data output by the meteorological radar is multiplied pixel by pixel with the spatial deviation correction factor matrix obtained in step 406. After this multiplication calculation, the corrected precipitation field not only retains the spatial distribution pattern of precipitation detected by the radar, but also aligns with the actual ground observations in terms of magnitude, thereby outputting a high-precision real-time fused precipitation field.

[0111] Step 408 involves accumulating the real-time fused rainfall field over time to track the cumulative rainfall indicators generated in the spatial distribution of various regions.

[0112] At each time step of the model simulation, the real-time fused precipitation field value at the current time step is added to the value in the historical accumulation register. This accumulation operation continues from the start timestamp of the precipitation event definition to the current time, generating a dynamically growing two-dimensional cumulative rainfall map, where the value of each grid cell is the cumulative rainfall index for the corresponding region.

[0113] Step 409: Obtain the height of the surface vegetation canopy in the corresponding area and extract the inundation depth corresponding to the pre-acquired model simulation water level.

[0114] For areas marked as inundated, their surface roughness characteristics evolve with water depth. The system first retrieves the surface vegetation canopy height for the corresponding area from a pre-configured vegetation parameter library based on the underlying land cover attributes before rainfall. Simultaneously, it calculates the inundation depth for the current time step by subtracting the ground elevation of the corresponding grid from the simulated water level output from the previous time step.

[0115] Step 410: Construct a weighted transition function with the ratio of floodwater depth to surface vegetation canopy height as the independent variable.

[0116] To smoothly simulate the resistance changes during vegetation submersion, the ratio of submerged water depth to the height of the surface vegetation canopy is calculated. Based on this ratio, a weighted transition function is constructed with values ​​truncated between 0.0 and 1.0. When the water depth is shallow, the ratio is less than 1.0, and the function outputs the corresponding decimal weight; as the water depth increases until the vegetation canopy is submerged, the ratio is greater than or equal to 1.0, and the function outputs a truncated upper limit value of 1.0.

[0117] In this embodiment, the weighted transition function w is implemented in a truncated form, specifically as follows:

[0118] w=min(d inun / h veg ,1.0);

[0119] Where d inun h represents the current submerged water depth. vegThis represents the height of the ground surface vegetation canopy. Alternatively, other monotonically increasing transition functions with a range in the [0,1] interval, such as the S-shaped smoothing function, can be used instead.

[0120] Step 411: Using a weighted transition function, a dynamic weighted calculation is performed between the in-situ calibration roughness value corresponding to the initial hydrological parameters and the pre-stored roughness value of the pure water body characteristics. The calculated weighted average value is used as the transitioned roughness parameter. The in-situ calibration roughness value corresponding to the initial hydrological parameters is the roughness parameter in the initial hydrological parameters.

[0121] The weighted calculation process is as follows:

[0122] n t =w*n water +(1-w)*n land ;

[0123] Where, n t Here, is the roughness parameter after transition, w is the output value of the weighted transition function calculated based on the ratio, and n is... water For the pre-stored characteristic roughness values ​​of pure water, n land This represents the in-situ calibration roughness value corresponding to the initial hydrological parameters before rainfall. Through this dynamic weighted calculation, the model can accurately simulate the physical evolution process from shallow water obstruction to smooth flow in deep water.

[0124] Step 412: Based on meteorological observation data, calculate the cumulative rainfall index for each region.

[0125] It should be noted that the macroscopic functions of this step have been specifically described in steps 405 to 408. The system achieves accurate statistics on the driving indicators of rainfall by fusing ground observations and radar estimations.

[0126] Step 413: Based on the flooding status label, for areas determined to be in a state of water flooding, the roughness parameter in the initial hydrological parameters corresponding to the area is gradually transitioned to the pre-configured water feature value, and its permeability parameter is constrained to a preset minimum constant.

[0127] The preset minimum constant can be a minimum constant that approaches zero. In this embodiment, the minimum constant is set to 0.001.

[0128] This step formally writes the transition roughness parameters calculated in steps 409 to 411 into the parameter field of the current time step. Simultaneously, it performs crucial physical constraint operations. Based on objective physical laws, in flood-inundated areas, because the underlying soil is already saturated, its saturated hydraulic conductivity is typically only 1 to 5 millimeters per hour, while the runoff of surface floodwaters can reach tens to hundreds of millimeters per hour. The infiltration component has lost its dominant position in the water balance equation. The system directly forces the permeability parameters of the inundated area to be covered by extremely small constants such as 0.001. This constraint not only conforms to physical reality but also reduces the ill-conditioned iteration risk of the underlying partial differential equation solver.

[0129] Step 414: For areas not identified as being in a state of water inundation, perform time-varying attenuation correction on the permeability and water storage capacity parameters in the initial hydrological parameters corresponding to the area based on the cumulative rainfall index.

[0130] For areas experiencing rainfall but not yet large-scale flooding, the continuous infiltration of rainwater into the soil causes a slow degradation of the underlying hydrological properties, prompting the system to activate a parameter decay engine based on cumulative rainfall.

[0131] Step 415: Based on the cumulative rainfall index and the previous water content of the corresponding area, calculate the dynamic relative saturation to characterize the degree of soil saturation; wherein, the previous water content is the previous water content included in the initial hydrological parameters.

[0132] The system obtains the cumulative rainfall index tracked in step 408, and combines it with the soil base moisture state before rainfall, i.e., the previous soil moisture state, to calculate the dynamic relative saturation using the following formula:

[0133] S r =(W0+F cum ) / W m_base ;

[0134] Among them, S r For dynamic relative saturation, W0 represents the anterior water-bearing state variable of the corresponding region, i.e., the anterior water-bearing state of the corresponding region, and F... cum W is the current cumulative infiltration amount estimated based on the cumulative rainfall index. m_base The baseline values ​​for the pre-configured water storage capacity parameters for the region.

[0135] Among them, F cum This is the current cumulative infiltration estimated based on the cumulative rainfall index tracked in step 408. In practical applications, when the soil has not yet reached saturation, F can be approximated. cum Equals the cumulative rainfall, assuming a small proportion of runoff at the beginning of heavy rainfall; when S r When the value approaches 1.0, the singularity protection mechanism (step 417) will directly lock the permeability parameter, at which point F cumThe accuracy of the estimation has a negligible impact on the result. Those skilled in the art can also estimate F using conventional infiltration calculation methods such as the Green-Ampt model and the Phillips infiltration formula. cum The specific method can be selected according to the soil texture of the target watershed and the accuracy requirements of the model.

[0136] Step 416: When calculating the exponential infiltration decay process of the permeability parameter, set a threshold for determining the saturation extreme value.

[0137] Traditional infiltration attenuation calculation models tend to have a division-by-zero denominator in their exponential terms when the soil is near saturation, causing the program to throw a division-by-zero exception. To ensure the relative stability of the business system, this embodiment forcibly sets a numerical boundary at the outer layer of the algorithm; for example, the preset saturation extreme value judgment threshold can be set to 0.999.

[0138] Step 417: When the dynamic relative saturation reaches or exceeds the saturation extreme value judgment threshold, the singularity protection mechanism is triggered, skipping the attenuation calculation and directly locking the permeability parameter to the pre-configured stable infiltration rate lower limit value.

[0139] Real-time monitoring of dynamic relative saturation S r The numerical value of S. r When S ≥ 0.999, the conventional attenuation equation calculation process is immediately blocked, and the pre-configured lower limit of stable infiltration rate is directly output. r When the value is strictly less than 0.999, an improved calculation formula with protection against minimal constants is used for processing. The formula is expressed as follows:

[0140] ;

[0141] Among them, f t f is the permeability parameter calculated at the current time step. c Here, f0 is the pre-configured lower limit of stable infiltration rate, e is the natural constant, k is the preset infiltration decay rate coefficient, and S is the initial permeability parameter before rainfall. r The calculated dynamic relative saturation is given by ε, which is a minimal constant used to prevent division-to-zero anomalies when the limit approaches.

[0142] The infiltration attenuation rate coefficient k can be determined by those skilled in the art based on the soil texture and infiltration characteristics of the target watershed using conventional historical flood event inversion methods.

[0143] In this embodiment, the lower limit of stable infiltration rate f c Set according to soil texture type, for example, set clay to 1.0 mm / hour and loam to 3.5 mm / hour.

[0144] Step 418: When the cumulative rainfall index exceeds the preset multiple of the water storage capacity parameter, a reduction coefficient is applied to the water storage capacity parameter according to the amount of rainfall exceeding the limit, and a degradation lower limit is set to maintain the physical reasonable range of the parameter.

[0145] Continuous heavy rainfall not only alters soil moisture but also causes topsoil pore filling and structural compaction, leading to a degradation of the soil's maximum water storage capacity. This physical process is calculated using the following controlled reduction formula:

[0146] W m_t =max(κ min ,1-c d *(ρ-ρ0))*W m_base ;

[0147] Among them, W m_t The parameter represents the water storage capacity after degradation correction for the current period, where max is the maximum value function, and κ is the maximum value function. min For the set lower limit of degradation, c d The attenuation rate coefficient is ρ, which is the rainfall excess factor. It is defined as the ratio of the current cumulative rainfall to the baseline value of the regional water storage capacity parameter, i.e., ρ = F. rainfall_cum / W m_base , where F rainfall_cum This refers to the rainfall value corresponding to the cumulative rainfall index tracked in step 408. ρ0 is a preset multiple threshold, and W... m_base This is the baseline value for the water storage capacity parameter before rainfall.

[0148] The controlled degradation calculation mentioned above is triggered only when the rainfall exceeds a multiple ρ by a preset multiple threshold ρ0; when ρ ≤ ρ0, the water storage capacity parameter remains at the baseline value, i.e., W. m_t =W m_base The two logical statements above can be uniformly expressed in a segmented form as follows:

[0149] When ρ>ρ0, W m_t =max(κ min ,1-c d *(ρ-ρ0))*W m_base ;

[0150] Otherwise W m_t =W m_base .

[0151] To illustrate the controlled degradation logic described above, a dimensionless numerical example is provided. For simplified calculations, a baseline value W for the water storage capacity parameter is used. m_base =1.0, the rainfall calculated for the current period exceeds the multiple ρ by 2.0, the preset multiple threshold ρ0 is 1.5, and the attenuation rate coefficient c d Set to 0.10, lower degradation limit κ minThe value is set to 0.80. Substituting into the formula, the reduction multiplier is calculated as: 1 - 0.10 * (2.0 - 1.5) = 0.95. Comparison shows that the calculated reduction multiplier of 0.95 is greater than the lower limit of degradation of 0.80, so the value is taken as 0.95. Finally, the degenerate parameter W is calculated. m_t =0.95*1.0=0.95. This calculation process ensures that the parameter degrades reasonably with the intensity of the rainstorm, and never breaches the lower limit of the physical rationality of the underlying layer.

[0152] Step 419: Integrate the results of each region after gradual transition, constraint and time-varying decay correction to form a physical correction parameter field.

[0153] Before the end of each time step, the roughness of the flooded zone after the transition in step 411, the permeability of the flooded zone constrained in step 413, and the parameters of the non-flooded zone after time-varying decay in steps 417 and 418 are globally spatially summarized. The matrix generated above through rigorous physical filtering and singularity protection is the physical correction parameter field for the current time step. This parameter field will be used as a high-confidence prior state benchmark and fed downstream for further data assimilation and correction.

[0154] Example 5 details how to perform online data assimilation correction using measured flow data after obtaining the physical correction parameter field. This addresses the technical problems of traditional ensemble Kalman filtering in distributed hydrological model assimilation, such as the tendency to disrupt hydrological physical laws, the susceptibility to divergence in high-dimensional matrix operations, and the lack of uncertain prior constraints.

[0155] Step 501: Using the physical correction parameter field as the prior central reference, construct a set of state parameter joint vectors containing the model state variables and the parameters to be corrected in each region. This step can also specifically involve using the physical correction parameter field as the prior central reference to construct a set of state parameter joint vectors containing the pre-maintained model state variables and the parameters to be corrected in the physical correction parameter field for each region.

[0156] Specifically, to utilize observational data for feedback correction of model parameters, a standardized data assimilation infrastructure needs to be constructed. The system extracts the current model state variables for all grid regions, including the current soil moisture content and surface water depth for each region. Simultaneously, it extracts the hydrological parameters to be corrected for each region, including roughness and permeability parameters. The state variable arrays and parameter variable arrays of all spatial nodes are sequentially concatenated along the data dimension to form a high-dimensional one-dimensional array, which is the ensemble state parameter joint vector. Based on this, using the corresponding values ​​in the physical correction parameter field as the central reference value, multiple parallel ensemble members are generated by superimposing random perturbations. In this embodiment, the number of ensemble members is set to 80.

[0157] Step 502: Analyze the classification quality label and assign differentiated initial perturbation scales to the joint vector of the set state parameters of different regions based on the confidence level represented by the label.

[0158] Traditional methods typically employ a uniform global variance when generating perturbation components for set members. This embodiment utilizes classification quality labels generated during the classification downgrade inference stage for spatial differentiation constraints. The system parses the two-dimensional classification quality label matrix to identify the confidence level of each region. For regions labeled with the first classification quality level, it indicates that the underlying surface classification result is reliable and the bias of the prior parameter estimation is small, thus assigning a smaller initial perturbation scale. For regions labeled with the second classification quality level, it indicates that they have undergone downgrade inference with missing modalities and the potential bias of the prior parameter estimation is large, thus assigning a larger initial perturbation scale to provide the data assimilation algorithm with greater search and correction space. The calculation formula for this setting is as follows:

[0159] σ i =η q *θ prior ;

[0160] Where, σ i η is the initial perturbation scale set for region i with respect to predetermined parameters. q θ is the relative perturbation coefficient obtained by mapping the classification quality label to the corresponding level. prior These are the prior parameter values ​​determined by the physically modified parameter field for this region.

[0161] For example, the relative perturbation coefficient can be set to 0.05 for the first category quality level region and to 0.15 for the second category quality level region.

[0162] Step 503: Synchronously obtain the flooding status labels of each area.

[0163] After establishing a differentiated initial perturbation scale framework, a physical state verification mechanism is further introduced. By accessing memory or a database, the flooding state labels output in the earlier dynamic change detection are obtained grid by grid.

[0164] Step 504: For areas where the physical properties of the underlying surface have undergone qualitative changes due to water inundation, the parameter components in the area that have undergone strong numerical switching are marked as assimilation locked.

[0165] Based on the acquired flooding state labels, areas currently in a flooded state are identified. During the previous physical mechanism correction phase, the roughness parameters of these areas underwent a gradual transition towards water-related characteristics, and their permeability parameters were also constrained to a minimal constant approaching zero. This strong switching of parameter values ​​triggered by a qualitative change in the surface state reflects a deterministic physical fact. To prevent the purely mathematical statistical optimization of the data assimilation algorithm from overriding or tampering with the rigid correction based on physical rules, the system marks the index positions of the roughness and permeability parameters of the predetermined areas with predetermined Boolean values ​​in the data structure of the set state parameter joint vector; this is the assimilation locking state.

[0166] Step 505: For parameter components marked as assimilation locked, their initial perturbation scale is forcibly assigned to zero so that the corresponding covariance update gain is zero in subsequent assimilation update steps, thus protecting the executed hydrophysical mechanism correction results from modification by the assimilation process.

[0167] For parameter components marked with assimilation lock status, the system performs mandatory intervention when generating ensemble member perturbations, forcibly setting their initial perturbation scale to 0.0. In the matrix operation mechanism of ensemble Kalman filtering, if the variance of a parameter across all ensemble members is 0.0, then when calculating the prediction error covariance matrix, the covariance elements of the row and column containing that parameter must also be 0.0. Since the Kalman gain matrix is ​​directly derived from the prediction error covariance matrix, the update gain corresponding to that parameter will be forcibly calculated as 0.0. Through this underlying algebraic logic truncation, the correction amount for that parameter in the assimilation update equation is necessarily 0.0, achieving lock protection for deterministic physical states without altering the standard assimilation computation pipeline.

[0168] For example, at a certain point in the flood's evolution, a low-lying area in the middle reaches of the basin, formerly farmland, was submerged. During the physical correction phase, the permeability of this area was constrained to 0.001, and the roughness had transitioned to a water-specific value. Upon entering the data assimilation phase, the system detected that the submerged state label for this area was true and forcibly set the perturbation variance of its roughness and permeability parameters to zero. Although there were discrepancies between the measured flow data and the model forecast, the update increment of the parameters for this area was zero after the assimilation update was calculated using the Kalman gain matrix, and the physical correction results were fully preserved. Meanwhile, the upstream unsubmerged areas underwent normal data assimilation adjustments, and the slight deviation in permeability was corrected through measured flow feedback.

[0169] Step 506: For areas where no qualitative change in the physical properties of the underlying surface has occurred, random disturbances are generated normally according to the differentiated initial disturbance scale to accept data assimilation adjustments.

[0170] For non-submerged, non-qualitative regions not marked as being submerged by water, a Gaussian random number generator is invoked to generate normally distributed perturbation values ​​for each ensemble member, based on the differentiated initial perturbation scale calculated in step 502. These perturbation values ​​are then superimposed onto the prior central baseline to form a complete initial field for the ensemble members, providing a statistical sample containing reasonable physical uncertainties for subsequent calculations.

[0171] Step 507: Drive the pre-built hydrodynamic model to perform forward propagation to map it to the observation space, and calculate the forecast error covariance between parameters in each region.

[0172] Each generated set member is input into the watershed hydrological and hydrodynamic coupled model. The model performs a numerical extrapolation for one time step to calculate the state of the entire watershed at the next time step. The system extracts the simulated flow values ​​located at the control section coordinates from the output states of each member, forming the model observation vector. Using the state parameter vectors of each set member and their corresponding simulated observation vectors, the covariance matrix is ​​calculated to quantify the linear correlation between the parameter changes of each patch and the flow deviation at the control section, generating the forecast error covariance.

[0173] The hydrodynamic model can be a one-dimensional / two-dimensional coupled hydrodynamic model based on the Saint-Venant equations, or it can be the Xin'anjiang model based on the conceptual runoff generation and confluence equations, or other distributed hydrological models. The specific model type can be selected by those skilled in the art based on the watershed characteristics and accuracy requirements.

[0174] Step 508: When calculating the forecast error covariance, a localized truncation operator based on the pre-configured hydrological distance of the water network flow path is introduced. The localized truncation operator suppresses covariance elements that exceed the preset spatial influence radius to eliminate long-distance spurious correlations caused by finite set samples.

[0175] Since the number of set members is usually much smaller than the total dimension of the parameters to be corrected, the covariance matrix obtained by direct estimation inevitably contains sampling noise, leading to spurious statistical correlations between remote areas that have no direct hydraulic connection. Therefore, this step constructs a spatial distance filtering mechanism. Specifically, a localized truncation operator is constructed using a fifth-order piecewise polynomial function. The hydrological distance along the river network flow path between any grid region and the observation section is calculated, rather than the straight-line Euclidean distance. When this hydrological distance exceeds the set spatial influence radius, the localized truncation operator outputs 0.0, forcibly setting the corresponding covariance element to zero. For example, the spatial influence radius can be set to between 1.5 and 2.0 times the equivalent radius of the control basin of the observation section. The above calculation logic is implemented through the following matrix operations:

[0176] P loc =P f *C matrix ;

[0177] Among them, P loc P is the forecast error covariance matrix after localization suppression. f The initial forecast error covariance matrix is ​​calculated using ensemble samples. * represents the element-wise Hadamard product operation. C matrix It is a localized truncation operator matrix generated based on a fifth-order piecewise polynomial function and hydrological distance calculation.

[0178] Step 509: After the forward propagation of the joint vector of set state parameters is performed and before the assimilation update is entered, a pre-configured set expansion coefficient is applied to the perturbation components of set members that deviate from the mean, in order to compensate for the set divergence contraction phenomenon that may occur during the assimilation cycle.

[0179] As the forecast observation assimilation cycle progresses, the differences between ensemble members typically decrease gradually, leading to an underestimation of actual uncertainty in the covariance. This can cause filter divergence or rejection of new observation data. To maintain the ensemble's divergence distribution, the difference vectors of each member's deviation from the ensemble mean are extracted and multiplied by a pre-configured scalar strictly greater than 1.0, known as the ensemble inflation factor. For example, the ensemble inflation factor can be set between 1.02 and 1.10, effectively adding the amplified difference vector back to the ensemble mean.

[0180] Step 510: Using the measured flow data as the true observation input, perform a set uniform update based on the forecast error covariance, and extract the updated optimal estimate mean to assemble the dynamic correction parameter field.

[0181] The system acquires the actual monitored flow rate transmitted from the station at the current moment as the measured flow rate data. Combining the observation error covariance matrix with the previously processed forecast error covariance matrix, the optimal Kalman gain matrix is ​​calculated. The deviation between the simulated flow rate and the measured flow rate is multiplied by the gain matrix to calculate the analysis increment of the joint vector of state parameters. This analysis increment is superimposed onto the forward forecast field of the ensemble members to complete the assimilation update. The arithmetic mean of all updated ensemble member parameter components is calculated and extracted as the optimal estimated mean for the current time step.

[0182] Step 511: Extract the physical feasible domain boundary of each pre-configured hydrological parameter under the corresponding underlying surface type.

[0183] Although the assimilation process introduces covariance localization and classification quality constraints, purely mathematical optimization updates may still produce numerical results lacking physical meaning when dealing with extreme observational noise, such as negative roughness or permeability exceeding porosity. The system retrieves and obtains the physical theoretical lower and upper limits of each type of hydrological parameter from the basic database built during the offline calibration phase, serving as the boundaries of the physically feasible domain.

[0184] Step 512: Perform a region-by-region validity check on the extracted optimal estimated mean to determine whether the parameter values ​​that have shifted after assimilation and update are out of bounds.

[0185] Traverse the grid data structure of the entire watershed, and perform numerical comparison judgment on the optimal estimated mean of the aforementioned assimilation update output one by one.

[0186] Step 513: When the parameter value of the corresponding region exceeds the boundary of the physical feasible domain, it is forcibly truncated to the boundary extreme value. After truncation, it participates in the assembly to generate the final dynamic correction parameter field.

[0187] If the updated roughness parameter value for a certain region is found to be less than the set physical lower limit, the parameter value is directly rewritten to the physical lower limit value; if it is greater than the upper limit, it is rewritten to the upper limit value. This legality truncation operation ensures that all parameters entering the hydrodynamic model are within a safe calculation range. All legal parameters are merged with the state, and a dynamically corrected parameter field is encapsulated and output, completing the full data flow closed loop for a single time step.

[0188] In another alternative implementation, as a replacement for the boundary truncation method in step 513, a nonlinear mapping-based alternative method is provided to address the boundary overflow truncation problem described in step 513. Specifically, the system can pre-construct a Gaussian transformation operator. Before performing the assimilation update in step 510, a nonlinear monotonic mapping function, such as logarithmic transformation or arctangent transformation, is used to map the original parameter space with a defined physical boundary to an unbounded standard normal Gaussian space. Performing the assimilation update matrix operation within the unbounded Gaussian space avoids the truncation operation from disrupting the statistical properties of the set. After the update is complete, the corresponding inverse mapping function is used to transform the analysis field back into the physical parameter space. This alternative strictly avoids the risk of parameter overflow through mathematical space transformation, further maintaining the theoretical rigor of the filtering algorithm.

[0189] In a typical watershed flood control simulation verification, the method of this invention was used to retrospectively simulate historical rainstorm and flood events in a certain watershed. Compared with the traditional static parameter scheme, the method of this invention reduced the relative error of peak discharge simulation at the control section from 12.5% ​​to 6.8%, improved the Nash efficiency coefficient of the flood process from 0.82 to 0.91, and reduced the relative error of flood volume from 8.3% to 4.1%. Under the condition of missing modes, the overall accuracy of underlying surface classification in the downgraded inference region still remained above 83%, ensuring the continuity of the simulation process compared with the traditional scheme that caused the calculation to be interrupted due to missing data.

[0190] The experimental results of the above embodiments demonstrate that the method of the present invention significantly improves accuracy and robustness compared to traditional methods. Those skilled in the art will understand that the specific improvement may vary depending on the application scenario and dataset.

[0191] This application introduces an adaptive classification and quality label transfer mechanism for modal missing data. When images are severely obscured by clouds, it adaptively cuts off missing signal paths through low-level gating and seamlessly switches to degraded inference with available modalities. Simultaneously, it outputs classification confidence levels downstream, ensuring the continuity of underlying surface spatial resolution under extreme conditions and effectively avoiding flood control calculation stagnation caused by data gaps. It solves the problem of insufficient robustness caused by the easy loss of remote sensing data under complex weather conditions.

[0192] This application employs a dynamic reclassification and physical linkage correction method that alternates between remote sensing state perception and hydrodynamic prediction verification. By frequently tracking the expansion of inundated water bodies and triggering infiltration attenuation with extreme value protection and roughness-weighted transition, the static parameter fixation assumptions of traditional models are improved. This allows the hydrological extrapolation process to accurately adapt to the actual hydrodynamic evolution of the land surface, effectively reducing the risk of computational collapse caused by mathematical model singularities. It also solves the problem of delayed response of surface parameters during flood evolution.

[0193] This application utilizes a two-layer dimensionality reduction and dimensionless regularization method to extract high-dimensional baseline parameters, and applies physical state locking and prior quality constraints during the data assimilation process. By forcing the covariance update gain of regions experiencing physical qualitative changes to zero, the purely data-driven parameter correction is strictly limited to the boundaries allowed by hydrophysical laws and classification confidence levels, eliminating filter divergence and numerical out-of-bounds phenomena, and improving the stability and forecast accuracy of watershed flood control numerical simulation. It solves the problems of physical law violation and high-dimensional divergence that are easily caused during parameter extraction and assignment.

[0194] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A rapid and intelligent method for extracting multiple parameters of the underlying surface of a watershed, characterized in that, include: Acquire multimodal remote sensing data, meteorological observation data, and measured flow data for the target watershed; Based on multimodal remote sensing data before rainfall, feature fusion and classification reasoning are performed to obtain the underlying surface classification results of the target watershed and the corresponding classification quality label; Based on the underlying surface classification results and the pre-configured baseline hydrological parameter field, the initial hydrological parameters of the target watershed are determined. Based on multimodal remote sensing data acquired during rainfall and pre-acquired model simulation water levels, change detection is performed on the underlying surface classification results to obtain inundation status labels; Based on the inundation status label and meteorological observation data, the initial hydrological parameters are corrected using hydrophysical mechanisms to obtain the physical correction parameter field; By combining classification quality labels and measured flow data, a data assimilation correction constrained by physical state is performed on the physical correction parameter field, and a dynamic correction parameter field is output. Combining classification quality labels and measured flow data, a data assimilation correction constrained by physical state is performed on the physical correction parameter field, outputting a dynamic correction parameter field, including: Using the physical correction parameter field as the prior central benchmark, a set of state parameter joint vectors containing model state variables and parameters to be corrected for each region is constructed. The classification quality labels are analyzed, and based on the confidence level represented by the labels, different initial perturbation scales are assigned to the joint vector of the set state parameters of different regions. The pre-built hydrodynamic model is driven to propagate forward to map to the observation space, and the forecast error covariance between parameters in each region is calculated. Using measured flow data as the true observation input, a set uniform update is performed based on the forecast error covariance, and the mean of the updated optimal estimate is extracted to assemble the dynamic correction parameter field. Differentiated initial perturbation scales are assigned to the joint vector of ensemble state parameters for different regions, including: Simultaneously acquire flooding status labels for each region; For regions where the physical properties of the underlying surface undergo qualitative changes due to water inundation, parameter components in the region that experience strong numerical switching are marked as assimilation locked. For parameter components marked as assimilation locked, their initial perturbation scale is forcibly assigned to zero so that the corresponding covariance update gain is zero in subsequent assimilation update stages, thus protecting the already executed hydrophysical mechanism correction results from modification by the assimilation process. For regions where no qualitative changes have occurred in the physical properties of the underlying surface, random disturbances are generated normally according to the differentiated initial disturbance scale to accept data assimilation adjustments.

2. The method according to claim 1, characterized in that, Based on feature fusion and classification inference using pre-rainfall multimodal remote sensing data, the underlying surface classification results for the target watershed are obtained, including: Using a pre-built independent feature encoder, branch features of each modality in multimodal remote sensing data are extracted separately; Diagnose the spatial availability status of each modality's data region by region; Based on the availability of space, for regions with missing modal data, the feature signal paths of the corresponding missing modalities are masked, and downgraded classification reasoning is performed based on the branch features of the remaining available modalities. For regions with complete modal data, perform full-modal fusion classification inference using branch features from all modalities; By merging the inference outputs from each region, a classification result of the underlying surface covering the entire target watershed is obtained.

3. The method according to claim 2, characterized in that, Obtain the corresponding classification quality identifier, including: Based on the spatial availability status of each region, determine the actual reasoning path experienced by the corresponding region; For regions that perform full-modal fusion classification reasoning, assign a first classification quality level that represents high confidence. For regions that perform downgraded classification reasoning, assign a second classification quality level that represents low confidence. By integrating the first and second classification quality levels of each region, a classification quality label is formed that spatially corresponds one-to-one with the classification results of the underlying surface.

4. The method according to claim 1, characterized in that, Based on multimodal remote sensing data acquired during rainfall and pre-acquired model simulation water levels, changes in the underlying surface classification results are detected to obtain inundation status labels, including: Extracting temporal abrupt change features from multimodal remote sensing data during rainfall; By comparing the temporal abrupt change characteristics with preset change thresholds, the first type of inundation area where water body expansion occurs can be identified; During the observation blind period when multimodal remote sensing data is not available, the simulated water level of the model is compared with the pre-acquired regional elevation data to predict the second type of flooding area where water body expansion will occur. By integrating the measured results of the first type of inundation area with the predicted results of the second type of inundation area, the classification attributes of the corresponding areas in the underlying surface classification results are dynamically updated to generate inundation status labels.

5. The method according to claim 1, characterized in that, Based on the inundation status label and meteorological observation data, the initial hydrological parameters were corrected using hydrophysical mechanisms to obtain the corrected physical parameter field, including: Based on meteorological observation data, cumulative rainfall indicators were statistically analyzed for each region; Based on the submerged state label, for areas determined to be in a submerged state, the roughness parameter in the initial hydrological parameters is gradually transitioned to a pre-configured water feature value, and the permeability parameter is constrained to a preset minimum constant. For areas not identified as being in a state of water inundation, time-varying attenuation corrections are applied to the permeability and water storage capacity parameters in their initial hydrological parameters based on the cumulative rainfall index. The results of gradual transition, constraint and time-varying decay correction of each region are integrated to form the physical correction parameter field.

6. The method according to claim 5, characterized in that, Based on the cumulative rainfall index, time-varying attenuation corrections are applied to the permeability and water storage capacity parameters in the initial hydrological parameters, including: Based on the cumulative rainfall index and the previous water content of the corresponding area, the dynamic relative saturation degree used to characterize the degree of soil saturation is calculated. When calculating the exponential infiltration decay process of permeability parameters, a threshold for determining saturation extreme values ​​is set. When the dynamic relative saturation reaches or exceeds the saturation extreme value judgment threshold, the singularity protection mechanism is triggered, skipping the attenuation calculation and directly locking the permeability parameter to the pre-configured stable infiltration rate lower limit value. When the cumulative rainfall index exceeds the preset multiple of the water storage capacity parameter, a reduction factor is applied to the water storage capacity parameter based on the amount of rainfall exceeding the limit, and a degradation lower limit is set to maintain the physical reasonable range of the parameter.

7. The method according to claim 1, characterized in that, The process of calculating the forecast error covariance between parameters in each region and performing aggregated updates also includes: When calculating the forecast error covariance, a localized truncation operator based on the pre-configured hydrological distance of the water network flow path is introduced. The localized truncation operator suppresses covariance elements that exceed the preset spatial influence radius, thereby eliminating long-distance spurious correlations caused by finite set samples. After the forward propagation of the joint vector of set state parameters is performed and before the assimilation update is performed, a pre-configured set expansion coefficient is applied to the perturbation components of set members that deviate from the mean.

8. The method according to claim 1, characterized in that, The pre-configured baseline hydrological parameter field is constructed by performing the following offline parameter calibration steps in advance: Extracting a database of historical flood events, the hydrological parameters to be calibrated in each region are decomposed into land use baseline values ​​determined by both underlying surface type and soil texture, as well as spatial disturbance components characterizing local deviations; The classification quality labels corresponding to the historical underlying surface classification results are analyzed, and differentiated perturbation boundaries are set for regions with different confidence levels to constrain the feasible domain of spatial perturbation components. With fixed spatial disturbance components, distributed hydrological simulation is performed driven by a database of historical flood events to optimize the global objective function for land use benchmark values. After the land use benchmark value is optimized, the spatial disturbance components of each region are converted into dimensionless relative disturbance rates. Under the constraint of the differentiated disturbance boundary, a spatial smoothing regularization term is introduced to locally fine-tune the spatial disturbance components. By integrating the optimized land use benchmark values ​​with the locally fine-tuned spatial disturbance components, a pre-configured benchmark hydrological parameter field is generated.

9. The method according to claim 5, characterized in that, The roughness parameter in its initial hydrological parameters is gradually transitioned to the water body characteristic value, including: Obtain the height of the surface vegetation canopy in the corresponding area and extract the inundation depth corresponding to the pre-acquired model simulation water level; Construct a weighted transition function with the ratio of floodwater depth to surface vegetation canopy height as the independent variable; Using a weighted transition function, a dynamic weighted calculation is performed between the in-situ calibration roughness value corresponding to the initial hydrological parameters and the pre-stored characteristic roughness value of pure water area. The calculated weighted average value is then used as the roughness parameter after the transition.

Images