Flood evolution dynamic reenactment method for floodplain based on multi-source data fusion

Abstract

The present application relates to the technical field of water dynamics numerical simulation, and discloses a flood evolution dynamic reenactment method for a flood area based on multi-source data fusion, which comprises the following steps: obtaining a basic data set, generating an operation data set and a discrete hydrological observation sequence set; extracting a surface water flow connection network based on the operation data set to measure a connection distance and construct a distribution weight matrix; processing the observation sequence by using the weight matrix to execute assimilation and generate a physical constraint initial field data set; building a one-dimensional and two-dimensional coupled shallow water dynamics calculation framework containing Manning roughness coefficients and lateral weir flow coefficients, constructing a multi-objective evaluation function to complete parameter calibration; and setting a topographic feature boundary to execute forward numerical simulation and output quantitative evaluation results. The present application avoids assimilation of a water level lower than a surface elevation by introducing elevation data of a digital elevation model to execute physical bottom boundary constraint, and completes parameter calibration by using reverse optimization iteration, so as to convert a flow field sequence into a flood control index value to realize objective quantitative evaluation of comprehensive flood control benefits.

Description

Technical Field

[0001] This invention relates to the field of hydrodynamic numerical simulation technology, specifically to a method for dynamic simulation of flood evolution in floodplains based on multi-source data fusion. Background Technology

[0002] Numerical simulations of flood evolution in floodplains typically rely on discrete station observation data to construct the initial computational field. Traditional data assimilation interpolation methods assign weights solely based on spatial linear distances, neglecting the obstructive effect of complex topography on surface water flow paths. In traditional interpolation calculations, the spatial weight distribution deviates from the actual physical topology, causing the generated water level data to deviate from reality. This often results in the assimilated water level being lower than the elevation data from the digital elevation model, failing to provide a physically meaningful initial computational state for shallow water dynamics calculation frameworks.

[0003] When constructing a coupled one-dimensional and two-dimensional shallow water dynamics calculation model, the Manning roughness coefficient and the lateral weir flow coefficient directly affect the calculation accuracy. Currently, the model parameter calibration process relies heavily on repeated trial calculations and adjustments based on manual experience, lacking a systematic reverse optimization iterative mechanism. Existing calibration methods do not divide the discretized computational domain into different spatial functional areas and set multi-objective evaluation functions for each, and a single fitness criterion cannot comprehensively reflect the hydrodynamic characteristics within the complex space. Reliance on manual parameter adjustment leads to a lack of objectivity in the parameter calibration results, making the calibrated model unable to accurately reflect the dynamic evolution process of water flow in the floodplain.

[0004] After performing forward numerical simulations on multi-scenario geomorphic boundary features, it is necessary to evaluate the flood control benefits of the engineering project. Current evaluation techniques mostly remain at the qualitative description level, rarely extracting the instantaneous effective water depth and spatial planar area within the effective flood storage and detention grid set for integral calculations. Existing methods do not convert the flow field sequence output from numerical simulations into specific flood control index values, and lack a mechanism for obtaining a comprehensive flood control benefit score through weighted calculations. This results in a lack of quantitative data support for the evaluation results of the engineering flood control benefits under different geomorphic boundary conditions. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a dynamic simulation method for flood evolution in floodplains based on multi-source data fusion. This method solves the problems of traditional data assimilation and interpolation not incorporating topographic relationships, resulting in assimilated water levels being lower than surface elevations; hydrodynamic model parameter calibration relying on manual parameter adjustment, leading to an inability to objectively reflect the flow evolution process; and the lack of objective quantification of flood control benefits by converting flow field sequences into flood control index values.

[0006] To achieve the above objectives, the present invention provides the following technical solution: This invention provides a method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion, comprising the following steps: A basic dataset for the target area is obtained, which includes basic datasets of data sources obtained by different observation methods. Spatiotemporal alignment is performed on the data sources obtained by different observation methods to generate a basic computational dataset containing multiple spatial computing grids and a discrete hydrological observation sequence set containing multiple discrete hydrological observation stations. Based on the basic computational dataset, the surface water flow connectivity network is extracted, the hydrological connectivity distance between the spatial computing grid and the discrete hydrological observation station is calculated, and a spatial distribution weight matrix is ​​constructed. The spatial distribution weight matrix is ​​used to process the discrete hydrological observation sequence set to perform spatial assimilation, and a boundary mask function is introduced to constrain the calculation region of the interpolation to generate a physically constrained initial field dataset. A one-dimensional and two-dimensional coupled shallow water dynamics calculation framework containing Manning roughness coefficient and lateral weir flow coefficient is constructed by loading the physical constraint initial field dataset and dynamic boundary input conditions. A multi-objective evaluation function is constructed, and the parameters are calibrated by using reverse optimization to iteratively modify the Manning roughness coefficient and the lateral weir flow coefficient. Define the boundaries of multi-scene geomorphic features, call the one-dimensional and two-dimensional coupled shallow water dynamics calculation framework that has completed the parameter calibration to perform forward numerical simulation, extract the spatiotemporal evolution flow field sequence and output quantitative evaluation results.

[0007] The technical advantage of this invention lies in its departure from conventional interpolation methods based on Euclidean geometric distance. This scheme extracts the surface water flow connectivity network along the actual terrain and calculates the hydrological connectivity distance, ensuring that the spatially distributed weight matrix conforms to the physical topological laws of water flow evolution. During the data assimilation stage, physical bottom constraints are implemented using terrain elevation data and boundary mask functions to prevent the assimilated water level from falling below the surface elevation, thus providing an accurate initial field for the hydrodynamic model. Simultaneously, a multi-objective evaluation function and a reverse optimization algorithm are introduced to collaboratively calibrate one-dimensional and two-dimensional resistance coefficients, improving the accuracy of numerical simulations of flood evolution under various scenarios and supporting the scientific quantitative evaluation of flood control benefits.

[0008] Furthermore, in the steps of generating the basic computational dataset and the discrete hydrological observation sequence set, performing spatiotemporal alignment on the data sources acquired by the different observation methods includes: uniformly mapping the spatial topographic data of the target area, the boundary vector data of water-related projects, and the planar coordinates of the optical satellite remote sensing image sequence to the target reference projection coordinate system, and correcting the elevation node values ​​to a unified vertical elevation reference surface; setting a globally unified computational grid size and time calculation step size of the spatial computational grid in combination with the Courant number constraint conditions of hydrodynamic solution; extracting the time nodes corresponding to the rise points of typical historical flood events and setting them as global zero time; and reconstructing the continuous numerical sequence of the target monitoring location under the time calculation step using a linear interpolation function for the discrete hydrological observation sequence set with inconsistent sampling frequencies.

[0009] Furthermore, in the step of constructing the spatial distribution weight matrix, extracting the surface water flow connectivity network and calculating the hydrological connectivity distance includes: locking the isolated elevation minima of the digital elevation model in the basic computational dataset for local depression filling correction; extracting the adjacent grid directions with the largest spatial elevation gradient to establish a global flow direction matrix; and deriving a global confluence accumulation matrix; using the polygon with the largest historical flood water body coverage to prune the global confluence accumulation matrix; extracting grid node connections where the confluence accumulation is greater than the connectivity threshold to generate the surface water flow connectivity network; mapping the discrete hydrological observation stations as effective mapping anchors on the surface water flow connectivity network; using the connections between adjacent grid nodes as undirected weighted graph structure edges; performing graph theory pathfinding operations along the surface water flow connectivity network; and accumulating the spatial physical length of the undirected weighted graph structure edges on the shortest topological path to generate the hydrological connectivity distance.

[0010] Furthermore, in the step of constructing the spatial distribution weight matrix, the construction of the spatial distribution weight matrix includes: constructing a distance decay function that takes into account terrain obstacles and energy dissipation; using the hydrological connectivity distance, the energy decay index characterizing the surface roughness of the floodplain, and a preset minimum smoothing constant, calculating the initial weight of water level control for each of the effective mapping anchor points to the spatial computing grid; setting a local interpolation search radius; accumulating the initial weights of water level control corresponding to the effective mapping anchor points within the local interpolation search radius as the weight sum; dividing the initial weight of water level control for a single effective mapping anchor point by the weight sum to derive a standardized weight value that satisfies the sum of constants 1, which is used to generate the spatial distribution weight matrix.

[0011] Furthermore, in the step of generating the initial physical constraint field dataset, generating the initial physical constraint field dataset includes: establishing a two-dimensional spatial boundary mask function; using a ray crossing algorithm and the grid overlap area ratio to jointly determine whether the spatial computing grid falls within the two-dimensional closed space formed by the polygon covering the largest water body of the historical flood, for dividing the effective computing area and the invalid shielding area; for each spatial computing grid within the effective computing area, multiplying the standardized weight value of the effective mapping anchor point associated with the spatial computing grid with the corresponding measured water level value and performing accumulation to obtain the initial assimilated water level value; introducing the elevation data of the digital elevation model to perform physical bottom constraint verification; when it is detected that the initial assimilated water level value is less than the elevation data of the digital elevation model, forcibly correcting the water level value of the corresponding spatial computing grid to the elevation data of the digital elevation model superimposed with the ultrathin water layer constant, generating a water level surface processed by physical constraints.

[0012] Furthermore, in the step of generating the physical constraint initial field dataset, generating the physical constraint initial field dataset also includes: subtracting the elevation data of the corresponding digital elevation model from the water level surface after physical constraint processing to obtain the effective water depth value at the current time step, generating an effective water depth distribution matrix; calculating the effective water depth gradient between adjacent spatial computing grids, and when the effective water depth gradient exceeds a preset gradient threshold, constructing a spatial computing window containing the surrounding neighborhood for the spatial computing grid and performing mean smoothing filtering; forcibly setting the initial flow velocity component of the entire spatial computing grid to 0, and setting a no-slip solid boundary condition outside the effective computing region to establish an initial momentum field; and encapsulating the physical constraint water level surface, the effective water depth distribution matrix, and the initial momentum field into the physical constraint initial field dataset.

[0013] Furthermore, in the step of completing the parameter calibration, the construction of the one-dimensional and two-dimensional coupled shallow water dynamics calculation framework includes: extracting the main channel boundary line and arranging lateral connecting weirs along the main channel boundary line; topologically binding the spatial calculation grids of the one-dimensional channel nodes and the two-dimensional floodplain area; based on the explicit finite volume method, using the one-dimensional and two-dimensional water level difference to deduce the interface flow transfer flux, and setting an interactive flow limiter to limit the maximum overflow rate when the preset water level difference threshold is exceeded; establishing a two-dimensional continuity control equation to calculate the water mass conservation state of the flow field, and using the interface flow transfer flux of the lateral connecting weirs as the source term input; and establishing a momentum conservation equation to verify the physical constraints of water convection, water level gravity gradient, and bed friction on the momentum change of the water body.

[0014] Furthermore, in the step of completing the parameter calibration, constructing the multi-objective evaluation function and completing the parameter calibration include: dividing the discretized computational domain into spatial functional zones including the main channel floodplain, densely populated urban residential areas, and agricultural ecological floodplains based on the land use attributes of the underlying surface; setting the root mean square error index of the peak water level for the measured and simulated water levels for the densely populated urban residential areas, setting the relative error index of inundation duration for the agricultural ecological floodplains, setting the spatial intersection-combination error index of the simulated inundation range and the remotely sensed measured water area for the overall physical connectivity, and setting the Nash efficiency coefficient index of the flow process line for the main channel floodplains; combining the range normalization results of each index with the preset protection importance weights to form the multi-objective evaluation function; using the particle swarm optimization algorithm to perform swarm intelligent reverse optimization iteration on the Manning roughness coefficient and the lateral weir flow coefficient, using the convergence value of the multi-objective evaluation function as the fitness criterion, and extracting the optimal solution vector that meets the convergence verification requirements to complete the parameter calibration.

[0015] Furthermore, in the step of outputting the quantitative evaluation results, the forward numerical simulation includes: locating the spatial computing grid affected by engineering intervention based on the preset multi-scenario geomorphic feature boundary; updating and correcting the elevation data of the digital elevation model by introducing spatial location indicator variables and engineering design elevations; and performing local depression filling verification; loading the corresponding flow process line driving sequence into the one-dimensional channel node, and extracting the physical constraint initial field dataset as the global starting time base state for the calculation; calling the one-dimensional and two-dimensional coupled shallow water dynamics calculation framework that has completed the parameter calibration to perform forward integral calculation; monitoring the maximum flow velocity of the spatial computing grid across the entire domain and calculating the local Corant number; automatically halving the time step when the local Corant number exceeds the stability critical value; and adding a numerical artificial viscosity term to the spatial computing grid that has exceeded the boundary if the halved time step is lower than the minimum step threshold.

[0016] Furthermore, in the step of outputting the quantitative evaluation results, the output quantitative evaluation results include: filtering out the spatial computing grids below the effective water depth threshold based on the spatiotemporal evolution flow field sequence output by the positive numerical simulation, and using a hydraulic connectivity topology traversal mechanism to remove isolated spatial computing grids that are not physically connected to the main channel flood discharge area, thereby determining the effective flood storage and detention grid set; for each time step node, summing and integrating the product of the instantaneous effective water depth and spatial plane area of ​​each spatial computing grid in the effective flood storage and detention grid set to obtain the instantaneous effective flood storage and detention total volume, and extracting the highest value from the time series as the transient peak effective flood storage and detention total volume; comparing the numerical simulation output characteristics under different geomorphic scenarios, extracting the flood peak water level reduction, the increase in the transient peak effective flood storage and detention total volume, and the diffusion duration delay, and performing range normalization processing, and using the entropy weight method to calculate the comprehensive flood control benefit score.

[0017] This invention provides a method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion. It has the following beneficial effects: 1. This invention constructs a spatially distributed weight matrix by extracting the hydrological connectivity distance from the surface water flow network and introducing elevation data from the digital elevation model during the spatial assimilation stage to perform physical bottom constraint verification. This replaces the interpolation calculation based on pure spatial straight-line distance, making the spatial weight distribution conform to the physical topological laws of water flow and preventing the water level value obtained from the assimilation calculation from being lower than the surface elevation. This provides a physically meaningful initial computational state for the shallow water dynamics calculation framework.

[0018] 2. This invention establishes a one-dimensional and two-dimensional coupled shallow water dynamics calculation framework, combines spatial functional zone division to establish a multi-objective evaluation function, and uses a reverse optimization iteration method to modify the Manning roughness coefficient and the lateral weir flow coefficient and complete parameter calibration. This reduces the reliance on manual parameter adjustment experience in the model calibration process, combines multiple error indicators into a unified fitness criterion, and enables the calibrated model to objectively reflect the dynamic evolution process of floodplain water flow.

[0019] 3. This invention performs forward numerical simulation by setting boundaries of geomorphic features in multiple scenarios. It multiplies the instantaneous effective water depth of each spatial calculation grid within the effective flood storage and detention grid set with the spatial planar area and integrates the result to obtain the transient peak effective flood storage and detention volume. It then uses the entropy weight method to calculate the comprehensive flood control benefit score and transforms the flow field sequence output by the numerical simulation into flood control index values, thereby achieving an objective quantitative evaluation of the flood control benefits of the project under different geomorphic boundary conditions. Attached Figure Description

[0020] Figure 1 This is a flowchart of the dynamic simulation method for flood evolution in floodplains based on multi-source data fusion, as described in this invention. Figure 2 This is a diagram illustrating the data flow and collaboration relationships between modules in this invention. Figure 3 This is a flowchart of the multi-source basic spatial and hydrological data acquisition process of the present invention; Figure 4 This is a flowchart illustrating the extraction process of the real surface water flow connectivity network in this invention. Figure 5 This is a flowchart of the physical constraint initial field boundary condition transformation process of the present invention; Figure 6 Flowchart for setting the multi-objective evaluation function for spatial functional differences in this invention; Figure 7 This is a line graph comparing the flow process curves at the main channel control section of the present invention. Figure 8 This is a bar chart for the quantitative comparative evaluation of the flood storage and detention functions of this invention. Detailed Implementation

[0021] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] See attached document Figure 1 This invention provides a method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion, comprising the following steps: S10, the basic data acquisition module collects spatial topographic data, water-related engineering boundary vector data, remote sensing image interpretation data, land use type data or national spatial planning vector data of the target area, and hydrological observation sequences. Multidimensional data acquisition integrates data sources obtained through different observation methods. Spatial topographic data, combined with remote sensing extraction results and water level and flow measurement records, constitutes the basic computational dataset. The preprocessing process eliminates time scale and spatial coordinate differences between different data types and aligns the underlying data reference surface. S20, the topology weighting module, calculates the hydrological connectivity within the study area based on spatial topographic data and remote sensing image interpretation data. It extracts the actual surface water flow path network using the principle of hydrodynamic spread and transforms this network into an undirected weighted graph structure. Considering the positional relationship between the spatial computation grid and various ground observation stations, it extracts bidirectional reachable connectivity distances based on the water flow paths. A spatially distributed weight matrix is ​​established based on the law of physical energy decay. The connectivity distance parameter directly measures the actual resistance effect of water flow along complex terrain surfaces. S30, the initial field generation module uses a spatially distributed weight matrix to process discrete hydrological observation sequences. A boundary mask function is introduced, specifying the maximum inundation range. This boundary mask function constrains the spatial range of the water level numerical interpolation calculation, blocking the path of numerical calculation divergence into uninundated areas. A continuous water level distribution state field for the floodplain is generated sequentially over time. The output continuous water level distribution field constitutes the initial field and dynamic boundary input conditions for hydrodynamic calculations that meet the requirements of continuous physics. Here, a clear distinction is made: the global water level, depth, and momentum fields output by the initial field module only represent the initial state at the simulation start time (global zero time); while the dynamic boundary input conditions specifically refer to the upstream flow process line and downstream water level process line substituted into the model. The spatially interpolated water level field generated in subsequent time steps is only used as a model verification reference and not as a fixed spatial hard constraint input for positive forced driving. S40, the model building and verification module loads the initial field and dynamic boundary input conditions to build a one-dimensional and two-dimensional coupled shallow water dynamics calculation framework. It extracts channel water conveyance characteristic parameters and regional inundation area parameters to construct a multi-objective evaluation function. It performs reverse optimization iteration by combining the measured time series set and spatial image overlap. The reverse optimization iteration modifies the distribution of the Manning roughness coefficient and the lateral weir flow coefficient, completing the parameter calibration operation within the calculation framework. The S50 dynamic simulation and evaluation module sets up three types of boundary elevation base maps: natural state geomorphology, construction state geomorphology, and post-construction engineering geomorphology. Depending on the evaluation purpose, historically measured hydrological processes are input when used for historical reenactment, and the design return period flood process is input when used for engineering flood control capacity assessment, to conduct numerical simulation calculations. Data on inundated area, water depth distribution, water retention duration, and channel diversion ratios are extracted. The output index results are compared under different geomorphic conditions to quantitatively calculate the hydrological diversion changes caused by construction intervention, and to output quantitative conclusions on the regional water distribution pattern.

[0023] See attached document Figure 2 This invention provides a method for dynamic simulation of flood evolution in floodplains based on multi-source data fusion, comprising the following module collaborative flow process: The basic data acquisition module transmits the extracted digital elevation model (DEM) file and the historical flood maximum water coverage polygon to the topology weighting construction module. The topology weighting construction module uses the DEM file to perform surface runoff analysis, extracts surface connectivity network paths, and calculates path connectivity resistance using the historical flood maximum water coverage polygon, generating a spatial distribution weight matrix. Simultaneously, the basic data acquisition module sends the discrete hydrological observation sequence set to the initial field generation module.

[0024] The initial field generation module receives the spatial distribution weight matrix and performs matrix multiplication operations in combination with the discrete hydrological observation sequence set. The initial field generation module outputs a continuous field of water level distribution across all grid nodes. This generated continuous field of water level distribution is then passed to the model building and verification module, serving as the starting background condition input parameter for solving the shallow water dynamic equations and as auxiliary comparison data for evaluating the spatial overlap error of the subsequent model.

[0025] The basic data acquisition module synchronously inputs the retained measured flow sequence and observed water body area contour into the model construction and verification module. The model construction and verification module compares the measured flow sequence with the model's internal output flow, and adjusts the Manning roughness coefficient and lateral weir flow coefficient using a reverse optimization iterative mechanism. The model construction and verification module completes the parameter locking operation for the one-dimensional coupled shallow water dynamics calculation framework, and packages the one-dimensional coupled shallow water dynamics calculation framework with locked parameters for transfer to the dynamic simulation and evaluation module.

[0026] The dynamic simulation and evaluation module receives a one- or two-dimensional coupled shallow water dynamics calculation framework with locked parameters, and updates the topographic elevation attributes, Manning roughness partitions, and water-blocking entity boundaries of the bottom-level grid cells of the framework. The module then calls the one- or two-dimensional coupled shallow water dynamics calculation framework with locked parameters to perform forward evolution calculations, recording the water depth evolution values ​​and flux transfer status of each grid node. Finally, the module summarizes the operational output data under various topographic elevation boundary conditions, extracting cross-sectional flow distribution statistics and the calculated effective flood storage and detention volume.

[0027] See attached document Figure 3 In this embodiment, the basic data acquisition module performs the basic data acquisition operation, which includes the following steps: S11, the basic data acquisition module collects spatial topographic data of the target area and boundary vector data of water-related projects. The multi-dimensional data acquisition operation retrieves the digital elevation model (DEM) file generated by airborne lidar scanning and mapping. The DEM file contains gridded ground elevation node values, forming a bottom-layer three-dimensional physical space depicting the undulating state of the natural landform. The basic data acquisition module simultaneously imports the artificial building boundary coordinate file output by computer-aided design software and land cover type vector data obtained from remote sensing image interpretation. The artificial building boundary coordinate file covers three-dimensional vector polygons of three types of water-related projects: canal entities, flood control embankments, and construction cofferdams. The basic data acquisition module uses a vector rasterization algorithm to convert the three-dimensional vector polygons of the aforementioned water-related projects into grid data with the same resolution as the original DEM, and updates the elevation attributes of the corresponding grids according to the polygon design elevations, thereby constructing a spatial elevation base map with composite artificial intervention conditions.

[0028] S12, the basic data acquisition module processes multi-temporal satellite observation images to delineate the maximum visible inundation range during historical flood events. Multi-dimensional data acquisition operations acquire multi-period optical satellite remote sensing image sequences covering the target area. For cloud obstruction in the images, multi-temporal pixel synthesis technology is used to prioritize extracting cloud-free pixels to fill obstructed areas. The basic data acquisition module calculates an improved normalized differential water index for the multi-period optical satellite remote sensing image sequences, separating the pure water pixel set. Conventional normalized water indices are easily interfered with by background noise from buildings in complex engineering areas. Based on this phenomenon, this step uses an improved index that includes the shortwave infrared band. Utilizing the physical mechanism of strong absorption of shortwave infrared light by water, it suppresses the reflection characteristics of bare land and buildings around water-related projects, thereby highlighting the actual water body information.

[0029] The basic data acquisition module extracts the outer contour boundary of the pure water body pixel set to generate the polygon representing the maximum water body coverage during historical floods. By overlaying thermal infrared and visible light band data, the basic data acquisition module employs a joint judgment logic of brightness threshold and thermal infrared low temperature threshold (a pixel is judged as a cloud when it simultaneously meets the characteristics of high brightness in visible light and low temperature in thermal infrared), further eliminating pseudo-water body pixels in cloud-covered areas and mountain shadow areas. The calculation formula for the improved normalized differential water body index is configured as follows: ; In the formula, The numerical value of the improved normalized differential water index representing the target pixel; This represents the reflectivity value received by the target pixel in the green light band. This represents the reflectivity value of the target pixel in the shortwave infrared band.

[0030] When separating a set of pixels from a pure water body, a classification threshold needs to be set for the determination. This classification threshold is extracted from the pixel histogram distribution of the target area image, typically using the bimodal histogram method or Otsu's algorithm, and fluctuates around a value of 0. When a single pixel... A pixel whose calculated result is greater than the threshold is identified as a pure water pixel.

[0031] S13, the basic data acquisition module retrieves historical water level and flow records from relevant hydrological observation stations within the study area. This module accesses discrete hydrological observation sequence sets transmitted from ground-based flow measurement terminals and water level gauge sensors. These discrete hydrological observation sequence sets contain cross-sectional flow rates and elevation / water level values ​​output by each monitoring location at each measurement instant. The multidimensional data acquisition operation extracts continuous test records of the entire process of typical historical flood events based on timestamp labels. If the sensor data return frequencies from different monitoring stations are inconsistent, the sampling frequency of the main channel hydrological control station is used as a basis, and a linear interpolation algorithm is employed to fill in missing time nodes, establishing the hydrological boundary sequence conditions driving the numerical simulation process.

[0032] In this embodiment, the basic data acquisition module performs heterogeneous data spatiotemporal reference alignment preprocessing, which includes the following steps: In step S14, the basic data acquisition module performs spatial coordinate projection transformation and elevation datum unification operations. The multi-source basic data are initially distributed across different coordinate frames. The module selects the National Geodetic Coordinate System 2000 as the plane datum and the 1985 National Elevation Datum as the vertical elevation datum. For geotagged optical satellite remote sensing image sequences and man-made building boundary coordinate files, the module uses a seven-parameter coordinate transformation model to uniformly map the plane coordinates to the National Geodetic Coordinate System 2000. For raster layers with different spatial resolutions, the module sets a globally unified computational grid size based on the Courant number constraints obtained from subsequent hydrodynamic model solutions.

[0033] After determining a uniform computational grid size, the basic data acquisition module resamples the continuously distributed digital elevation model using bilinear interpolation and the discretely distributed remote sensing water body boundary grid using nearest neighbor interpolation to ensure a smooth transition of topographic data and faithful reproduction of water-related boundaries. For elevation values ​​from different sources, the basic data acquisition module corrects the associated elevation node values ​​to the 1985 National Height Datum based on the regional geoid control point network, ensuring seamless vertical alignment between the topographic base map and water-related engineering structures.

[0034] S15, the basic data acquisition module performs time scaling synchronization and step size discretization. When solving shallow water partial differential equations, the hydrodynamic model requires continuous boundary inputs with a fixed step size in the time dimension to maintain the stability of the numerical calculation. There are differences in the time base recorded by hydrological observation sequences and optical satellite remote sensing image sequences. The basic data acquisition module extracts the time nodes corresponding to the initiation points of typical historical flood events and sets these time nodes as the global zero time in the simulation process. The basic data acquisition module sets a unified time calculation step size according to the convergence requirements of the hydrodynamic model. For discrete hydrological observation sequence sets with inconsistent frequencies, the basic data acquisition module uses a linear interpolation function to reconstruct the continuous numerical sequence of the corresponding monitoring locations under a unified time calculation step size. The interpolation function formula is set as follows: ; In the formula, Represents the target time step after reconstruction. The water level value is calculated in real time; This represents the target time step variable after discretization. Represents the time step from the target The most recent preceding observation time point; Represents the time step from the target The most recent post-observation time point; This represents the previous observation time node. The actual measured water level value; Represents the later observation time node The actual measured water level value.

[0035] The application of this interpolation formula must satisfy the following conditions: Time interval constraints. When the target time step... When the data exceeds the coverage of the measured sequence, the basic data acquisition module maintains the boundary condition values ​​at the measured values ​​of the nearest endpoint to avoid data distortion caused by extrapolation calculations.

[0036] After completing time-scale synchronization, the basic data acquisition module matches the time label corresponding to the instant of imaging in the optical satellite remote sensing image sequence to the nearest time calculation step node, forming the specific time trigger anchor point for the spatial constraint boundary conditions to take effect. Through the above-mentioned unified preprocessing of spatial reference and time scale, the basic data acquisition module outputs the basic computational dataset aligned with the underlying data reference plane.

[0037] See attached document Figure 4 This invention provides a method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion, comprising the following steps: S21, in this embodiment, the topology weighting module calls the basic computational dataset that has completed spatiotemporal benchmark alignment and reads the internal digital elevation model to perform terrain smoothing and depression filling operations. Due to noise interference from surface vegetation cover and building shadows during airborne lidar scanning and mapping, local closed elevation depressions exist on the natural terrain surface. Such local closed elevation depressions can cause interruptions in hydrological calculations, leading to dead zones in the numerical model's pathfinding. To establish continuous surface flow paths, the topology weighting module identifies and locks the elevation minima within the digital elevation model grid. The elevation minima referred to here are isolated digital noise depressions where the elevations of all surrounding adjacent grids are higher than the central grid. The topology weighting module increases the elevation value of the corresponding elevation minima until the elevation value reaches the lowest overflow elevation value of the surrounding adjacent grid cells. In particular, to ensure physical accuracy, the filling and lifting of depressions are not performed on known river channels, natural lakes, ponds, and flood storage and detention areas designed in the engineering. Only isolated depressions caused by surveying noise are corrected to generate a corrected digital elevation base map without confluence interruption dead zones.

[0038] S22, the topology weighting module calculates the spatial water flow vector direction of surface grid cells based on the corrected digital elevation base map. The module uses a deterministic eight-direction flow convergence algorithm (D8) combined with a multiple flow direction algorithm and a minimum resistance path algorithm to assign weights, calculating the spatial elevation gradient from each central grid to its eight adjacent grids. This addresses the topology interruption caused by the unreachability of lateral floodplain grids due to a single D8 directed graph. Elevation gradient calculation is implemented using three-dimensional spatial location information. The module extracts the direction of the adjacent grid with the largest spatial elevation gradient value and sets this direction as the target water flow output direction of the central grid. The module integrates the water flow output directions of the central grid to establish a global flow direction matrix. The formula for calculating the spatial elevation gradient is configured as follows: ; In the formula, Represents the central grid reaching the first Spatial elevation gradient results of adjacent grids; The terrain elevation value representing the center grid; Representing the The terrain elevation values ​​of each adjacent grid; Represents the geometric center of the central grid reaching the first The actual spatial distance between the geometric centers of adjacent grids. When the... When adjacent grids are in orthogonal positions Assign the value as the grid side length; when the first When adjacent grid cells are located on the diagonal. The value is assigned as the square root of two times the grid side length, and the grid side length is determined based on the globally uniform grid size of the basic computation dataset.

[0039] S23, the topology weighting module extracts the hydrodynamic connectivity network by combining the global flow direction matrix and the polygon with the largest water body coverage in historical floods. The module uses the global flow direction matrix to deduce the water transfer process from upstream to downstream grids, accumulating the total number of upstream nodes for each grid cell along the path, and outputting a global runoff accumulation matrix. The runoff network generated solely based on digital terrain is detached from the actual flood spread boundary constraints and cannot truly reflect the floodplain diversion state caused by water level rise during flood season; therefore, actual boundary control needs to be introduced. The topology weighting module converts the polygon with the largest water body coverage in historical floods into a binary space mask layer (this step is only used for initial field constraint extraction of the connectivity network and is never used as the outer boundary for subsequent model calculations).

[0040] The topology weighting module prunes the global flow accumulation matrix based on the binary space mask layer, eliminating redundant flow paths outside the binary space mask layer. This pruning operation retains the flow accumulation values ​​of the grids within the binary space mask layer, while clearing the flow accumulation values ​​of the grids in the outer non-submerged areas to zero, thus preventing waterless areas from participating in subsequent topology construction. The topology weighting module sets a connectivity threshold, extracting and retaining grid node connections located within the binary space mask layer with flow accumulation values ​​greater than the connectivity threshold, generating a surface water flow connectivity network that conforms to real hydrodynamic characteristics. The connectivity threshold is set according to the number of grids covered by the actual average width of the upstream river channel in the study area, typically set to an integer range of 50 to 200, used to filter out minor surface capillary noise and ensure that the extracted surface water flow connectivity network spatially matches the actual physical channel development characteristics.

[0041] In this embodiment, the topology weight construction module performs the hydrological connectivity distance calculation model establishment operation, which includes the following steps: In step S24, the topology weighting module reads the spatial coordinates of discrete hydrological observation stations and performs node projection matching operations in conjunction with the surface water flow connectivity network. Due to coordinate measurement errors and grid discretization, the original coordinates of hydrological observation stations are often difficult to directly coincide with network nodes, and forced topology calculations can lead to pathfinding dead zones. The topology weighting module establishes a preset buffer zone around the hydrological observation station, with the search radius set according to the grid diagonal size, typically 1.0 to 2.0 times the diagonal length. If no node is found within the buffer zone, the search radius is expanded outward in steps of a single grid width until a node is captured, avoiding matching interruptions. The topology weighting module searches for network nodes within the buffer zone, calculates the straight-line distance between coordinates, and marks the network node with the smallest straight-line distance as an effective mapping anchor point, thereby reducing the probability of pathfinding interruptions caused by the point deviating from the water flow network.

[0042] S25, the topology weighting module establishes a connectivity distance calculation model to measure the path length of the spatial grid along the resistance. Flow on the surface of floodplains is constrained by topographic relief; traditional Euclidean distance does not consider the obstruction of dikes and the guiding effect of river channels, leading to misjudgments of hydraulic connections. Hydrological connectivity distance can more realistically represent the trajectory of water flow around obstacles. The topology weighting module sets the effective mapping anchor point as the pathfinding endpoint and the spatial grid within the mask as the pathfinding starting point, performing graph theory pathfinding operations along the connected network. Since the directed graph has been converted into an undirected weighted graph in the previous stage, this pathfinding process allows for real bidirectional distance detection of water flow nodes in the reverse direction or in laterally flat floodplains. The topology weighting module configures the connections between adjacent grids as graph structure edges, their lengths as basic weights, extracts the shortest topological path, and accumulates the basic weights to generate the hydrological connectivity distance. The calculation formula is configured as follows: ; In the formula, Representing the The spatial computing grid reaches the [number]th [number]. The hydrological connectivity distance of each effective mapping anchor point; This represents the total number of graph edges on the shortest topological path. Representing the The spatial physical length of the edges of a graph structure.

[0043] S26, the topology weighting module performs connectivity distance numerical compensation correction for isolated grid cells. Due to the influence of micro-topography, some shallow-water-inundated grids exist within the mask layer that are not directly associated with the network. Ignoring these would result in holes in the interpolation field. For isolated grid cells, the topology weighting module locates the nearest connected network access node, calculates the local straight-line distance, and adds it to the already calculated hydrological connectivity distance of that access node to form the compensated hydrological connectivity distance. The weight coefficient for the local straight-line distance is set to 1.0. This operation ensures that all associated spatial computational grids within the maximum inundation range can obtain continuous distance scales, avoiding numerical blanking anomalies during interpolation assimilation.

[0044] In this embodiment, the topology weight construction module performs the energy decay constraint spatial distribution weight generation operation, which includes the following steps: S27, the topology weighting module introduces an energy decay coefficient, transforming hydrological connectivity distance into an initial weight for the degree of impact on spatial nodes. When water flows across a complex bed surface, hydrodynamic energy dissipates nonlinearly with the accumulation of frictional resistance. Conventional spatial interpolation uses the reciprocal of the straight-line distance as the weight, neglecting the actual energy loss mechanism of the physical flow. The topology weighting module constructs a distance decay function that considers both topographic obstacles and energy dissipation, using this function to quantify the water level control effect of each effective mapping anchor point on a specific spatial computational grid. The calculation formula for the distance decay function is configured as follows: ; In the formula, Representing the The effective mapping anchor point pair of the th The initial weights of each spatial computation grid; Representing the The spatial computing grid reaches the [number]th [number]. The hydrological connectivity distance of each effective mapping anchor point; Represents an energy decay index; This represents the local minimum smoothing constant.

[0045] S28, the topology weighting module performs value calibration and anomaly prevention for the control parameters within the distance decay function. This is to avoid hydrological connectivity distance issues arising when the spatial computation grid and the spatial location of the effective mapping anchor point coincide in a discretized grid system. The denominator is 0, which triggers a program crash dead zone due to the denominator being 0. The topology weight construction module smooths the minimum value using a constant. Set as a constant greater than 0, with a typical value range of 0.001 to 0.01. Energy decay index. This is used to characterize the rate at which surface roughness in floodplains affects the energy consumption of water flow. The topology weighting module is determined based on the substrate roughness distribution characteristics of the study area provided by the basic computational dataset. The value of this parameter is typically set between 1.5 and 3.0. When the surface vegetation is dense or the water-blocking effect of buildings is strong, a value closer to the upper limit is selected to accelerate the decay of the initial weight with distance, thus better reflecting the physical law of water surface line drop under the actual roughness. It should be clarified here that this parameter... As a distance decay parameter used in the initial field space assimilation stage, it no longer interferes with or participates in the control and solution process of the hydrodynamic physical equations after the model completes parameter calibration and enters the forward hydrodynamic solution.

[0046] In step S29, the topology weight construction module performs global normalization on the initial weights to generate the final spatially distributed weight matrix. During the spatial interpolation assimilation stage, if weight values ​​are not normalized, the interpolated intermediate water level values ​​may exceed the boundary physical constraints of the measured stations, violating the principle of mass conservation in fluid dynamics. Furthermore, forcing all stations in the region to participate in the calculation of a single grid will cause matrix dimensionality expansion and introduce noise distortion from distant, weakly correlated stations.

[0047] The topology weighting module sets a local interpolation search radius and extracts only the set of valid mapping anchor points within this radius. The search radius is defined based on the average width of the historical flood inundation cross-section of the target watershed. The module accumulates the initial weight values ​​corresponding to each valid mapping anchor point within the search radius to generate the total weight of the spatial computation grid. The module then divides the initial weight of a single valid mapping anchor point by this total weight to derive a standardized weight value that sums to a constant 1. The corresponding normalized calculation formula is expressed as: ; In the formula, This represents the normalized weight values; This represents the total number of valid mapping anchor points extracted within the search radius. Representing the The initial weights of each effective mapping anchor point for the spatial computation grid.

[0048] After normalization, the topology weighting module integrates the standardized weight values ​​of the global spatial computational grid relative to each effective mapping anchor point, outputting a spatial distribution weight matrix that encompasses spatial connectivity and energy attenuation characteristics. This spatial distribution weight matrix is ​​used for subsequent underlying spatial assimilation operations to derive the global continuous spatial water level distribution values.

[0049] In this embodiment, the initial field generation module performs the boundary mask function space mapping setting operation, which includes the following steps: In S31, the initial field generation module loads a globally unified computational grid system aligned with the underlying data baseline and simultaneously extracts the polygon with the largest historical flood water coverage as a geographical constraint boundary. Forcibly performing spatial assimilation interpolation on the entire grid would cause matrix operation dimensionality expansion and excessive memory consumption, and could easily generate false water levels in non-floodplain areas, leading to numerical divergence dead zones in the solver due to excessively large local water surface elevation gradients. Therefore, the initial field generation module converts this polygon into a control boundary (this polygon serves only as a spatial constraint mechanism for historical verification and interpolation, and is never used as a forced deadlock boundary for subsequent forward simulations). It then performs spatial overlay and intersection analysis between the polygon vector layer and the rasterized globally unified computational grid system, using the actual inundation range to define the effective computational area.

[0050] S32, the initial field generation module establishes a two-dimensional spatial boundary mask function, quantifies the spatial affiliation attributes of each computational grid, and delineates the wet and dry boundaries of the computational domain. The initial field generation module extracts the coordinates of the spatial geometric center point of the grid and uses a ray crossing algorithm to determine whether it falls within the two-dimensional closed space formed by the aforementioned polygons. To avoid discretization truncation errors causing edge-crossing grids to be missed, when the geometric center point of a grid is located outside the polygon, but its overlap area with the interior of the polygon reaches a preset area threshold, the initial field generation module performs forced merging of the out-of-bounds grid. For the overlap ratio calculation, the initial field generation module uses a grid subdivision algorithm to divide the target grid into a sub-pixel array and approximate the solution by counting the number of sub-pixels falling inside the polygon. The preset area threshold is set to 0.4–0.6 based on the grid size. The mathematical expression is configured as follows: ; In the formula, Representing the Line 1 The output of the two-dimensional spatial boundary mask function of the computational grid is 1, which means it belongs to the valid computational region and 0 means it belongs to the invalid masking region. Represents the coordinates of the spatial geometric center point of the corresponding grid; Represents a two-dimensional closed space composed of polygons; This represents the percentage of overlapping area. This represents the preset area threshold.

[0051] In step S33, the initial field generation module performs a spatial mapping isolation operation based on the results of the two-dimensional spatial boundary mask function. This result is transformed into a spatial mapping Boolean matrix of the same dimension as the computational grid, and this matrix is ​​used to perform dot product filtering on the grid system. The module removes invalid masked regions with an output of 0, preventing non-aqueous grids from obtaining interpolated water levels, and assigns a no-data flag (e.g., -9999) to the corresponding matrix node values ​​to prevent non-numeric error dead zones in matrix operations. The module extracts valid computational regions with an output of 1 to generate a discretized computational domain, removing redundant terrain background to ensure that subsequent numerical assimilation operations and hydrodynamic solutions have spatial constraints.

[0052] In this embodiment, the initial field generation module performs the spatiotemporal dynamic assimilation operation of discrete hydrological sequences, which includes the following steps: S34, the initial field generation module initiates the time-series extrapolation process based on the discretized sequence of time steps. The hydrodynamic model extrapolation relies on the initial flow field input, which dynamically evolves over time. The initial field generation module extracts the measured water level values ​​corresponding to each effective mapping anchor point at the current time step node. Sensor communication delays or failures can easily lead to data loss, resulting in assimilation interruption dead zones. If relying solely on linear interpolation, algorithm anomalies may occur when sequence endpoints are missing or when consecutive long time intervals are lost. If the current time step lacks corresponding measured data, the initial field generation module detects the temporal position of the missing data. When the missing data is located in the middle of the sequence and the consecutive missing span is less than the preset step-size threshold, it is filled using measured data from adjacent time points through a linear interpolation mechanism; when the missing data is located at the end of the sequence, it is filled using the effective measured values ​​from the nearest neighbor time point; when the consecutive missing span is greater than or equal to the preset step-size threshold, the historical average of observations from the same period is introduced for replacement. The preset step-size threshold is set based on the stability of hydrodynamic evolution and is typically configured as 3 to 5 time steps. Through the above-mentioned multi-level fault tolerance mechanism, a basic hydrological driving sequence with continuous time dimension is formed.

[0053] S35, the initial field generation module combines the spatially distributed weight matrix and the basic hydrological driving sequence to generate a spatially continuous water level distribution field. In physical hydrology, the water level state at any location within a floodplain is jointly controlled by multiple surrounding inflow nodes, and the control intensity decays along the flow energy. For each computational grid within the discretized computational domain, the initial field generation module multiplies the standardized weight values ​​of the effective mapping anchor points associated with the computational grid with the corresponding measured water level values, performs an accumulation operation on the product results, and derives the initial assimilated water level value of the target computational grid at the current time step. The mathematical expression for the spatially continuous assimilation operation is configured as follows: ; In the formula, Representing the The computational grid at the _th ... Initial assimilation water level values ​​at each time step; This represents the total number of valid mapping anchor points associated within the local interpolation search radius. Representing the Each effective mapping anchor point is compared to the normalized weight value of the computational grid; Representing the The effective mapping anchor point is at the _th _th The measured water level values ​​at each time step.

[0054] In S36, the initial field generation module incorporates topographic elevation data to perform physical bottom-edge constraint verification on the initial assimilated water level values. Since spatial assimilation interpolation is essentially a smooth surface fitting within a purely mathematical dimension, it lacks the physical perception of underlying surface undulations. The generated water level is prone to being lower than the surface elevation in local topographic uplift areas. Directly inputting such data into the hydrodynamic solver will produce negative water depths, leading to divergence errors and dead zones in the partial differential equation calculations. The initial field generation module extracts the topographic elevation corresponding to each computational grid from the corrected digital elevation map and compares the initial assimilated water level values ​​with the corresponding topographic elevations. When the initial assimilated water level value is lower than the topographic elevation, the initial field generation module forcibly corrects the water level value of that computational grid to the topographic elevation superimposed with a very thin water layer constant. This extremely thin water layer constant is only used as a transitional processing of numerical stability for the dry-wet boundary computational unit; for inland grids far from the current wet boundary and not yet meeting the flooding criteria, the model maintains its pure dry unit state, and is dynamically activated through a rigorous dry-wet conversion mechanism during subsequent physical evolution calculations to avoid the non-physical anomaly of the screen being filled with water film.

[0055] In the numerical solution domain of the shallow water equation, the extremely thin water layer constant is typically configured to be 0.001–0.005 m to maintain the continuity of the spatial flow field calculation while avoiding negative water depth anomalies. If the initial assimilated water level is greater than or equal to the topographic elevation, the initial field generation module keeps the initial assimilated water level unchanged.

[0056] After physical bottom constraint verification, the initial field generation module traverses all time steps within the preset time period in the order of time evolution, and saves the global grid water level corresponding to each time step as a multi-dimensional spatiotemporal matrix to generate a discrete hydrological sequence spatiotemporal dynamic assimilation field that supports the evolution of flood dynamics in the floodplain.

[0057] See attached document Figure 5 In this embodiment, the initial field generation module performs the physical constraint initial field boundary condition transformation operation, which includes the following steps: S37, the initial field generation module extracts the water level surface within the spatiotemporal dynamic assimilation field of the discrete hydrological sequence and generates an effective water depth distribution matrix by combining it with a corrected digital elevation base map. In solving the two-dimensional shallow water dynamic equations, directly using absolute water level elevation as the driving variable is susceptible to interference from datum differences and typically needs to be converted to relative water depth. For each computational grid within the discretized computational domain, the initial field generation module subtracts the corresponding topographic elevation from the initial assimilated water level value to obtain the water depth value at the current time step. To prevent abrupt water depth differences between adjacent grids, which could lead to numerical oscillations or even divergence dead zones in the solver's internal flow convection term calculations, the initial field generation module calculates the water depth gradient between adjacent grids. When a local water depth gradient exceeds a preset gradient threshold, the initial field generation module performs a mean smoothing filter operation on that local grid. Specifically, the initial field generation module constructs a 3×3 spatial computation window containing the target computational grid and its eight neighboring grids. It extracts the water depth values ​​from the effective computational grids within this window, calculates their arithmetic mean, and replaces the original target computational grid's water depth value, thereby mitigating computational instability caused by local abrupt changes. The preset gradient threshold is set according to the grid size and terrain undulation, typically ranging from 0.1 to 0.5. The formula for calculating the effective water depth value is configured as follows: ; In the formula, Representing the The computational grid at the _th ... Effective water depth values ​​at each time step; Representing the The computational grid at the _th ... Initial assimilation water level values ​​at each time step; Representing the The terrain elevation corresponding to each computational grid.

[0058] S38, the initial field generation module sets the boundary conditions of the initial momentum field based on the effective water depth distribution matrix. Flood dynamics simulation requires scalar-dimensional water depth and vector-dimensional flow velocity as the starting state for the time integration of the equations. Since conventional hydrological observations cannot provide a high-precision continuous velocity field across the entire region, manually specifying the initial velocity can easily lead to shock wave effects in the initial stage due to the non-conservation of mass and momentum between grids. The initial field generation module employs a cold-start mechanism, forcibly setting the initial velocity components of the entire grid to zero. This ensures that the water flow dynamics in the initial stage are mainly driven by the gravitational component generated by the water level gradient, gradually approximating the real flow field through the natural conversion of gravitational potential energy into kinetic energy. Meanwhile, the initial field generation module sets a no-slip solid boundary condition on the outer perimeter of the polygon with the largest water body coverage in historical floods. That is, both the normal flow velocity and the tangential flow velocity are configured to be 0, which restricts the abnormal diffusion of water into the ineffective shielding area. This no-slip setting is only used as an initial stability protection at the start time (assimilation stage). The subsequent forward simulation process will be automatically promoted based on the dynamic evolution mechanism of the global dry and wet boundary, and will no longer be controlled by this rigid dead zone boundary.

[0059] In step S39, the initial field generation module converts the physically constrained water level surface, effective water depth distribution matrix, and initial momentum field into a format that encapsulates them into a physically constrained initial field dataset. Different underlying hydrodynamic solvers have specific differences in input file formats. To ensure compatibility with standard computational frameworks, the initial field generation module maps the aforementioned multidimensional matrices to a standard network-wide data format. During the encapsulation process, the initial field generation module writes timestamps and grid projection coordinate system information into the dataset to ensure that the output driver file is aligned with the hydrodynamic solver in both the time series and spatial positioning dimensions, avoiding program dead zones caused by misalignment of spatiotemporal references.

[0060] The encapsulated physical constraint initial field dataset is directly used as the underlying initial calculation condition for subsequent flood dynamic replay calculations in the floodplain (without making dynamic boundary inputs that cover every time step), in order to maintain the continuity of data flow in the overall technical solution.

[0061] In this embodiment, the model building and verification module performs the following steps to build a one- or two-dimensional coupled shallow water dynamics calculation framework: S41, the model building and verification module establishes the topological interaction between the one-dimensional channel and the two-dimensional floodplain. In real flood evolution, frequent mass and momentum exchanges occur between the main channel and the floodplains on both sides. The module extracts the main channel boundary line, arranges lateral connecting weirs along this boundary line, and spatially binds the one-dimensional channel nodes to the two-dimensional computational grid. Based on the explicit finite volume method, the module calculates the flow transfer flux between the grids and the interface using the Roe or HLLC scheme. For flow calculation, the module uses a broad-crested weir flow mechanism to extrapolate the interactive flow using the one-dimensional and two-dimensional water level differences. To prevent abnormally large increases in the exchange flow due to excessive topographic elevation differences, which could lead to negative water depths and computational divergence dead zones, the module configures an interactive flow limiter. When the water level difference between nodes exceeds a preset threshold, the module uses a critical flow mode to limit the maximum overflow rate. The preset water level difference threshold is set based on the main channel depth, with a range of 0.5–2.0 m.

[0062] S42, the model building and verification module configures a two-dimensional shallow water dynamic equation set including source terms, establishing a mass and momentum conservation system. The two-dimensional continuity governing equations are used to calculate the mass conservation state of the flow field, reflecting the water recharge process from one-dimensional overflow to the two-dimensional floodplain. The model building and verification module uses the exchange flow of the lateral connecting weir as the source term input into the bound two-dimensional grid. The mathematical expression of the two-dimensional continuity governing equations is configured as follows: ; In the formula, The effective water depth representing the two-dimensional spatial computational grid (compared to the aforementioned discrete state) equivalence); The continuous fluid dynamic time represented by the model derivation; Representative at Flow rate per unit width in the direction; Representative at Flow rate per unit width in the direction; This represents the unit area interactive source term flow transported from a one-dimensional channel to a two-dimensional spatial computational grid.

[0063] Meanwhile, the model construction and verification module incorporates the momentum conservation equation. The momentum conservation calculation explicitly includes the coupled expression of the convection term transmission mechanism, the gravity component driven by the water level gradient, and the Manning resistance friction term of the substrate. This is used to verify the physical constraints of water flow convection, gravity gradient, and substrate friction on the change of water momentum, thus constructing the underlying dynamic calculation architecture.

[0064] S43, the model building and verification module sets a dynamic time step constraint mechanism to maintain the stability of the numerical solution. In the explicit finite volume method, if the time step is too large, the physical wave velocity of the water flow will cross multiple grids within a single step, causing the Courant number to exceed the limit and numerical oscillation dead zone. The model building and verification module introduces the Courant-Friedrich-Lévy condition to perform step updates. In shallow water systems, the disturbance propagation velocity is the sum of the current velocity and the gravity wave velocity. The model building and verification module extracts the maximum current velocity across the entire grid, calculates the local gravity wave velocity by combining gravitational acceleration and water depth, and uses the minimum grid size to deduce the allowable time step that satisfies stability. The Courant number safety threshold range is configured to 0.5–0.9 to avoid solver crashes while maintaining computational efficiency.

[0065] See attached document Figure 6 In this embodiment, the model building and verification module performs the spatial functional difference multi-objective evaluation function setting operation, which includes the following steps: S44, the model construction and verification module performs spatial functional zoning of the computational grid based on the underlying land use attributes. Traditional hydrodynamic model parameter verification often uses globally unified error evaluation standards, rarely considering the differences in the sensitivity of different geographic spaces to flood inundation. The model construction and verification module extracts land cover type vector data of the study area and divides the discretized computational domain into the main channel floodplain, densely populated urban residential area, and agricultural ecological floodplain. Specifically, the model construction and verification module uses a spatial overlay cross-analysis algorithm to compare the position of the geometric center point of the discrete computational grid with the vector layer. The model construction and verification module assigns independent spatial identifiers to each spatial functional zone and establishes a mapping matrix between the computational grid and functional attributes.

[0066] S45, the model building and verification module sets differentiated basic error assessment indicators for different spatial functional areas. Densely populated urban areas have low tolerance for excessively high water levels. The module compares measured peak flood levels from hydrological stations with simulated levels, setting the root mean square error of peak flood levels as a core indicator to strictly constrain the accuracy of water surface line fitting. Agricultural ecological floodplains focus on assessing flood retention time. The module compares the actual inundation period monitored by remote sensing with the simulated receding time, setting the relative error of inundation duration as an assessment indicator. The main channel floodplain combines measured flow sequences to configure the Nash efficiency coefficient of the flow process line, assessing the dynamic trend of water transport. Through the configuration of differentiated indicators, the core of the assessment of overall physical connectivity and the fidelity of water spread restoration is the spatial intersection-combination ratio (SCR) between the remotely sensed inundation range and the simulated inundation range. The core evaluation indicators and their calculation logic are as follows: ; In the formula, Represents the intersection and union ratio of spatial representations; The simulated flooded area represents the output of the model; This represents the measured water area as interpreted from remote sensing images.

[0067] The main channel floodplain uses a Nash efficiency coefficient for the flow process line, based on measured flow sequences, to assess the dynamic trend of water transport. By configuring differentiated indicators, the model evaluation dead zone, where global indicators mask local fitting biases, is prevented.

[0068] S46, the model building and verification module constructs a multi-objective evaluation function through a weight combination mechanism. Due to the differences in the physical dimensions of the various indicators, direct summation can lead to a computational distortion dead zone where large numerical indicators mask smaller numerical indicators. The model building and verification module extracts the historical statistical extreme values ​​of each indicator and performs range normalization. To prevent a program error dead zone caused by the maximum and minimum values ​​being equal (resulting in a denominator of 0), when the range value is less than a preset small tolerance (usually configured as 1×10⁻⁶), a range normalization process is implemented. 6 When this occurs, the model building and verification module forcibly assigns a default normalized baseline value to the indicator. Subsequently, the model building and verification module introduces the analytic hierarchy process (AHP) to assign corresponding weights based on the importance of the flood protection objects, with the weight constant for densely populated urban residential areas configured to be 0.5–0.7. The mathematical expression of the multi-objective evaluation function is configured as follows: ; In the formula, The output value of the multi-objective evaluation function representing the overall model; This represents the total number of functional zones in the space. Representing the The weight constants configured for each spatial functional area are constrained to a constant 1. Representing the The basic error assessment index values ​​of each spatial functional area after normalization (the index values ​​include the aforementioned root mean square error of water level, Nash efficiency coefficient, and spatial IoU overlap, etc.).

[0069] The model building and verification module uses the output values ​​of the multi-objective evaluation function as the optimization benchmark for automatic calibration of the model's underlying parameters, ensuring that the calculation and solution prioritize the flood control needs of highly sensitive areas.

[0070] In this embodiment, the model building and verification module performs the reverse optimization iterative parameter calibration operation, which includes the following steps: S47, the model building and verification module establishes the space of parameters to be determined and configures the constraint boundaries. Parameter calibration is essentially an inversion process that maps multi-objective simulation errors to spatial physical baseline characteristics. In hydrodynamic principles, the resistance and overflow efficiency of a hydrodynamic model depend on physical parameters, and traditional manual parameter tuning based on experience is usually insufficient to achieve the optimal solution for multiple objectives. The model building and verification module extracts the Manning roughness coefficient and the lateral weir flow coefficient of each spatial functional zone within the discrete computational domain and assembles them into a vector of parameters to be optimized. If the algorithm generates resistance parameters less than or equal to 0, it will trigger an error dead zone where the denominator of the friction slope is 0. Therefore, the model building and verification module configures physical boundaries based on vegetation attributes. The Manning roughness coefficient is configured to range from 0.015 to 0.120, and the lateral weir flow coefficient is configured to range from 0.30 to 0.85. The model building and verification module uses Latin hypercube sampling to generate the initial population, ensuring the uniformity of parameter distribution within the multidimensional search space.

[0071] In step S48, the model building and verification module initiates a swarm intelligence reverse optimization loop based on the particle swarm optimization algorithm, calling the hydrodynamic framework to perform trial calculations. The module substitutes the initial population into the framework to extrapolate the hydrological sequence and calculates the fitness value using a multi-objective evaluation function. During the iterative solution process, unreasonable parameter combinations can easily lead to Corange number exceeding the limit, causing the solver to diverge, crash, and interrupt the process in a dead zone. The model building and verification module incorporates an anomaly capture mechanism; when solver divergence is detected, a large penalty value is assigned to the solution (e.g., configured as 1×10). 10 This forces the algorithm to discard dangerous regions. The optimization speed evolution formula is configured as follows: ; In the formula, and Representative individuals in the first Next and first The speed of each iteration; This represents the inertia weight in the particle swarm optimization algorithm, and its value range is configured to be 0.4 to 0.9. and The cognitive and social learning factor is typically configured as a constant of 2.0. and Represents a random value that is uniformly distributed within the interval 0 to 1; and Representing up to the Individual and global optimal solutions; This represents the current location.

[0072] S49, the model building and verification module performs convergence verification, extracts the optimal solution, and calibrates it. To avoid getting trapped in local optima and creating an infinite loop (a control dead zone), the model building and verification module monitors changes in the fitness of the optimal solution. When the absolute value of the change in consecutive iterations is less than a preset tolerance (usually configured as 1×10⁻⁶), the module will perform a calibration. 5 The loop terminates when the cumulative iterations reach the maximum iteration threshold (configured to 100-200 generations). The model building and verification module extracts the optimal solution vector, deconstructs it into functional area physical parameters, and writes it back to the underlying configuration file to complete the calibration.

[0073] In this embodiment, the dynamic simulation evaluation module performs multi-scene terrain configuration and forward numerical simulation operations, including the following steps: S51, the dynamic simulation evaluation module configures multi-scenario geomorphic feature files based on preset model evaluation requirements. In surface hydrology, geomorphic changes directly alter the path of flood spread. The dynamic simulation evaluation module extracts digital elevation model files and topographic change data under different engineering scenarios, uses a spatial location overlay algorithm to locate the computational grid affected by engineering intervention, and simultaneously updates the roughness partition map, physical boundaries of water-blocking entities, and connectivity exchange parameters of lateral connecting weirs for the corresponding areas according to engineering boundaries. Subsequently, it introduces scenario indication matrices and engineering design elevations to quantitatively correct the topography. The mathematical expression for multi-scenario geomorphic configuration is configured as follows: ; In the formula, Representing the The computational grid at the _th ... Corrected terrain elevation for each scenario; This represents a spatial location indicator variable, taking the value 1 when the grid is inside the intervention area and 0 when it is outside. Represents the original terrain elevation; Represents the engineering design elevation.

[0074] To prevent the creation of false depressions and the resulting local catchment dead zones in the flow field calculation, the dynamic simulation evaluation module performs a rationality check on the corrected terrain elevation. The dynamic simulation evaluation module scans the eight neighboring elevations of the target grid. When the corrected terrain elevation of the target grid is lower than the lowest elevation of all surrounding neighboring grids, the depression filling logic is triggered, forcibly raising the elevation of that grid to the lowest neighboring elevation value to maintain the gravity continuity of the flow calculation.

[0075] S52, the dynamic simulation evaluation module loads the corresponding boundary condition driving sequence. Based on the multi-purpose flood return period standard and reenactment requirements, the historical measured process line is called when used for verification of real historical events, and the flood return period design flow process line is used when evaluating the capacity of the preceding engineering.

[0076] The dynamic simulation evaluation module inputs the flow process curve into the one-dimensional channel node. Simultaneously, it retrieves the initial physical constraint dataset as the starting state for the numerical simulation. To address temporal misalignment anomalies caused by asynchronous observation times, the dynamic simulation evaluation module uses the initial field timestamp as a reference to perform time-axis alignment interpolation on the driving sequence, ensuring data continuity and avoiding dead zones caused by abrupt boundary changes that could lead to model read errors.

[0077] S53, the dynamic simulation evaluation module invokes a one- or two-dimensional coupled shallow water dynamics calculation framework to perform forward numerical simulations in parallel for each evaluation scenario. The numerical solver, based on corrected topographic elevations and driving sequences, solves the water depth and velocity across the entire grid within the computational domain of the potential wide floodplain, jointly determined by the watershed outline and the potential design flood impact surface. It does not lock the boundary mask, allowing for the evolution of wet and dry states across the entire domain, and recursively solves for the water depth and velocity across the entire grid step-by-step. To prevent sudden strong source inputs from causing a sharp increase in local velocity, which could lead to Corange number exceedances and solver divergence dead zones, the dynamic simulation evaluation module enables an adaptive step-size fine-tuning procedure. When a local Corange number exceeds the stability critical value, the dynamic simulation evaluation module automatically halves the computational time step. To prevent the program from entering an infinite loop, the dynamic simulation evaluation module configures a minimum step-size threshold (usually 0.01s). If the halved time step is lower than this threshold, the halving operation is stopped, and a numerical artificial viscosity term is introduced into the locally exceeding the boundary grid. This forcefully suppresses numerical oscillations through a kinetic energy dissipation mechanism, and then the time integration continues.

[0078] In this embodiment, the dynamic simulation evaluation module performs the quantitative comparative evaluation and analysis of flood storage and detention functions, including the following steps: S54, the dynamic simulation evaluation module extracts the spatiotemporal evolution flow field sequence from the positive numerical simulation output and identifies characteristic physical quantities within the floodplain. This step rigorously quantifies key parameters of hydrodynamic analysis: extracting the absolute time difference between the moment when the water depth of a single grid first crosses the effective water depth threshold and the moment when it finally falls below that threshold during the receding period, which is determined as the flood retention time; at preset key flood channels or control sections, extracting the percentage of the cumulative flow over the entire period relative to the total inflow of the system, and calculating the diversion ratio of the generated channels.

[0079] For each assessment scenario, the dynamic simulation assessment module traverses the spatial computational grid within the computational domain to extract the maximum inundation depth and extreme flow velocity. Shallow water dynamics often introduces small numerical residual depths due to truncation errors when solving for wet-dry boundaries. To prevent these residuals from causing algorithmic dead zones that distort the determination of effective flood storage area, the dynamic simulation assessment module configures an effective depth threshold (0.05–0.1 m) to filter out grids below this threshold. Simultaneously, localized low-lying areas are susceptible to isolated water accumulation due to minor overflows; including these areas would lead to inflated volumetric indices. To prevent false catchment dead zones, the dynamic simulation assessment module introduces a hydraulic connectivity topology search mechanism. Starting from the one-dimensional main channel overflow point, it performs an eight-neighbor connectivity traversal outwards, eliminating isolated grids not physically connected to the main channel to determine the actual set of effective flood storage grids participating in peak flow reduction.

[0080] S55, the dynamic simulation evaluation module calculates the actual total effective flood storage and detention volume under various scenarios based on the effective flood storage and detention grid set. To avoid the overestimation error of physical conservation caused by the overlay of the envelope of a traditional single grid throughout the entire time period, the dynamic simulation evaluation module implements the calculation with a rigorous global time-domain integral optimization mechanism: first, it independently calculates the instantaneous total effective flood storage and detention volume at each time step node within the simulation period, and then extracts the highest peak value from the time-series volume array as the final quantification result. That is: it combines the grid spatial area to perform the summation and integration of the instantaneous area product within the time step, and performs water level elevation truncation verification to deduct the ineffective dry land volume within the grid that is higher than the surrounding water surface. The mathematical expression for the total effective flood storage and detention volume is strictly configured as follows: ; In the formula, This represents the total effective flood storage and detention volume at the transient peak value under a specific assessment scenario. Represents the entire continuous simulation time domain; This represents the total number of computational grids within the effective flood storage and detention grid set; Representing the A specific effective flood storage and detention grid in a given... Instantaneous effective water depth at a given time step; Representing the The spatial planar area of ​​an effective flood storage and detention grid.

[0081] S56, the dynamic simulation evaluation module performs multi-scenario spatial flow field comparison to generate flood storage and detention function evaluation indicators. Using the original natural landform as the baseline scenario, the module calculates the flood characteristic attenuation rate under engineering intervention scenarios, encapsulating the peak flood level reduction, the increase in effective flood storage and detention volume, and the duration of flood spread into a comparative evaluation matrix. To prevent numerical distortion dead zones caused by direct comparison of heterogeneous physical dimensions, the module employs a range normalization method to perform dimensionless processing on the matrix data. Subsequently, the module introduces the entropy weighting method, calculating objective weights based on the data information entropy of each indicator (where the weight constant for peak flood level reduction is typically configured to 0.4–0.6), and then performs weighted calculations to output a quantified comprehensive flood control benefit score.

[0082] The comprehensive score and comparison matrix output by the dynamic simulation evaluation module are written into the evaluation report, providing underlying physical quantitative data support for the optimization of the layout of watershed flood control projects.

[0083] Specific application example: Dynamic simulation and quantitative evaluation of comprehensive benefits of flood storage and detention in a typical watershed flood control project. To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments and related drawings. It should be understood that the specific embodiments described herein are merely illustrative and are not intended to limit the scope of protection of this invention.

[0084] This embodiment takes a typical flood-prone watershed (referred to as the M watershed) as the research object, extracts a historical real rainfall and flood event with a 50-year return period (lasting 72 hours) for dynamic reenactment, and evaluates the flood control benefits of the proposed flood storage and detention area based on this, aiming to illustrate the specific calculation process and implementation effect of the method provided in this application.

[0085] Basic data processing and spatial topology weight assimilation: When constructing a watershed hydrodynamic model, spatial interpolation is required because discrete hydrological station data cannot directly characterize the spatially continuous flow state of the entire watershed.

[0086] The basic data acquisition module collects spatial topographic data of the M watershed and resamples it into a uniform 30m×30m computational grid. In the initial field generation stage, in order to deduce the spatially continuous water level distribution field, the topology weight construction module calculates the hydrological connectivity distance and spatial distribution weight of discrete hydrological stations.

[0087] Numerical computation example: Assume a target space computational grid within a discretized computational domain. There are two valid mapping anchor points (hydrological station A and hydrological station B) within the local interpolation search radius.

[0088] Calculations show that the grid... Hydrological connectivity distance to hydrological station A =300m, hydrological connection distance to hydrological station B =600m. Set the energy decay index. =2.0, (this index is usually empirically determined based on the watershed topographic relief between 1.5 and 3.0; in this embodiment, 2.0 is used to characterize moderately undulating terrain), minimum smoothing constant. =0.001 (used to prevent calculation overflow caused by a denominator of 0).

[0089] According to the initial weight formula The calculation shows that: ≈1.11×10 -5 , ≈2.78×10 -6 .

[0090] Using the normalization calculation formula The standardized weight values ​​are obtained as follows: =0.8, =0.2.

[0091] In the At a given time step, if the measured water level =25.0m, =24.0m, then the initial assimilation water level value =0.8×25.0+0.2×24.0=24.8m. This value has been verified by the physical bottom constraint of the terrain elevation (i.e., to determine if the calculated water level...). If the surface elevation is lower than that of the grid, the water level of that grid will be forcibly assigned the surface elevation to conform to the hydrophysical laws, and then the initial momentum field will be applied.

[0092] Multi-objective evaluation function configuration and reverse optimization: Due to the complex surface cover types within the watershed, a single roughness coefficient is insufficient to accurately reflect flood evolution resistance. Therefore, after loading the initial field, the model construction and verification module performs parameter calibration. The M watershed is divided into densely populated urban areas and agricultural ecological floodplains. Weights are assigned based on the differences in the importance of different areas in flood control safety. In this embodiment, the weight for the urban area is configured as follows: =0.6, the weight allocation for the floodplain area is =0.4.

[0093] For agricultural ecological floodplains, the measured inundation area is extracted using remote sensing imagery. =120ha, the simulated inundation area output by the model. =130ha, the spatial overlap area (intersection) of the two is 100ha, and the combined total area (union) is 150ha.

[0094] The core evaluation indicators are calculated based on the formula: =150 / 100=0.667 Will After normalizing the root mean square error of water levels in urban areas, the result is substituted into the multi-objective evaluation function. (in The total number of evaluation indicators, For the first The index deviation value after range normalization). This evaluation function... As the fitness function of the particle swarm optimization algorithm, when The algorithm stops when it converges continuously to the preset extreme value or reaches the maximum number of iterations (e.g., 100 times), thus completing the reverse optimization iteration of the Manning roughness coefficient (final calibration result: 0.025 for the main channel and 0.055 for the floodplain).

[0095] Experimental verification and effect comparison: To verify the effectiveness of the spatial assimilation and multi-objective verification method proposed in this invention in practical applications, a control experiment was conducted. Two simulation modes were set up: Traditional model: The inverse distance weighting method is used to interpolate the initial field, and a single uniform Manning roughness (0.040) is used for calibration using the traditional trial-and-error method across the entire domain.

[0096] The model of this invention: adopts the extraction of surface water flow connectivity network, setting of differentiated objective function and reverse optimization iteration.

[0097] like Figure 7 As shown, Figure 7 A line graph comparing the flow process curves at the main channel control section. Figure 7 The horizontal axis represents the flood evolution time (h), and the vertical axis represents the cross-sectional flow rate (m³). 3 / s). Those skilled in the art can intuitively distinguish different simulation results by the line types in the figure: the thick black solid line represents the measured hydrological sequence, the black dashed line represents the simulated value of the model of this invention, and the black dotted line represents the simulated value of the traditional model.

[0098] Combination Figure 7 The data trends show that during the initial rise phase (0-24 hours), the three curves largely overlap. However, during the peak flood period (30-45 hours), the traditional model fails to accurately capture the momentum decay caused by lateral overflow from the floodplain, resulting in a lower calculated peak flood level (approximately 1850m). 3 The flood peak ( / s) was significantly higher than the measured peak flood (1620m). 3The peak flow rate was 1645 m³ / s, and the peak time was advanced by approximately 2.5 hours. Conversely, the dashed line output by the model of this invention has a higher degree of agreement with the measured thick solid line, and more objectively reproduces the peak flow rate (1645 m³ / s). 3 / s (with an error of approximately 1.5%) and the complete receding tail process. Quantitative calculations show that the Nash efficiency coefficient of the flow process line configured in the main channel floodplain using the method of this invention reaches 0.92, which is a significant improvement over the 0.74 of the traditional model, indicating that this scheme has better applicability in terms of numerical simulation accuracy and dynamic process reflection.

[0099] Quantitative comparative evaluation of multi-scenario terrain configuration and flood storage and detention functions: After the model verification is passed, the dynamic simulation evaluation module evaluates the intervention effect of the proposed Project X (including flood control dikes and flood storage areas) in the watershed based on the design return period flood.

[0100] Two scenarios are set up: natural terrain (Scenario A) and post-construction terrain (Scenario B). The formula is used... Update grid elevation (where, The original elevation of the grid. To design the surface elevation of the project, For spatial Boolean variables, when the grid... (Take 1 if it is located in the engineering area, otherwise take 0), and run parallel forward numerical simulation.

[0101] The calculation process for flood storage and detention volume is as follows: Under scenario B, the dynamic simulation evaluation module identified an effective flood storage and detention grid set (a total of 5000 grids, with an area of ​​[missing information] per grid). =900m 2 ). In the The transient peak value was detected at time 40h, and the average effective water depth was measured. =2.5m.

[0102] The total effective flood storage volume at peak value is calculated using the integral optimization formula: =5000×2.5×900=11,250,000m 3 (In the formula, The total duration of the flood's evolution, The total effective flood storage and detention grid is calculated to be 11.25 million cubic meters.

[0103] To further quantify the differences in flood control benefits before and after the project implementation, please refer to [link / reference needed]. Figure 8 . Figure 8The chart shows the bar chart results of the quantitative comparative evaluation of flood storage and detention functions. The horizontal axis represents three types of core evaluation indicators, and the vertical axis reflects the relative distribution magnitude of each indicator. The actual evaluation and calculation values ​​are labeled at the top of each gray bar, with the light gray bar and the dark gray bar corresponding to the aforementioned scenario A and scenario B, respectively.

[0104] analyze Figure 8 The bar chart and numerical labels show that, in terms of the total effective flood storage and detention volume, the value marked in Scenario B (dark gray bar) reaches 11.25 million m³. 3 Compared to the natural overflow state of scenario A (light gray column, marked as 3.2 million m³), 3 The efficiency was increased by more than 2.5 times, a result consistent with the calculated values ​​using the aforementioned volume formula. Under engineering intervention, the peak flood level reduction indicator showed that the peak flood level at key downstream town sections was reduced by 1.25m from the uninterrupted state, causing it to fall below the warning safety level; the flood spread delay indicator showed that the time for floodwaters to reach highly sensitive residential areas was delayed by 4.5 hours. The dynamic simulation evaluation module ultimately used the entropy weight method (where the peak flood reduction weight was configured to 0.55 based on local flood control needs) to extract the characteristics of the above indicators and output a comprehensive flood control benefit score. Scenario B scored 89.5, higher than Scenario A's 42.1. This quantitative comparison data can serve as an objective reference for the feasibility analysis of flood control projects.

Claims

1. A method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion, characterized in that, Includes the following steps: A basic dataset for the target area is obtained, which includes basic datasets of data sources obtained by different observation methods. Spatiotemporal alignment is performed on the data sources obtained by different observation methods to generate a basic computational dataset containing multiple spatial computing grids and a discrete hydrological observation sequence set containing multiple discrete hydrological observation stations. Based on the basic computational dataset, the surface water flow connectivity network is extracted, the hydrological connectivity distance between the spatial computing grid and the discrete hydrological observation station is calculated, and a spatial distribution weight matrix is ​​constructed. The spatial distribution weight matrix is ​​used to process the discrete hydrological observation sequence set to perform spatial assimilation, and a boundary mask function is introduced to constrain the calculation region of the interpolation to generate a physically constrained initial field dataset. A one-dimensional and two-dimensional coupled shallow water dynamics calculation framework containing Manning roughness coefficient and lateral weir flow coefficient is constructed by loading the physical constraint initial field dataset and dynamic boundary input conditions. A multi-objective evaluation function is constructed, and the parameters are calibrated by using reverse optimization to iteratively modify the Manning roughness coefficient and the lateral weir flow coefficient. Define the boundaries of multi-scene geomorphic features, call the one-dimensional and two-dimensional coupled shallow water dynamics calculation framework that has completed the parameter calibration to perform forward numerical simulation, extract the spatiotemporal evolution flow field sequence and output quantitative evaluation results.

2. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 1, characterized in that, In the steps of generating the basic computational dataset and the discrete hydrological observation sequence set, performing spatiotemporal alignment on the data sources obtained by the different observation methods includes: The spatial topographic data of the target area, the boundary vector data of water-related projects, and the planar coordinates of the optical satellite remote sensing image sequence are uniformly mapped to the target reference projection coordinate system, and the elevation node values ​​are corrected to a unified vertical elevation reference surface. The computational grid size and time step size of the spatial computational grid are set globally uniformly based on the Corange number constraint conditions of hydrodynamic solutions. The time nodes corresponding to the starting points of typical historical flood events are extracted and set as global time zero. For the discrete hydrological observation sequence set with inconsistent sampling frequencies, a linear interpolation function is used to reconstruct the continuous numerical sequence of the target monitoring location under the time calculation step.

3. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 1, characterized in that, In the step of constructing the spatial distribution weight matrix, extracting the surface water flow connectivity network and calculating the hydrological connectivity distance include: Local depression filling correction is performed on isolated elevation minima of the digital elevation model in the basic computational dataset. The direction of the adjacent grid with the largest spatial elevation gradient is extracted to establish a global flow direction matrix, and a global confluence accumulation matrix is ​​generated by deduction. The global runoff accumulation matrix is ​​clipped using the polygon with the largest water body coverage in historical floods, and the grid node connections with runoff accumulation greater than the connectivity threshold are extracted to generate the surface water flow connectivity network. The discrete hydrological observation stations are mapped as effective mapping anchor points on the surface water flow connectivity network. Adjacent grid nodes are connected as undirected weighted graph structure edges. Graph theory pathfinding is performed along the surface water flow connectivity network. The spatial physical length of the undirected weighted graph structure edges on the shortest topological path is accumulated to generate the hydrological connectivity distance.

4. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 3, characterized in that, In the step of constructing the spatial distribution weight matrix, constructing the spatial distribution weight matrix includes: A distance decay function that takes into account terrain obstacles and energy dissipation is constructed. The initial weights of water level control for each effective mapping anchor point to the spatial computing grid are calculated using the hydrological connectivity distance, the energy decay index that characterizes the surface roughness of the floodplain, and the preset minimum smoothing constant. A local interpolation search radius is set, and the initial weights of the water level control corresponding to the effective mapping anchor points within the local interpolation search radius are accumulated as the total weights. The initial weight of the water level control of a single effective mapping anchor point is divided by the total weights to derive a standardized weight value that satisfies the sum of 1, which is used to generate the spatial distribution weight matrix.

5. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 4, characterized in that, In the step of generating the physical constraint initial field dataset, generating the physical constraint initial field dataset includes: A two-dimensional spatial boundary mask function is established, and the ray crossing algorithm and the mesh overlap area ratio are used to jointly determine whether the spatial calculation grid falls inside the two-dimensional closed space formed by the polygon covering the largest water body of the historical flood, which is used to divide the effective calculation area and the invalid masking area. For each of the spatial computing grids within the effective computing area, the standardized weight values ​​of the effective mapping anchor points associated with the spatial computing grids are multiplied by the corresponding measured water level values ​​and accumulated to obtain the initial assimilated water level values. The elevation data of the digital elevation model is introduced to perform physical bottom constraint verification. When the initial assimilated water level value is detected to be less than the elevation data of the digital elevation model, the water level value of the corresponding spatial computing grid is forcibly corrected to the elevation data of the digital elevation model superimposed with the ultrathin water layer constant, thereby generating a water level surface that has been physically constrained.

6. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 5, characterized in that, In the step of generating the initial physical constraint field dataset, generating the initial physical constraint field dataset further includes: Subtracting the elevation data of the corresponding digital elevation model from the water level surface that has undergone physical constraint processing, the effective water depth value at the current time step is obtained, and an effective water depth distribution matrix is ​​generated. Calculate the effective water depth gradient between adjacent spatial computing grids. When the effective water depth gradient exceeds a preset gradient threshold, construct a spatial computing window containing the surrounding neighborhood for the spatial computing grid and perform mean smoothing filtering. The initial velocity components of the spatial computation grid are forced to zero, and an initial momentum field is established by setting no-slip solid boundary conditions outside the effective computation region. The physically constrained water level, the effective water depth distribution matrix, and the initial momentum field are encapsulated into the physically constrained initial field dataset.

7. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 1, characterized in that, In the step of completing the parameter calibration, the construction of the one-dimensional and two-dimensional coupled shallow water dynamics calculation framework includes: Extract the main channel boundary line and arrange lateral connecting weirs along the main channel boundary line, and topologically bind the spatial computing grid of the one-dimensional channel node and the two-dimensional floodplain area; Based on the explicit finite volume method, the interface water flow transfer flux is estimated by using one- or two-dimensional water level difference, and an interactive flow limiter is set to limit the maximum overflow rate when the preset water level difference threshold is exceeded. A two-dimensional continuity control equation is established to calculate the water mass conservation state of the flow field. The interface water flow transfer flux of the lateral connecting weir is used as the source term input to establish a momentum conservation equation to verify the physical constraints of water flow convection, water level gravity gradient and bed friction on the momentum change of the water body.

8. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 7, characterized in that, In the step of completing the parameter calibration, constructing the multi-objective evaluation function and completing the parameter calibration include: Based on the land use attributes of the underlying surface, the discretized computational domain is divided into spatial functional zones including the main river channel floodplain, densely populated urban residential area, and agricultural ecological floodplain. The root mean square error index of flood peak water level is set for the densely populated urban area and the measured water level and the simulated water level. The relative error index of inundation duration is set for the agricultural ecological floodplain. The spatial intersection ratio error index of simulated inundation range and remote sensing measured water area is set for the physical connectivity of the whole area. The Nash efficiency coefficient index of flow process line is set for the main channel flood discharge area. The range normalization results of each indicator are combined with the preset protection importance weights to form the multi-objective evaluation function; The particle swarm optimization algorithm is used to perform swarm intelligent reverse optimization iteration on the Manning roughness coefficient and the lateral weir flow coefficient. The convergence value of the multi-objective evaluation function is used as the fitness criterion, and the optimal solution vector that meets the convergence verification requirements is extracted to complete the parameter calibration.

9. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 8, characterized in that, In the step of outputting the quantitative evaluation results, performing the forward numerical simulation includes: Based on the preset multi-scene terrain feature boundary positioning of the spatial computing grid subject to engineering intervention, spatial location indicator variables and engineering design elevations are introduced to update and correct the elevation data of the digital elevation model, and local depression filling verification is performed. The corresponding flow process line driving sequence is loaded and input to the one-dimensional channel node, and the physical constraint initial field dataset is extracted as the basic state of the global calculation start time. The one-dimensional and two-dimensional coupled shallow water dynamics calculation framework that has completed the parameter calibration is invoked to perform forward integral calculation, monitor the maximum water flow velocity of the spatial calculation grid in the whole domain and calculate the local Corange number. When the local Corange number exceeds the stability critical value, the time step is automatically halved. If the halved time step is lower than the minimum step threshold, a numerical artificial viscosity term is added to the spatial calculation grid that has exceeded the limit.

10. The method for dynamic reconstruction of flood evolution in floodplains based on multi-source data fusion according to claim 9, characterized in that, In the step of outputting quantitative evaluation results, the output quantitative evaluation results include: Based on the spatiotemporal evolution flow field sequence output by the forward numerical simulation, the spatial computing grids below the effective water depth threshold are filtered out, and the isolated spatial computing grids that are not physically connected to the main channel flood discharge area are eliminated using the hydraulic connectivity topology traversal mechanism to determine the effective flood storage and detention grid set. For each time step node, the product of the instantaneous effective water depth and the spatial plane area of ​​each spatial computing grid in the effective flood storage and detention grid set is summed and integrated to obtain the instantaneous effective flood storage and detention total volume, and the highest value is extracted from the time series as the transient peak effective flood storage and detention total volume. By comparing the output characteristics of numerical simulations under different geomorphic scenarios, the reduction in flood peak water level, the increase in the effective total flood storage volume of the transient peak, and the delay in the spread duration are extracted and range normalization is performed. The comprehensive flood control benefit score is calculated by weighting the entropy weight method.

Images