Method and system for evaluating flood toughness of grid scale basin based on underlying surface function and endurance capability

Abstract

The invention discloses a grid scale basin flood toughness evaluation method and system based on underlying surface functions and endurance capability, and the method comprises the steps: obtaining high-resolution remote sensing image data and POI data, obtaining a fine classification result through multi-level classification processing, and carrying out the calculation through combining with social economic data, thereby obtaining various underlying surface function weight values; performing statistical analysis and gridding calculation on the fine classification result and historical disaster data to obtain a grid unit endurance capacity value; a coupling hydrodynamic model is constructed for simulation calculation, and grid submerging water depth and submerging time data are obtained; calculating a system performance curve, and extracting toughness characteristic values to obtain a grid toughness characteristic value set; and calculating comprehensive toughness based on the toughness characteristic value set, and performing spatial integration in combination with the functional weight values of the underlying surface at different moments to obtain toughness evaluation results of different spatial and temporal scales. Accurate decision support is provided for urban flood control and disaster reduction.

Description

Technical Field

[0001] The present invention belongs to the field of water resources management, and in particular to a grid-scale watershed flood resilience evaluation method and system based on underlying surface function and tolerance. Background Art

[0002] Urban flooding, a global natural disaster, is experiencing a significant increase in frequency and losses due to climate change and accelerated urbanization. In rapidly urbanizing areas, significant changes in the underlying surface area, leading to altered rainfall-runoff relationships, coupled with insufficient drainage systems, are increasing the risk of urban flooding. Therefore, establishing a scientific method for assessing watershed flood resilience and improving urban flood prevention and disaster reduction capabilities are crucial for ensuring safe urban operations and protecting residents' lives and property. Furthermore, resilience assessment results can provide a scientific basis for urban planning and infrastructure development, and have significant practical value in enhancing a city's overall resilience to disasters.

[0003] Currently, scholars at home and abroad have conducted extensive research on urban flood resilience assessment. Traditional assessment methods are primarily based on a single indicator system, quantifying resilience levels by constructing evaluation indicators and weighting systems. Data acquisition relies primarily on statistical data and field surveys, resulting in low spatial resolution. Flood simulation often employs single one- or two-dimensional hydrodynamic models, which struggle to accurately describe flood evolution under complex terrain conditions. Disaster-prone object identification relies primarily on manual interpretation and simple remote sensing classification methods, making it difficult to achieve refined functional identification. Resilience calculations often employ static assessment methods, failing to fully consider the system's dynamic recovery process.

[0004] However, existing research methods still have the following technical problems: First, in terms of underlying surface function identification, traditional methods find it difficult to accurately identify the subdivided functional types of urban land, especially in functional mixed areas, and the classification accuracy is low; second, in terms of tolerance capacity assessment, there is a lack of quantitative methods based on historical disaster data, and threshold determination often relies on expert experience; third, in terms of hydrodynamic simulation, there are stability problems in the exchange of boundary conditions and numerical solutions in the model coupling process, and the computational efficiency is low; fourth, in terms of performance curve analysis, the identification of key time nodes lacks objective mathematical support, especially when the system experiences multiple fluctuations, it is difficult to accurately identify the lowest performance point and the end of recovery; fifth, in terms of multi-scale resilience evaluation, the weight transfer between different scales lacks a scientific theoretical basis, and the reliability of the evaluation results needs to be improved; finally, in terms of temporal dynamics, existing methods find it difficult to reflect the changes in the importance of underlying surface functions in different time periods (such as working days and non-working days, and crop growing season), which affects the practicality of the evaluation results. Summary of the Invention

[0005] The purpose of the invention is to provide a grid-scale watershed flood resilience evaluation method and system based on underlying surface function and tolerance, in order to solve at least one technical problem existing in the prior art.

[0006] The technical solution, a grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance, includes the following steps:

[0007] S1. Obtain high-resolution remote sensing image data and POI data, and obtain detailed classification results through multi-level classification processing; combine with pre-stored socioeconomic data to calculate the functional weight values ​​of various underlying surfaces;

[0008] S2. Statistical analysis and grid calculation are performed on the sub-classification results and pre-stored historical disaster data to obtain the grid unit tolerance value;

[0009] S3. Perform simulation calculations based on a preconfigured coupled hydrodynamic model to obtain grid inundation data, including grid inundation depth and inundation time;

[0010] S4. Based on the grid flooding data and the grid unit tolerance value, calculate the system performance curve, extract the toughness characteristic value, and obtain the toughness characteristic value set of each grid;

[0011] S5. Calculate comprehensive resilience based on the set of resilience characteristic values; perform spatial integration based on the weight values ​​of various underlying surface functions to obtain resilience evaluation results at different spatiotemporal scales.

[0012] A grid-scale watershed flood resilience assessment system based on underlying surface function and tolerance includes:

[0013] at least one processor; and,

[0014] a memory communicatively connected to at least one of the processors; wherein,

[0015] The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the grid-scale watershed flood resilience evaluation method based on underlying surface function and tolerance.

[0016] Beneficial effect: The present invention realizes accurate underlying surface function identification through multi-source data fusion and multi-level classification, then establishes a tolerance assessment system based on historical disaster data, simulates the flood process through coupling hydrodynamic models, and finally extracts resilience characteristics based on system performance curves and realizes multi-scale evaluation; it not only takes into account the differences and spatiotemporal variation characteristics of regional functions, but also includes a comprehensive assessment of the system's resistance and recovery capabilities, improves the scientificity and reliability of the evaluation results, overcomes the shortcomings of traditional evaluation methods in spatial accuracy and temporal dynamics, and provides more accurate technical support for urban flood control and disaster reduction decision-making. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 Flowchart of the method of the present invention.

[0018] Figure 2 This is a flow chart of step S1 of the present invention.

[0019] Figure 3 This is a flow chart of step S2 of the present invention.

[0020] Figure 4 This is a flow chart of step S3 of the present invention.

[0021] Figure 5 This is a flow chart of step S4 of the present invention.

[0022] Figure 6 This is a flow chart of step S5 of the present invention. DETAILED DESCRIPTION

[0023] The following describes the present application in more detail with reference to specific embodiments. Figure 1 As shown, this application proposes a grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance, which includes the following steps:

[0024] S1. Obtain high-resolution remote sensing image data and POI data, and obtain detailed classification results through multi-level classification processing; based on the detailed classification results and pre-stored socioeconomic data, calculate the functional weight values ​​of various underlying surfaces;

[0025] S2. Statistically analyze the sub-classification results with pre-stored historical disaster data to identify the disaster-prone characteristics of various land uses and determine flood tolerance thresholds. Based on the flood tolerance thresholds, grid-based calculations are performed to obtain grid unit tolerance values.

[0026] S3. Read terrain data, pipe network data, and river data to build a coupling model; perform simulation calculations based on the coupling model and pre-stored design rainfall data to obtain grid inundation data, including grid inundation depth and inundation time;

[0027] S4. Based on the grid flooding data and the grid unit tolerance value, calculate the system performance curve, extract the toughness characteristic value, and obtain the toughness characteristic value set of each grid;

[0028] S5. Calculate the comprehensive resilience based on the set of resilience characteristic values; combine the comprehensive resilience with the weight values ​​of various underlying surface functions at different times, perform spatial integration, and obtain resilience evaluation results at different temporal and spatial scales.

[0029] like Figure 2 As shown, according to one aspect of the present application, step S1 is further:

[0030] S11, reading high-resolution remote sensing image data, performing radiation correction processing in sequence to obtain a radiation corrected image; performing geometric correction processing on the radiation corrected image to obtain a geometric corrected image; performing atmospheric correction processing on the geometric corrected image to obtain a corrected remote sensing image;

[0031] S12. Based on the corrected remote sensing image, combined with pre-stored support vector machine training sample data, a support vector machine classifier is trained to obtain a primary classification model; the primary classification model is used to classify the corrected remote sensing image to obtain a primary classification result including paddy fields, dry land, woodland, grassland, urban land, and unused land; urban land data in the primary classification result is read, and combined with pre-stored deep learning training sample data, the urban land in the primary classification result is divided into blocks to obtain urban image block data; data enhancement processing is performed on the urban image block data to obtain enhanced urban image data; the enhanced urban image data is standardized to obtain standardized urban image data; the standardized urban image data is input into a pre-configured deep learning model for training to obtain a refined classification model; the refined classification model is used to perform secondary classification processing on the urban land in the primary classification result to obtain a refined classification result;

[0032] S13. Reading original POI data from a map service, performing a coordinate duplication check on the original POI data to obtain duplicate coordinate data; deleting duplicate points based on the duplicate coordinate data to obtain cleaned POI data; reading a pre-stored POI function classification standard, performing function classification annotation on the cleaned POI data, and generating classified POI data including commercial, residential, industrial, and public facility function types; performing spatial distribution analysis on the classified POI data to obtain POI density data;

[0033] Step S14: Based on the refined classification results, the classified POI data, and the POI density data, data are superimposed according to spatial positions to obtain spatial superposition data; a pre-stored density threshold standard is read, and the POI point density is calculated for each grid in the spatial superposition data to obtain a grid density value; based on the grid density value and the classified POI data, the dominant functional type of each grid is calculated to obtain grid functional data; based on the grid functional data, urban land is divided into commercial land, residential land, industrial land, public facility land, green land, administrative and business land, warehousing and logistics land, transportation land, and cultural and entertainment land to obtain a detailed classification result;

[0034] S15. Obtain socioeconomic data, construct a hierarchical analysis judgment matrix based on the detailed classification results and socioeconomic data, calculate the eigenvectors, and obtain the initial weights; perform a consistency check based on the initial weights to obtain the revised weights; standardize the revised weights based on the pre-stored economic contribution values ​​of various types of land use to obtain the weight values ​​of various underlying surface functions.

[0035] In one embodiment of the present application, the radiation correction method is: radiation brightness L = a·(DN - Lmin) / (Lmax- Lmin) + b·cos(θ); wherein DN is the grayscale value of the original image; Lmin and Lmax are radiation calibration parameters; a and b are atmospheric correction coefficients, a = 1 / (1 + c·τ), b = Ls / (1 + c·τ), τ is the atmospheric optical thickness, Ls is the atmospheric scattered radiance; θ is the solar zenith angle, θ = arccos(sin(φ)·sin(Δ) + cos(φ)·cos(Δ)·cos(ω)), φ is latitude, Δ is the solar declination, and ω is the hour angle.

[0036] POI density calculation method: POI density value D(x, y) = Σ(Ki·Wi·fi(x, y)) / (πr 2 ); where Ki is the kernel function, Ki = (1 - (di / r) 2 ) 2 , di is the distance, r is the search radius; Wi is the POI weight, Wi = ln(1 + Ni / N0), Ni is the number of POIs of this type, N0 is the benchmark number; fi(x, y) is the spatial distribution function, fi(x, y) = exp(-((x-xi) 2 +(y-yi) 2 ) / (2σ 2 )), (xi, yi) are the POI coordinates, σ is the smoothing parameter; the spatial clustering index C = Σ(Di·ln(Di / D*)), Di is the local density, and D* is the average density.

[0037] This embodiment achieves accurate underlying surface function identification and weight calculation through multi-source data fusion and multi-level classification processing, improving the accuracy and spatiotemporal adaptability of underlying surface function identification and providing reliable basic data for subsequent resilience evaluation.

[0038] According to one aspect of the present application, step S11 is further as follows:

[0039] S111, reading high-resolution remote sensing image data, extracting image metadata, and obtaining image parameter data; performing integrity check on the image parameter data to obtain parameter integrity data; performing data completion processing based on the parameter integrity data to obtain complete image parameters;

[0040] S112. Based on the complete image parameters, the solar zenith angle and azimuth angle are calculated to obtain solar angle data; the radiance of the top of the atmosphere is calculated by combining the solar angle data and pre-stored sensor calibration parameters to obtain initial radiance data; noise estimation is performed on the initial radiance data to obtain noise level data; selective denoising is performed based on the noise level data to obtain reduced-noise radiance data; angle correction is performed on the reduced-noise radiance data using the solar angle data to obtain a radiometrically corrected image;

[0041] S113. Based on the radiometric correction image and pre-stored ground control point data, feature matching is performed on the control points to obtain matching point pair data; based on the matching point pair data, geometric transformation parameters are calculated to obtain transformation parameter data; error evaluation is performed on the transformation parameter data to obtain transformation accuracy data; it is determined whether the transformation accuracy data meets a preset threshold requirement; if not, feature matching is performed again; if so, the radiometric correction image is resampled using the transformation parameter data to obtain a geometrically corrected image;

[0042] S114. Read the geometrically corrected image and the meteorological observation data acquired in real time, calculate the atmospheric transmittance using the preconfigured 6S radiation transfer model, and obtain transmittance data; calculate the atmospheric scattered radiance based on the transmittance data, and obtain scattered radiance data; extract the dark pixels in the geometrically corrected image based on the scattered radiance data, and obtain dark image metadata; estimate the aerosol optical thickness based on the dark image metadata, and obtain optical thickness data; construct an atmospheric correction model by combining the transmittance data, scattered radiance data, and optical thickness data, and obtain correction model parameters; use the correction model parameters to correct the geometrically corrected image to obtain a corrected remote sensing image.

[0043] This example achieves high-quality image data correction through a systematic remote sensing image preprocessing process. This not only improves the spectral and spatial accuracy of the images, but also effectively eliminates atmospheric influences through dark pixel extraction and aerosol optical depth estimation, providing a high-quality data foundation for subsequent land use classification. Through rigorous quality control and precise correction processing, the application value of remote sensing images is enhanced and the reliability of subsequent analysis results is guaranteed.

[0044] According to one aspect of the present application, step S13 is further as follows:

[0045] S131. Read the original POI data from the map service, extract the coordinate information to obtain the POI coordinate data; perform projection transformation on the POI coordinate data to obtain unified projection data; read the attribute information in the unified projection data to obtain attribute field data; perform integrity check on the attribute field data to obtain field integrity data; perform data completion processing based on the field integrity data to obtain complete POI data.

[0046] S132. Read the complete POI data, calculate the spatial proximity relationship to obtain proximity data; perform cluster analysis on the proximity data to obtain spatial cluster data; identify duplicate points based on the spatial cluster data to obtain duplicate coordinate data; perform spatial distance calculation on the duplicate coordinate data to obtain distance matrix data; set a merging threshold based on the distance matrix data to obtain merging threshold data; merge the duplicate points according to the merging threshold data to obtain cleaned POI data.

[0047] S133. Read the pre-stored POI function classification standard and cleaned POI data, construct a function classification dictionary to obtain classification dictionary data; perform word segmentation processing on the name field of the cleaned POI data to obtain word segmentation result data; perform semantic matching on the word segmentation result data based on the classification dictionary data to obtain semantic matching data; perform rule filtering on the semantic matching data to obtain function type data; associate the function type data with the cleaned POI data to obtain classified POI data.

[0048] S134. Read the classified POI data, construct a kernel density estimation model to obtain density model parameters; classify and summarize the classified POI data to obtain category statistical data; calculate the density of each type of POI based on the density model parameters to obtain category density data; perform spatial interpolation on the category density data to obtain density interpolation data; clip the density interpolation data according to the scope of the study area to obtain POI density data.

[0049] This example develops a method for cleaning and functional analysis of point-of-interest (POI) data. This method not only addresses quality issues in the original data but also accurately identifies functional types through semantic analysis and spatial statistics, providing reliable auxiliary data for calculating underlying surface functional weights. In particular, in terms of density calculation and spatial interpolation, scientific parameter settings and spatial analysis methods ensure the rationality and continuity of density distribution results.

[0050] like Figure 3 As shown, according to one aspect of the present application, step S2 is further:

[0051] S21. Based on the detailed classification results and pre-stored historical disaster data, disaster impact statistics are collected for each type of land use, and information on the main disaster-affected objects in the disaster records is extracted; based on the information on the main disaster-affected objects, a correspondence table between land use types and disaster-affected objects is constructed, and a list of disaster-affected object types containing the main disaster-affected entities of each type of land use is generated;

[0052] S22. Based on the list of disaster-prone body types and pre-stored historical disaster records, data preprocessing is performed on each type of disaster-prone body to obtain a disaster impact data set; based on the disaster impact data set, a survival analysis model is constructed to obtain a disaster prediction model; the disaster prediction model is used to analyze and process the disaster impact data set, and the survival probability of each type of disaster-prone body under different flooding depths is calculated to obtain a water depth impact probability value; the disaster prediction model is used to analyze and process the disaster impact data set, and the survival probability of each type of disaster-prone body under different flooding durations is calculated to obtain a time impact probability value; based on the water depth impact probability value, the minimum impact water depth and maximum bearing water depth of each type of disaster-prone body are determined to obtain a flood tolerance depth threshold; based on the time impact probability value, the shortest impact time and the longest bearing time of each type of disaster-prone body are determined to obtain a flood tolerance time threshold;

[0053] S23. Based on the detailed classification results, the flooding depth threshold and the flooding time threshold, the study area is divided into grid units of equal size to obtain basic grid data; the detailed classification results are spatially divided according to the basic grid data to obtain grid classification data; based on the grid classification data, a corresponding flooding depth threshold is configured for each grid unit to obtain grid water depth threshold data; a corresponding flooding time threshold is configured for each grid unit to obtain grid time threshold data; based on the grid water depth threshold data and the grid time threshold data, a comprehensive calculation is performed to obtain the grid unit tolerance value.

[0054] In one embodiment of the present application, the survival analysis model is specifically: the hazard function h(t) = h0(t)·exp(Σβi·Xi(t)); where h0(t) is the baseline hazard function, h0(t) = λ·κ·(λt) κ-1 , λ is the scale parameter, κ is the shape parameter; Xi(t) is the time-varying covariate, including the flooding depth D(t) and duration T(t); βi is the regression coefficient, obtained by maximum likelihood estimation, L(β) = Π[h(ti) Δi ·S(ti)], Δi is the censoring indicator variable, Π is the product operator; survival function S(t) = exp(-∫h(u)du); flooding tolerance threshold Dc = inf{D: S(t|D) < Sc}, Sc is the critical survival probability.

[0055] Grid tolerance calculation method: Grid tolerance value C(i, j) = Σ(Wk·Sk·Tk·Ek); where Wk is the land use type weight, Wk = (ak·bk) / (Σak·bk), ak is the area ratio, and bk is the importance coefficient; Sk is the terrain correction factor, Sk = 1 + μ·(▽z / z0) α, ▽z is the elevation gradient, z0 is the reference elevation; Tk is the time influence coefficient, Tk= 1 - exp(-t / τk), t is the flooding duration, τk is the characteristic time; Ek is the environmental adaptability index, Ek = (1 +θ·Dk / D0)·(1 + ω·Pk / P0), Dk is the density of drainage facilities, Pk is the protection level, D0 and P0 are reference values, and θ and ω are adjustment parameters.

[0056] This example combines historical disaster data with land use types to develop a resilience assessment method based on survival analysis. Incorporating survival analysis theory, it considers not only the impact of inundation depth but also the duration of inundation, making the resilience assessment more comprehensive and objective. Through gridding, the assessment results are refined to the grid scale, improving the spatial accuracy and practicality of the assessment and providing more accurate disaster resilience parameters for resilience evaluation.

[0057] According to one aspect of the present application, step S22 is further as follows:

[0058] S221. Read the list of disaster-prone body types and pre-stored historical disaster records, extract the time of disaster occurrence to obtain disaster time series data; extract the disaster extent information to obtain loss extent data; extract the inundation characteristics to obtain inundation feature data; correlate and match the disaster time series data, loss extent data and inundation feature data to obtain a disaster impact data set.

[0059] S222. Read the disaster impact data set, construct the Cox proportional risk model to obtain the basic risk function; calculate the covariate influence coefficient to obtain the influence coefficient data; perform the likelihood ratio test to obtain the model significance data; screen the influence coefficient data based on the model significance data to obtain the optimized coefficient data; substitute the optimized coefficient data into the basic risk function to obtain the disaster prediction model.

[0060] S223. Read the disaster prediction model, perform stratification processing on the water depth data in the disaster impact data set to obtain water depth stratification data; calculate the cumulative risk function of each water depth layer to obtain water depth risk data; calculate the survival function based on the water depth risk data to obtain water depth survival data; perform probability conversion on the water depth survival data to obtain the water depth impact probability value.

[0061] S224. Use the disaster prediction model to stratify the flooding time data in the disaster impact data set to obtain time-stratified data; calculate the cumulative risk function of each time layer to obtain time risk data; calculate the survival function based on the time risk data to obtain time survival data; and perform probability conversion on the time survival data to obtain a time impact probability value.

[0062] S225. Read the water depth impact probability value, construct a probability threshold curve to obtain a water depth threshold curve; perform inflection point analysis on the water depth threshold curve to obtain water depth inflection point data; determine a critical probability value based on the water depth inflection point data to obtain a water depth critical value; extract the minimum impact water depth and the maximum bearing water depth according to the water depth critical value to obtain a flood-resistant water depth threshold.

[0063] S226. Read the time impact probability value, construct a probability threshold curve to obtain a time threshold curve; perform inflection point analysis on the time threshold curve to obtain time inflection point data; determine the critical probability value based on the time inflection point data to obtain a time critical value; extract the shortest impact time and the longest tolerance time according to the time critical value to obtain a flood resistance time threshold.

[0064] This example incorporates survival analysis theory to establish a disaster resilience assessment model based on historical data. This model not only considers the combined effects of flooding depth and duration but also reflects the varying vulnerabilities of different hazard-sustaining bodies through a probabilistic model, enhancing the scientific nature of the resilience assessment. In particular, the objective determination of critical values, through hierarchical analysis and inflection point identification, allows for an objective threshold, avoiding the subjectivity of traditional empirical methods.

[0065] According to one aspect of the present application, step S23 is further:

[0066] S231. Read the detailed classification results, extract the boundary of the study area to obtain boundary range data; calculate the appropriate grid size based on the boundary range data to obtain grid size parameters; use the grid size parameters to regularly divide the study area to obtain initial grid data; perform boundary processing on the initial grid data to obtain boundary grid data; merge the boundary grid data with the internal grid to obtain basic grid data.

[0067] S232. Read the basic grid data and detailed classification results, construct a spatial index to obtain grid index data; calculate the coverage relationship between the grid and the land type based on the grid index data to obtain coverage area data; perform percentage statistics on the coverage area data to obtain type proportion data; determine the dominant type of the grid based on the type proportion data to obtain grid classification data.

[0068] S233, read the grid classification data and the flooding depth threshold, construct the classification mapping relationship to obtain the water depth mapping data; assign a threshold to each grid based on the water depth mapping data to obtain the initial water depth threshold; read the pre-stored terrain undulation data, calculate the elevation change of the grid to obtain the elevation change data; correct the initial water depth threshold based on the elevation change data to obtain the grid water depth threshold data.

[0069] S234, read the grid classification data and the flood resistance time threshold, construct the classification mapping relationship to obtain the time mapping data; assign a threshold to each grid based on the time mapping data to obtain the initial time threshold; read the pre-stored drainage condition data, analyze the drainage capacity of the grid to obtain the drainage capacity data; correct the initial time threshold based on the drainage capacity data to obtain the grid time threshold data.

[0070] S235. Read the grid water depth threshold data and the grid time threshold data, construct a water depth-time tolerance curve to obtain tolerance curve data; normalize the tolerance curve data to obtain normalized tolerance; read the pre-stored disaster tolerance level standard, perform graded assessment on the normalized tolerance to obtain tolerance level data; calculate the comprehensive tolerance based on the tolerance level data to obtain the grid unit tolerance value.

[0071] This example develops a grid-based resilience assessment method. This not only improves the spatial accuracy of the assessment but also makes the results more realistic by considering local factors such as topography and drainage. Specifically, in terms of resilience calculation, a two-dimensional resilience curve is established to express disaster resilience under the dual constraints of water depth and time, providing more comprehensive basic data for resilience assessment.

[0072] like Figure 4 As shown, according to one aspect of the present application, step S3 is further:

[0073] S31. Read the topographic data of the study area, i.e., DEM data, and generate a topographic grid; read the pipe network GIS data and construct a one-dimensional pipe network model; read the river channel data, including river channel cross-section data and river channel plane data, and construct a one-dimensional river channel model; read the surface cover data and construct a two-dimensional surface model; based on the topographic grid, couple the one-dimensional pipe network model, the one-dimensional river channel model, and the two-dimensional surface model through exchange boundary conditions to generate a one- and two-dimensional coupled hydrodynamic model;

[0074] S32. Using pre-stored design rainfall data as input boundary conditions to obtain rainfall boundary data; extracting underlying surface parameters based on the detailed classification results to obtain underlying surface parameter data; constructing a calculation time step based on the one- and two-dimensional coupled hydrodynamic model to obtain calculation step data; constructing an output time interval based on the one- and two-dimensional coupled hydrodynamic model to obtain output interval data; inputting the rainfall boundary data, underlying surface parameter data, calculation step data, and output interval data into the one- and two-dimensional coupled hydrodynamic model to perform model calculation to obtain calculation results; based on the calculation results, recording the waterlogging depth for each calculation grid at each output time to obtain grid inundation depth data; and based on the calculation results, counting the waterlogging duration for each calculation grid to obtain grid inundation time data.

[0075] S33. Based on the grid inundation depth data and the grid inundation time data, perform outlier detection to obtain outlier data; based on the outlier data, eliminate abnormal records to obtain verified inundation depth; perform spatiotemporal continuity analysis on the verified inundation depth to obtain continuity analysis data; based on the continuity analysis data, perform data smoothing to obtain corrected grid inundation data.

[0076] In one embodiment of the present application, the rainfall runoff calculation method is: unit grid flow Q(t) = C·i(t)·A- f(t); where C is the comprehensive runoff coefficient, C = Σ(ci·ai), ci is the underlying surface classification runoff coefficient, and ai is the area ratio; i(t) is the rainfall intensity, i(t) = i0·(1 + k·ln(T / t0)), i0 is the design rainfall intensity, T is the recurrence period, and t0 is the rainfall duration; f(t) is the infiltration rate, f(t) = fc + (f0 - fc)·exp(-kt), fc is the stable infiltration rate, f0 is the initial infiltration rate, and k is the attenuation coefficient; A is the catchment area.

[0077] This example constructs a one- and two-dimensional coupled hydrodynamic model to achieve unified simulation of surface runoff, pipe network drainage, and river channel evolution. This not only accounts for the complexity of urban drainage systems but also enables unified simulation of multi-scale hydrological processes, providing more accurate inundation characteristic data for flood risk assessment.

[0078] According to one aspect of the present application, step S31 is further as follows:

[0079] S311. Read pre-stored pipe network GIS data, perform a topology check, and obtain pipe network connectivity data; perform a breakpoint check on the pipe network connectivity data to obtain pipe network breakpoint data; perform connectivity repair based on the pipe network breakpoint data to obtain corrected pipe network data; perform an attribute integrity check on the corrected pipe network data to obtain pipe network attribute data; perform parameter supplementation and verification on the pipe network attribute data based on pre-stored pipe network parameter standards to obtain complete pipe network data;

[0080] S312. Based on the complete pipe network data, construct the Saint-Venant equations to obtain the initial pipe network equations; discretize the initial pipe network equations to obtain the discrete equations; use the eigendecomposition method to solve the discrete equations to obtain the pipe network eigenvalues; based on the pipe network eigenvalues, construct a discrete format for the finite volume method to obtain numerical format parameters; substitute the numerical format parameters into the discrete equations for solution verification to obtain the pipe network model;

[0081] S313, reading pre-stored river channel cross-sectional data and river channel plane data, performing cross-sectional interpolation calculation to obtain dense cross-sectional data; performing curve fitting on the dense cross-sectional data to obtain a river channel cross-sectional curve; reading pre-stored river channel roughness data, performing segmented assignment to obtain cross-sectional roughness data; constructing a one-dimensional river channel hydrodynamic equation based on the river channel cross-sectional curve and cross-sectional roughness data to obtain a river channel model;

[0082] S314: Read pre-stored DEM data and perform depression-filling processing to obtain depression-filling DEM data; perform slope analysis on the depression-filling DEM data to obtain surface slope data; perform roughness calculation based on the surface slope data and pre-stored surface cover data to obtain surface roughness data; construct a two-dimensional shallow water equation system based on the depression-filling DEM data and the surface roughness data to obtain a surface model;

[0083] S315. Based on the pipe network model, the river model, and the surface model, a grid node correspondence is constructed to obtain node mapping data; based on the node mapping data, a water exchange equation is constructed to obtain exchange flux data; physical constraint processing is performed on the exchange flux data to obtain effective flux data; based on the effective flux data, a boundary condition equation is constructed to obtain boundary condition data;

[0084] S316. Based on the pipe network model, river channel model, surface model and boundary condition data, a model coupling matrix is ​​constructed to obtain coupling matrix data; the coupling matrix data is sparsely processed to obtain sparse matrix data; based on the sparse matrix data, an interactive iterative format is constructed to obtain an iterative solution format; the iterative solution format is applied to the pipe network model, river channel model and surface model to obtain a coupling solver; a case study is performed on the coupling solver to obtain verification result data; based on the verification result data, the solution parameters are optimized to obtain optimized parameter data; the optimized parameter data is configured in the coupling solver to obtain a one- and two-dimensional coupled hydrodynamic model.

[0085] In one embodiment of the present application, the exchange flux Q(t) = Σ(Qij·ηij·φij); where Qij is the inter-node flow rate, Qij = μ·A·(2g·ΔH) 0.5 , μ is the flow coefficient, A is the water flow area, g is the acceleration of gravity, ΔH is the water level difference; ηij is the coupling coefficient, ηij = (1 + exp(-k·(h-hc))) -1 , h is the water depth, hc is the critical water depth, k is the conversion coefficient; φij is the stability factor, φij = min(1, (Δt·sqrt(g / h)) / Δx), Δt is the time step, and Δx is the space step.

[0086] This example establishes a complete one- and two-dimensional coupled hydrodynamic model by integrating pipe network, river channel, and surface models. This model not only overcomes the limitations of a single model but also improves computational efficiency through sparse matrix technology and an interactive iterative format, providing high-precision numerical support for flood simulations. In particular, solver optimization, through case studies and parameter optimization, ensures the computational stability and accuracy of the coupled model.

[0087] According to one aspect of the present application, step S33 is further as follows:

[0088] S331. Perform statistical feature analysis on the grid inundation depth data to obtain water depth mean data and water depth standard deviation data; calculate water depth upper and lower thresholds based on the water depth mean data and water depth standard deviation data to obtain water depth threshold data; perform statistical feature analysis on the grid inundation time data to obtain time mean data and time standard deviation data; calculate time upper and lower thresholds based on the time mean data and time standard deviation data to obtain time threshold data;

[0089] S332: Based on the water depth threshold data, the time threshold data, and the pre-stored terrain elevation data, a physical constraint check is performed on each grid to obtain physical limit violation data; based on the grid submerged water depth data of adjacent grids, the water depth gradient between grids is calculated to obtain water depth gradient data; the water depth gradient data is compared with a preset maximum allowable gradient value to obtain gradient anomaly data;

[0090] S333. Based on the grid flooded water depth data and the grid flooded time data, calculate the water depth change rate at adjacent moments to obtain water depth change rate data; compare the water depth change rate data with a preset maximum allowable change rate to obtain time-varying anomaly data; calculate a spatial correlation index based on the grid flooded water depth data of adjacent grids to obtain spatial correlation data; compare the spatial correlation data with a preset spatial correlation threshold to obtain spatial anomaly data;

[0091] S334. Based on the water depth threshold data, time threshold data, physical overlimit data, gradient anomaly data, time-varying anomaly data and spatial anomaly data, a comprehensive analysis is performed to obtain anomaly level data; the anomaly level data is divided into normal data, slightly anomaly data and severely anomaly data; the slightly anomaly data is subjected to spatiotemporal interpolation processing to obtain interpolation correction data; the severely anomaly data is marked as an invalid value to obtain invalid value marking data; the interpolation correction data and the normal data are merged to obtain the corrected grid inundation data.

[0092] This embodiment improves the accuracy, spatiotemporal continuity, and overall quality of data through multi-dimensional anomaly detection and data correction, enhances the precision and stability of model calculations, and provides a solid data foundation for flood simulation and risk assessment.

[0093] like Figure 5 As shown, according to one aspect of the present application, step S4 is further:

[0094] S41. Based on the corrected grid flooding data and the grid unit tolerance value, calculate the system performance at each moment to generate a system performance time series; extract the time when the performance starts to decline, the time when the performance reaches its lowest point, and the time when the performance recovers from the system performance time series to obtain system performance characteristic time data; wherein the value range of the system performance is 0-1;

[0095] S42. Based on the system performance time series and the system performance characteristic time data, extract the lowest performance value to obtain a robustness value; based on the robustness value, calculate the performance difference before and after the backup drainage facility is activated to obtain a redundancy value; based on the robustness value and the system performance characteristic time data, calculate a resilience value; combine the robustness value, redundancy value, and resilience value to form a set of resilience characteristic values. Wherein the resilience value R3 = (1-Pmin) / (P(te)-Pmin), where Pmin is the lowest performance value and P(te) is the system performance value at the time te when performance recovery ends.

[0096] In one embodiment of the present application, the system performance calculation method is as follows: system performance value P(t) = 1 - Σ[Wi×Fi(t) ×Si×Ki×Ri(t)]; where Wi is the function weight value of the i-th grid, ranging from [0 to 1]; Fi(t) is the flooding impact factor of the i-th grid at time t, Fi(t) = α·(Hi(t) / HTi) β + (1-α)·(Ti(t) / TTi) γ , Hi(t) is the actual flooding depth, HTi is the flooding depth threshold, Ti(t) is the actual flooding time, TTi is the flooding time threshold, α, β, and γ are influence coefficients; Si is the spatial association weight, Si = exp(-λ·∑(dij / D0)), dij is the distance between grid i and adjacent grid j, D0 is the characteristic distance, and λ is the attenuation coefficient; Ki is the key node influence factor, Ki = 1 + μ·∑(Ii·Vi), Ii is the facility importance index, and Vi is the vulnerability index; Ri(t) is the resilience coefficient, Ri(t) = 1 - exp(-θ·t / Tc), θ is the recovery rate parameter, and Tc is the characteristic recovery time.

[0097] Calculation method for toughness characteristic value: comprehensive toughness value R = α·Rob + β·Red + γ·Rec; where Rob is the robustness value, Rob = 1 - |Pmin - P0| / P0, Pmin is the minimum performance value, and P0 is the initial performance value; Red is the redundancy value, Red = ∫(Pr(t) - Pb(t))dt / T, Pr(t) is the performance curve after the backup facility is activated, Pb(t) is the baseline performance curve, and T is the evaluation period; Rec is the recovery value, Rec = ΣVi·exp(-λ·ti), Vi is the recovery rate, Vi = (Pi+1- Pi) / (ti+1 - ti), ti is the time of the recovery stage, and λ is the time attenuation coefficient; α, β, and γ are characteristic weight coefficients, and α+β+γ=1.

[0098] This embodiment can identify the weak links of the system in flood disasters and propose corresponding improvement measures, which helps to enhance the overall flood prevention and disaster reduction capabilities and reduce the losses caused by flood disasters; the generated system performance time series and resilience characteristic value sets provide decision makers with scientific decision support data, which can help decision makers formulate more effective flood prevention and disaster reduction strategies and improve the scientificity and effectiveness of flood response.

[0099] According to one aspect of the present application, step S41 is further as follows:

[0100] S411. Read the corrected grid inundation data and grid unit tolerance value, calculate the inundation water depth ratio of each grid to obtain water depth ratio data; perform spatial autocorrelation analysis on the water depth ratio data to obtain water depth correlation data; construct a spatial weight matrix based on the water depth correlation data to obtain spatial weight data; perform weighted calculation on the spatial weight data and the water depth ratio data to obtain weighted water depth data.

[0101] S412, read the corrected grid inundation data, calculate the inundation time ratio of each grid to obtain time ratio data; read the pre-stored time weight standard, perform graded weighting on the time ratio data to obtain time weight data; perform weighted calculation on the time weight data and the time ratio data to obtain weighted time data.

[0102] S413. Read the weighted water depth data and weighted time data, and perform area normalization processing on each grid to obtain area ratio data; read the weight values ​​of various underlying surface functions, perform multi-dimensional weighted calculation on the area ratio data, weighted water depth data, and weighted time data to obtain the grid performance value.

[0103] S414. Read the grid performance value and reconstruct the time series to obtain performance time series data; perform wavelet transform on the performance time series data to obtain wavelet coefficient data; identify key time nodes based on the wavelet coefficient data to obtain key moment data; perform cluster analysis on the key moment data to obtain moment cluster data; extract the moment when performance starts to decline from the moment cluster data to obtain decline moment data.

[0104] S415. Read the performance time series data, construct a dynamic time warping model to obtain a time warping model; input the performance time series data into the time warping model for pattern matching to obtain pattern matching data; identify the lowest performance point based on the pattern matching data to obtain the lowest point moment data; perform confidence interval analysis on the lowest point moment data to obtain confidence interval data; determine the final performance lowest point moment based on the confidence interval data to obtain the corrected lowest point moment data.

[0105] S416. Read the performance time series data, use the change point detection algorithm to identify the performance recovery inflection point to obtain the inflection point moment data; perform trend analysis on the inflection point moment data to obtain the recovery trend data; determine the recovery end time based on the recovery trend data to obtain the recovery moment data; integrate the decline moment data, the corrected lowest point moment data and the recovery moment data to obtain the system performance characteristic moment data.

[0106] This embodiment proposes a system performance evaluation method based on time series analysis. It not only considers spatial correlation but also uses multiple mathematical methods to ensure the accuracy of time node identification, providing reliable data support for the calculation of resilience eigenvalues. In particular, in performance curve analysis, the combined application of time warping and change point detection effectively avoids the uncertainty in identifying critical moments associated with traditional methods.

[0107] According to one aspect of the present application, step S42 is further as follows:

[0108] S421. Read the system performance time series, smooth the time series to obtain smoothed sequence data; calculate the change rate of the smoothed sequence data to obtain the performance change rate; identify key change points based on the performance change rate to obtain key point data; extract the global minimum value from the key point data to obtain minimum value data; perform confidence assessment on the minimum value data to obtain a robustness value.

[0109] S422. Read the pre-stored drainage facility data, extract the facility activation time to obtain the facility time series data; time-align the facility time series data with the system performance time series to obtain aligned sequence data; calculate the system performance before the facility is activated to obtain the pre-activation performance value; calculate the system performance after the facility is activated to obtain the post-activation performance value; perform difference calculation on the pre-activation performance value and the post-activation performance value to obtain the redundancy value.

[0110] S423. Read the system performance characteristic time data, extract the lowest performance moment to obtain the lowest point time data; extract the performance recovery end moment to obtain the recovery time data; calculate the time difference between the two moments to obtain the recovery cycle data; read the performance recovery amount during this period to obtain the recovery amount data; calculate the ratio of the recovery amount data and the recovery cycle data to obtain the initial recovery rate.

[0111] S424. Read the initial recovery rate, extract the turning points in the recovery process to obtain turning point data; calculate the recovery rate of each stage to obtain segmented rate data; perform weighted averaging on the segmented rate data to obtain weighted rate data; calculate the comprehensive recovery rate based on the weighted rate data to obtain the recovery value.

[0112] S425. Read the robustness value, redundancy value and recovery value, construct a feature vector to obtain feature vector data; standardize the feature vector data to obtain a standardized feature value; read the pre-stored feature importance weight, perform weighted combination on the standardized feature value, and obtain a set of toughness characteristic values.

[0113] This example achieves a quantitative expression of resilience characteristics through a system's eigenvalue extraction method. This not only comprehensively reflects the system's resistance and recovery capabilities, but also enables a comprehensive assessment of resilience characteristics by introducing feature importance weights. In particular, the practical significance of the assessment results is enhanced by considering the nonlinear characteristics of the recovery process in the calculation of resilience.

[0114] like Figure 6 As shown, according to one aspect of the present application, step S5 is further:

[0115] S51. Read pre-stored historical case data, construct a fuzzy relationship matrix, and obtain a fuzzy evaluation matrix; perform membership calculation on the fuzzy evaluation matrix to obtain membership data; calculate subjective weights based on the fuzzy evaluation matrix and the membership data to obtain subjective weight values; calculate entropy values ​​based on the toughness characteristic value set to obtain objective weight values; perform arithmetic averaging on the subjective weight values ​​and the objective weight values ​​to obtain characteristic weight coefficients;

[0116] S52. Perform weighted calculation on each grid based on the characteristic weight coefficient and the toughness characteristic value set to obtain a grid toughness value;

[0117] S53. Based on the functional weight values ​​of various underlying surfaces and the pre-stored time period weight coefficients, the grid resilience values ​​in each sub-basin are weighted averaged to obtain the sub-basin resilience value; the weighted average resilience values ​​of all sub-basins are weighted averaged to obtain the overall resilience value of the basin, that is, the resilience evaluation result.

[0118] In one embodiment of the present application, the characteristic weight coefficient θi = Δ·θsi + (1-Δ)·θoi; where θsi is the subjective weight, θsi = (Πaij 1 / n ) / Σ(Πaij 1 / n ), aij is the element of the expert judgment matrix, and n is the matrix order; θoi is the objective weight, θoi = (1 - Hi / H) / Σ(1 - Hi / H), Hi is the entropy value of the i-th indicator, Hi = -k·Σ(pij·lnpij), pij is the standardized value of the j-th sample of the i-th indicator; Δ is the balance coefficient, Δ = ρ·(1 - σs / σmax), σs is the standard deviation of the subjective weight, σmax is the maximum allowable standard deviation, and ρ is the adjustment parameter.

[0119] This example integrates multi-scale resilience assessment results by combining fuzzy evaluation and entropy weighting. This not only considers differences in regional functional importance but also reflects the impact of temporal variations, making the assessment results more practical and instructive. Spatial integration and temporal weighting enhance the reliability and practicality of the assessment results, providing a scientific basis for urban flood control and disaster reduction decision-making.

[0120] According to one aspect of the present application, the steps for calculating the time period weight coefficient are:

[0121] S53a, obtaining population heat map data at different times (e.g., morning peak, evening peak, holidays, etc.) of the study area, and converting it into raster data consistent with the grid cells;

[0122] S53b. Based on the grid data, the population thermal value of each grid unit in different time periods is calculated as the population activity intensity index of the grid unit. All population activity intensity indices are normalized to determine the time period weight coefficient.

[0123] This embodiment reflects the intensity of human activities in urban areas at different times, and can reflect human demand for different land use types at different times. In comparison, the time scale of the crop growing period is larger, measured in months, while human activities are usually calculated in hours. By considering the crop cycle, the weights of land use types such as agriculture, forestry, animal husbandry, and fishery are increased in some time periods, not just in urban areas.

[0124] According to one aspect of the present application, step S51 is further as follows:

[0125] S511. Read pre-stored historical case data and perform data standardization processing to obtain standardized case data; perform outlier detection on the standardized case data to obtain abnormal case data; eliminate abnormal samples based on the abnormal case data to obtain cleaned case data; perform stratified sampling on the cleaned case data to obtain training sample data.

[0126] S512. Read the training sample data, construct a fuzzy evaluation index system to obtain evaluation index data; perform correlation analysis on the evaluation index data to obtain index correlation data; perform index screening based on the index correlation data to obtain preferred index data; perform fuzzy language variable conversion on the preferred index data to obtain fuzzy language data; construct a fuzzy relationship matrix based on the fuzzy language data to obtain a fuzzy judgment matrix.

[0127] S513, reading the fuzzy evaluation matrix, constructing a fuzzy membership function to obtain membership function data; performing parameter optimization on the membership function data to obtain optimized function parameters; performing membership calculation on the fuzzy evaluation matrix based on the optimized function parameters to obtain membership data.

[0128] S514. Read the toughness characteristic value set, calculate the sample entropy value to obtain characteristic entropy value data; calculate the information contribution based on the characteristic entropy value data to obtain information weight data; normalize the information weight data to obtain an objective weight value.

[0129] S515. Read the fuzzy judgment matrix, construct a judgment matrix to obtain a weight judgment matrix; perform a consistency test on the weight judgment matrix to obtain a consistency index; based on the consistency index, modify the weight judgment matrix to obtain a modified judgment matrix; perform eigenvalue calculation on the modified judgment matrix to obtain a subjective weight value.

[0130] S516: Read the objective weight value and the subjective weight value, perform arithmetic average operation, and obtain the characteristic weight coefficient.

[0131] This example combines fuzzy evaluation and entropy weighting to establish a combined objective and subjective weight determination method. This not only avoids the limitations of a single weighting method but also enhances the scientific nature of the evaluation process through the application of fuzzy theory. In particular, consistency checks and matrix corrections ensure the rationality of subjective weight determination, improving the reliability and practicality of resilience assessment results and providing strong technical support for urban flood control and disaster reduction decision-making.

[0132] In another embodiment of the present application, step S12 can also be: reading the corrected remote sensing image, constructing an SVM classification model, setting kernel function parameters, training the classifier to obtain the SVM classification model, and using the SVM classification model to classify the corrected remote sensing image to obtain a primary classification result including paddy fields, dry land, woodland, grassland, urban land and unused land.

[0133] Step S13 may also be: reading original POI data from a map service, deleting duplicate coordinate points to obtain cleaned POI data, classifying and labeling the cleaned POI data according to functional types such as commercial, residential, industrial, and public facilities, and generating classified POI data.

[0134] Step S15 can also be: reading the detailed classification results and pre-stored socio-economic data, constructing a hierarchical analysis judgment matrix to obtain an initial judgment matrix; performing eigenvector calculation on the initial judgment matrix to obtain an initial weight value; combining pre-stored time period basic data, reading the initial weight value and real-time collected time information data, and performing working day and non-working day time period calculation on the weight value of the urban area to obtain an urban time period weight value; reading the initial weight value and pre-stored growth cycle data, and performing growing period weight calculation on the weight value of the agricultural and forestry ecological area to obtain an ecological time period weight value; merging the urban time period weight value and the ecological time period weight value to obtain a time period weight coefficient; performing a product operation on the initial weight value and the time period weight coefficient to obtain the weight value of each type of underlying surface function.

[0135] By constructing an SVM classification model and classifying the corrected remote sensing images, this embodiment can accurately identify and classify different types of land use, such as paddy fields, dry land, woodlands, grasslands, urban land, and unused land, providing a reliable data foundation for subsequent analysis and decision-making. By deleting duplicate coordinate points and classifying and labeling the cleaned POI data with functional types, the accuracy and completeness of the POI data can be effectively improved, helping to more accurately reflect actual geographic information and provide support for urban planning and management. By constructing a hierarchical analysis judgment matrix and combining it with real-time collected time information data, the weight values ​​of urban areas and agricultural and forestry ecological areas can be dynamically calculated. This can better reflect the actual conditions of different time periods and growth periods and improve the scientificity and rationality of the weight values. By multiplying the initial weight value with the time period weight coefficient, the functional weight values ​​of various underlying surfaces are obtained, which can more comprehensively consider the influence of different factors and improve the accuracy and applicability of the weight values. This embodiment helps to fully understand and evaluate the actual conditions of the study area and provide strong support for scientific decision-making.

[0136] According to one aspect of the present application, step S15 can also be: reading the subclassification results and pre-stored socio-economic data, constructing a fuzzy comparison matrix to obtain an initial fuzzy matrix; performing triangular fuzzy number conversion on the fuzzy numbers in the initial fuzzy matrix to obtain a triangular fuzzy number matrix; reading the triangular fuzzy number matrix, and using the extent analysis method to calculate the fuzzy comprehensive extension value to obtain a fuzzy extension value; performing probability calculation on the fuzzy extension value to obtain a probability matrix; calculating the weight vector based on the probability matrix to obtain a fuzzy weight value; combining the pre-stored time period basic data, reading the fuzzy weight value and the real-time collected time information data, and calculating the weight value of the urban area for working days and non-working days to obtain the urban time period weight value; reading the fuzzy weight value and the pre-stored growth cycle data, and performing growing period weight calculation on the weight value of the agricultural and forestry ecological area to obtain the ecological time period weight value; performing fuzzy arithmetic operation on the urban time period weight value and the ecological time period weight value to obtain the time period weight coefficient; performing fuzzy product operation on the fuzzy weight value and the time period weight coefficient to obtain the weight value of each type of underlying surface function.

[0137] In one embodiment of the present application, the fuzzy weight calculation method is as follows: function weight W = ΣWi·Mi·Ti; where Wi is the basic function weight, Wi = (Σaij / Σ(Σaij)), aij is the fuzzy judgment matrix element; Mi is the time period adjustment coefficient, Mi = β1·exp(-((t-tp) / τ1) 2 ) + β2·exp(-((t-tp) / τ2) 2 ), t is the current time, tp is the peak time, τ1 and τ2 are time scale parameters, β1 and β2 are weight coefficients; Ti is the functional type coefficient, Ti = 1 + γ·(Pi / Pmax) α , Pi is the POI density value, Pmax is the maximum density value, γ is the density influence coefficient, and α is the nonlinear adjustment parameter.

[0138] This embodiment constructs a fuzzy comparison matrix and triangular fuzzy number conversion, adopts the extent analysis method to calculate the fuzzy comprehensive extension value, and performs probability calculation, which can more accurately calculate the weight vector, so as to effectively deal with uncertainty and ambiguity, and improve the accuracy and reliability of weight calculation. Combined with the pre-stored time period basic data and the real-time collected time information data, the weight values ​​of urban areas and agricultural and forestry ecological areas are dynamically calculated, which can better reflect the actual conditions of different time periods and growth periods, and improve the scientificity and rationality of the weight values. Through fuzzy arithmetic operations and fuzzy product operations, the influence of various factors on the weight values ​​can be comprehensively considered, which can more comprehensively reflect the actual conditions of various underlying surface functions and improve the accuracy and applicability of the weight values. Through accurate weight calculation and multi-dimensional analysis, scientific decision support data can be provided to decision makers, helping them to formulate more effective management and planning strategies and improve the ability to cope with complex environments and uncertainties.

[0139] According to one aspect of the present application, step S22 can also be: reading the list of disaster-prone body types and historical disaster records, counting the degree of disaster damage of each type of disaster-prone body at different flooding depths, determining the minimum impact water depth and maximum bearing water depth of each type of disaster-prone body to obtain the flood-resistant water depth threshold, counting the degree of disaster damage of each type of disaster-prone body at different flooding durations, determining the shortest impact time and longest bearing time of each type of disaster-prone body to obtain the flood-resistant time threshold.

[0140] Step S23 can also be: reading the detailed classification results, the flooding depth threshold and the flooding time threshold, dividing the study area into grid units of equal size, dividing the detailed classification results by grid units, assigning corresponding flooding depth and time thresholds to each grid unit, and calculating the grid unit tolerance value.

[0141] In another embodiment of the present application, step S3 may also be:

[0142] S3a, read the pre-stored pipe network GIS data and perform topological structure processing to obtain the pipe network topology data; construct a one-dimensional pipe network hydraulic model based on the pipe network topology data to obtain the pipe network model; read the pre-stored river section data and river plane data, construct a one-dimensional river hydraulic model to obtain the river model; read the pre-stored DEM data, and construct a two-dimensional surface model in combination with the real-time collected surface cover data to obtain the surface model; generate an adaptive calculation grid based on the pipe network model, the river model and the surface model to obtain grid data; read the grid data, and perform terrain feature analysis on the grid to obtain a terrain feature grid; read the terrain feature grid to perform water flow feature analysis to obtain a water flow feature grid; perform grid dynamic optimization based on the water flow feature grid to obtain optimized grid data; exchange the boundary conditions of the pipe network model, the river model and the surface model on the basis of the optimized grid data to obtain a one- and two-dimensional coupled hydrodynamic model.

[0143] S3b. Read the designed rainfall data as the input boundary condition, read the detailed classification results to extract the underlying surface parameters, read the one- and two-dimensional coupled hydrodynamic model, set the calculation time step and output interval, execute the model calculation, record the water depth for each grid cell at each output moment to obtain the grid flooding depth (H(i, t)), and count the water accumulation duration of each grid to obtain the grid flooding time (T(i, t)).

[0144] S3c, read the grid inundation depth (H(i, t)) and grid inundation time (T(i, t)), perform data verification on the calculation results, eliminate abnormal values ​​to obtain the verified inundation depth, perform spatiotemporal continuity analysis on the verified inundation depth, perform data smoothing, and finally obtain the corrected grid inundation data.

[0145] This example uses adaptive grid technology to dynamically optimize grid distribution based on topographic and flow characteristics, improving computational efficiency and accuracy. This example supports flood analysis and assessment at various scales, helping to fully understand and assess the actual conditions in the study area and providing strong support for scientific decision-making.

[0146] In another embodiment of the present application, step S4 may also be:

[0147] S4a. Read the corrected grid inundation data and grid unit tolerance value, calculate the inundation water depth ratio of each grid to obtain water depth ratio data; calculate the inundation time ratio of each grid to obtain time ratio data; calculate the area ratio of each grid to obtain area ratio data; read the weight values ​​of various underlying surface functions, and perform weighted calculation with the water depth ratio data, time ratio data and area ratio data to obtain the grid performance value; summarize the grid performance value P(t) to obtain the system performance time series; extract the performance start decline time from the system performance time series to obtain decline time data; extract the performance lowest point time to obtain the lowest point time data; extract the performance recovery end time to obtain the recovery time data; integrate the decline time data, the lowest point time data and the recovery time data to obtain the system performance characteristic time data. P(t) = 1- Σ(Hi(t) / HTi × Ti(t) / TTi × Ai / AT) × Wi; where P(t) is the system performance value at time t; Hi(t) is the actual flooding depth of the i-th grid at time t; HTi is the flooding depth threshold of the i-th grid; Ti(t) is the actual flooding time of the i-th grid at time t; TTi is the flooding time threshold of the i-th grid; Ai is the area of ​​the i-th grid; AT is the total area of ​​the study area; Wi is the functional weight value of the i-th grid; or P(t) = 1- Σ[(Hi(t) / HTi × Ti(t) / TTi × Ai / AT) × Wi × Si × Ki]; where the newly added parameter Si is the spatial correlation weight, and Ki is the key node influence factor.

[0148] S4b. Read the system performance time series and system performance characteristic moment data, extract the lowest performance value from the system performance time series to obtain the robustness value; calculate the performance difference before and after the standby drainage facility is activated to obtain the redundancy value; read the lowest performance point and the end of performance recovery in the system performance characteristic moment data, calculate the recovery rate to obtain the recovery value; combine the robustness value, redundancy value and recovery value to form a set of resilience characteristic values.

[0149] This example proposes a method for extracting resilience characteristics based on system performance curves. This method not only considers the system's ability to resist disturbances but also its ability to recover, making resilience evaluation more comprehensive and scientific. By accounting for spatial correlation and the influence of key nodes, the spatial representativeness of the evaluation results is enhanced.

[0150] According to another aspect of the present application, step S41 can also be: reading the corrected grid inundation data and grid unit tolerance value, reorganizing the data in time series to obtain time series feature data; constructing an LSTM neural network model to obtain a time series prediction model; inputting the time series feature data into the time series prediction model for training to obtain a trained LSTM model; reading the corrected grid inundation data, calculating the inundation water depth ratio of each grid to obtain water depth ratio data; calculating the inundation time ratio of each grid to obtain time ratio data; calculating the area ratio of each grid to obtain area ratio data; reading the weight values ​​of various underlying surface functions, combining the water depth ratio data, time ratio data and area ratio data to obtain performance characteristic data; inputting the performance characteristic data into the trained LSTM model for prediction to obtain a predicted performance sequence; resampling the predicted performance sequence by time step to obtain a system performance time series; extracting the moment when the performance starts to decline from the system performance time series to obtain decline moment data; extracting the moment when the performance is at its lowest point to obtain the lowest point moment data; extracting the moment when the performance recovers to obtain recovery moment data; integrating the decline moment data, the lowest point moment data and the recovery moment data to obtain system performance characteristic moment data.

[0151] Through accurate system performance prediction and detailed performance characteristic time data, this embodiment can identify the weak links of the system in flood disasters and propose corresponding improvement measures, which helps to improve the overall flood prevention and disaster reduction capabilities and reduce the losses caused by flood disasters.

[0152] In another embodiment of the present application, step S5 may also be:

[0153] S5a. Read historical case data to construct an expert rating matrix, use the hierarchical analysis method to calculate the subjective weight, read the toughness characteristic value set to calculate the entropy value to obtain the objective weight, and perform the arithmetic average of the subjective weight and the objective weight to obtain the characteristic weight coefficient (A1, A2, A3).

[0154] S5b. Read the toughness characteristic value set and characteristic weight coefficients, and calculate the weighted sum ResG = A1×R1 + A2×R2 + A3×R3 for each grid, where ResG is the grid comprehensive toughness value, R1 is the robustness value; R2 is the redundancy value; R3 is the recovery value; A1, A2, and A3 are characteristic weight coefficients; and obtain the grid toughness value.

[0155] S5c. Read the grid resilience value, the weight values ​​of various underlying surface functions, and the preset time period weight coefficient, perform weighted average of the grid resilience values ​​in each sub-basin to obtain the sub-basin resilience value, and perform weighted average of all sub-basin resilience values ​​to obtain the overall basin resilience value.

[0156] According to another aspect of the present application, step S51 can also be: reading pre-stored historical case data, establishing a fuzzy set evaluation factor system to obtain an evaluation factor set; constructing a fuzzy comment set to obtain comment set data; establishing a fuzzy evaluation matrix based on the evaluation factor set and the comment set data to obtain a fuzzy relationship matrix; reading a set of toughness characteristic values, calculating a membership function to obtain membership data; performing entropy calculation on the membership data to obtain an objective weight value; reading a fuzzy relationship matrix, and using the fuzzy Delphi method to calculate a subjective weight value; performing a fuzzy weighted average operation on the subjective weight value and the objective weight value to obtain a characteristic weight coefficient.

[0157] According to one aspect of the present application, a grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance includes the following steps:

[0158] Step S1: collect remote sensing image data of the watershed area, perform preliminary classification of the watershed underlying surface, and obtain paddy fields, dry land, woodland, grassland, urban land, and unused land. Perform secondary classification based on POI data to obtain more detailed land use of urban land, including commercial land, residential land, industrial land, public facility land, green land, administrative and business land, warehousing and logistics land, transportation land, and cultural and entertainment land;

[0159] Step S11: Acquire high-resolution remote sensing images (such as WorldView-3, QuickBird, etc.) covering the area of ​​interest (such as the entire city or a specific area), and perform radiometric correction, geometric correction, and atmospheric correction on the remote sensing images to ensure high accuracy and consistency of the images;

[0160] Step S12: using image processing technology (such as K-means clustering or SVM (support vector machine) classification algorithm) to perform preliminary classification of land objects on the image into major categories such as paddy fields, dry land, woodlands, grasslands, and urban land;

[0161] Step S13: Perform secondary classification of urban areas, using deep learning models (such as convolutional neural networks (CNN)) to extract more detailed spatial features. The trained model identifies buildings, roads, green spaces, etc., and generates preliminary classification results for urban areas, including the general urban land use type.

[0162] Step S14: Perform secondary classification and function analysis based on POI data:

[0163] Step S141: Acquire POI (point of interest) data, including coordinates and category information for commercial, residential, industrial, public facilities, transportation, cultural and entertainment functional areas. POI data can typically be obtained through channels such as OpenStreetMap (OSM), Baidu Maps, and Amap. Clean the POI data to remove duplicate and invalid data, and categorize them by category (e.g., commercial, residential, transportation, etc.).

[0164] Step S142: Perform spatial overlay analysis on the POI data and remote sensing image data. Utilize GIS software (e.g., ArcGIS, QGIS, etc.) to further classify the urban areas identified in the remote sensing image based on the location of the POI data. For example, all POI points labeled "commercial" or "shopping center" may correspond to commercial land in the remote sensing image. For areas in the city that cannot be directly identified through imagery (e.g., mixed commercial and residential areas in the suburbs), the functional identification of the area can be performed using surrounding POI data to ensure a more accurate classification result.

[0165] Step S143: Within the initially identified urban areas (e.g., urban land), further subdivide the land into commercial, residential, and industrial land categories. Combined with the functional information of the POIs, particularly for warehousing and logistics land, cultural and entertainment land, and administrative and commercial land, the specific uses of the corresponding areas are determined through spatial matching and functional analysis.

[0166] Step S15: Verify and optimize the remote sensing image classification results using auxiliary data such as urban planning maps and land use planning maps. These data usually mark detailed land use types in the city, such as commercial areas, residential areas, industrial areas, etc., and are highly authoritative and accurate.

[0167] Step S16: Utilize socioeconomic data, such as population density, traffic flow, and economic activity index, as auxiliary variables to optimize the classification model; for example, areas with high population density are usually residential or commercial land, while areas with low population density may be industrial land or green land;

[0168] Step S17: Classification result verification and accuracy assessment:

[0169] Step S171: Use ground-based measured data, existing land use maps, or relevant statistical data as a reference to evaluate classification accuracy. Verify the accuracy of the classification results by calculating classification accuracy indicators (such as user accuracy, producer accuracy, and Kappa coefficient).

[0170] Step S172: Post-process the classification results to eliminate isolated pixels and misclassified areas to ensure that the final classified image is smooth and coherent; generate GIS layers and visualization maps to display the spatial distribution of various land use types;

[0171] Step S2: Based on different land use types, clarify their functions in social and economic activities, quantify the contribution of each type of function to social and economic activities, and calculate the functional score of each underlying surface type;

[0172] Step S3: Comprehensively consider the main disaster-prone bodies within each underlying surface type, determine the tolerance of each type of disaster-prone body to flood disasters based on historical disaster events, including the flooding depth and flooding duration, and calculate the tolerance capacity of each grid unit in the watershed;

[0173] Step S4: construct a one-dimensional and two-dimensional hydrodynamic coupling model, including a one-dimensional pipe network model, a one-dimensional river channel model, and a two-dimensional surface model, to simulate flood disasters and calculate the inundation depth and inundation time of each grid cell in the watershed area;

[0174] Step S5: Construct a basin flood resilience evaluation method that considers the underlying surface function and tolerance. Based on the inundation depth and inundation time, the system performance curve, and the tolerance of each grid, calculate the robustness, redundancy, and resilience of each grid. Furthermore, subjective and objective weights are used to determine the weight of each characteristic, and the flood resilience of the grid cell is calculated (A1*R1+A2*R2+A3*R3). (Robustness R1 is equal to the minimum performance value, redundancy R2 is equal to the difference between the minimum performance values ​​before and after the backup drainage capacity is used, and resilience R3 is equal to (1-minimum performance value) / (system performance at the time of simulation termination-minimum performance value)).

[0175] Step S6: When further calculating the flood resilience of sub-basin units or the entire basin, it is necessary to systematically consider the functional weight of each underlying surface, and consider the different levels of attention paid to the functions of different underlying surfaces in different time periods. By dynamically adjusting the weight of each underlying surface, the resilience of the basin in different time periods can be obtained, such as working days and non-working days, daytime and nighttime, and different growth stages of crops.

[0176] In one embodiment of the present application, the system performance curve is a classic theory for evaluating system resilience, and has also been widely used in the field of flood disasters. The system performance curve reflects the changes in system performance when the system faces the impact of flood disasters, and the value of system performance varies from 0 (i.e., total loss of system performance) to 1 (i.e., no loss). The initial service level is 1, and it starts to decline from ts after the extreme rainfall event, reaches the minimum value of pf at tps, recovers from tpe, and fully recovers to the full service level at te. The system performance curve is used to quantify the flood resilience of the system. The change process of flood resilience under an extreme rainfall event is as follows: the initial system performance is 1, and it starts to decline from ts after the extreme rainfall event, reaches the minimum value of p(t) at tps, and gradually recovers until it is fully recovered at tps.

[0177] Flood severity, which represents the total amount of system performance loss during the entire flood event, is calculated as follows: Sev = 1 / t n ∫ t=0 tn [1- p(t)]dt; where Sev is the severity of flood, t n is the simulation duration, p(t) is the system performance value at time t; the flood severity Sev is approximated as a rectangular area, and the calculation formula of the system resilience is: Sev=(V TF / V TI )×(t f / t n ); Res0=1-Sev=1-(V TF / V TI )×(t f / t n ); where V TF is the total amount of accumulated water, V TI is the total amount of water flowing into the drainage system, t f is the average waterlogging time, t n is the total simulation time, and Res0 is the system resilience.

[0178] The existing solution provides a fixed value for the elasticity of the system, which fails to reflect the spatiotemporal changes of elasticity. This embodiment evaluates the spatiotemporal process of elasticity in a more detailed way, specifically: Res w =∑ t=1 T ∑ i=1 N Res w (i, t) = ∑ t=1 T ∑ i=1 N [1- sev w (i, t)]; where Res wis the flood resilience of the basin; w (i, t) is the flood resilience of the i-th grid point at time t; Sev w (i, t) is the flood severity at the i-th grid point at time t; T is the total simulation time; N is the number of grids.

[0179] When calculating grid system resilience indices, many existing studies use a waterlogging depth threshold as the decision threshold. That is, when the inundation depth of grid i is less than the inundation threshold, the system performance of grid i is 1; when the inundation depth of grid i exceeds the inundation threshold, the system performance of grid i becomes 0. However, in many cases, grid system performance begins to decline even at shallow inundation depths. For example, a waterlogging depth of 2 cm on urban roads can cause vehicles to slip, while shallow waterlogging (e.g., 5 cm) on pedestrian walkways can impact pedestrians. Furthermore, existing studies do not consider the impact of waterlogging duration. However, different land use types have different requirements for waterlogging duration due to their different functions. For example, in urban areas, the lower threshold for the maximum allowable water receding time is 1 hour, while the upper threshold is 24 hours. For cultivated land, the flood tolerance threshold ranges from 2 to 6 days, depending on the crop type. Therefore, traditional resilience calculation methods that consider a single threshold result in a sudden change in system resilience at the water depth threshold. Furthermore, the failure to consider the impact of flood tolerance duration leads to an overestimation of resilience, which fails to reflect the true performance of the system. In view of this, this embodiment proposes two thresholds based on flooding depth and flooding time to quantify resilience, namely upper threshold and lower threshold. The improved system resilience calculation formula is as follows:

[0180] Sev(i, t) = Sev H (i, t) × Sev T (i, t);

[0181]

[0182] Among them Sev H (i, t) is the flood severity caused by the inundation depth of the i-th grid point at time t; Sev T (i, t) is the flood severity caused by the flooding time at the i-th grid point at time t; H(i, t) is the flooding depth of the i-th grid point at time t; H min (i) is the upper threshold of the flooding depth of the i-th grid point; H max (i) is the lower threshold of the flooding depth of the i-th grid; T(i, T) is the flooding time of the i-th grid point at time T; T min (i) is the lower threshold of the time when the i-th grid enters the water; T max (i) is the upper limit threshold of flooding time of grid i.

[0183] At the watershed scale, land use types were divided into five categories: urban land, paddy fields, dry land, grassland, and forest land. The flooding tolerance depth and flooding tolerance time of each land use type were determined. The results are as follows:

[0184] When the land use type is urban land, the flooding depth is 15~50cm and the flooding time is 1~24h; when the land use type is paddy field, the flooding depth is 30~60cm and the flooding time is 48~144h; when the land use type is dry land, the flooding depth is 5~10cm and the flooding time is 48~96h; when the land use type is grassland, the flooding depth is 2~6cm and the flooding time is 144~432h; when the land use type is forest land, the flooding depth is 10~30cm and the flooding time is 240~960h.

[0185] In addition, in order to compare the recovery capabilities of different entities under the influence of flood disasters, this embodiment proposes the recovery capacity (R c ) calculation index, which is defined as the ratio of the resilience recovery value to the resilience loss value. Resilience recovery refers to the resilience at the end of the simulation (Res E ) and minimum toughness (Res min ), the toughness loss value refers to the toughness at the beginning of the simulation (Res s ) and the difference between the minimum toughness. The formula is as follows: R c =( Res s - Res min ) / ( Res E - Res min ).

[0186] This embodiment takes into account the functions of different underlying surfaces, calculates functional weights, and calculates basin resilience, which is more reasonable than conventional mean statistics. It takes into account the tolerance of different underlying surfaces and introduces inundation time, making the calculation results more reasonable and improving the accuracy. Through the improved resilience calculation method, the resilience at different stages is used as the characteristic value of resilience, and then the resilience is obtained by weighting, taking into account different characteristics.

[0187] According to one aspect of the present application, a grid-scale watershed flood resilience assessment system based on underlying surface function and tolerance includes:

[0188] at least one processor; and,

[0189] a memory communicatively connected to at least one of the processors; wherein,

[0190] The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the grid-scale watershed flood resilience evaluation method based on underlying surface function and tolerance as described in any of the above embodiments.

[0191] This application addresses the issue of accurate recognition of underlying surface functions by employing two-stage classification and POI-assisted recognition. A support vector machine is used for primary classification, dividing land into broad categories such as paddy fields and dry land. A deep learning model is then used for secondary, refined classification of urban land. In particular, data augmentation and standardization improve the generalization capabilities of the deep learning model. By introducing POI data for assisted recognition, and through spatial overlay analysis and density calculation, combined with preset density threshold standards, accurate identification of urban functions such as commercial and residential areas is achieved, effectively resolving the classification challenge of mixed-function areas.

[0192] Survival analysis theory was introduced to quantitatively assess resilience. By constructing a Cox proportional hazards model and calculating the probability of survival under different inundation conditions based on historical disaster data, the authors then objectively determined the flooding depth and time thresholds for different hazard-bearing bodies using probability threshold curves and inflection point analysis. This approach overcomes the subjectivity of traditional empirical judgments and establishes a data-driven, scientific assessment system.

[0193] To address the computational stability issues associated with model coupling, a coupling strategy was employed. Specifically, by constructing node mappings, establishing water exchange equations, and applying physical constraints, the rationality of the exchange flux was ensured. Furthermore, computational efficiency was improved by constructing a sparse matrix and employing an interactive iterative format. In particular, through case studies and parameter optimization, the computational stability of the coupled model was further ensured.

[0194] To identify characteristic points on performance curves, a time series feature recognition method based on multiple mathematical methods was developed. This method reconstructs the time series using wavelet transforms and identifies key nodes using wavelet coefficients. A dynamic time warping model is used for pattern matching, combined with confidence interval analysis to determine the lowest performance point. A change point detection algorithm is used to identify recovery inflection points, improving the accuracy and reliability of identifying key time nodes.

[0195] To address the multi-scale weight transfer problem, a combined objective and subjective weight determination method was established. Objective weights were calculated using the entropy method, subjective weights were determined using fuzzy hierarchical analysis, and the final characteristic weight coefficients were obtained using fuzzy weighted averaging. During the scale conversion process, the underlying surface function weight and time period weight coefficients were considered, achieving scientific weight transfer and ensuring consistency of evaluation results across different scales.

[0196] To address the problem of expressing temporal dynamics, a time period weight coefficient was introduced. By combining basic time period data with temporal information data, urban time period weights for weekdays and non-weekdays, as well as ecological time period weights for different growing seasons, were calculated. Dynamic adjustment of the time period weight coefficients allows for the expression of the time-varying characteristics of the underlying surface functional importance, improving the practicality of the evaluation results.

[0197] This invention integrates multiple technical fields, including remote sensing image processing, POI data analysis, hydrodynamic simulation, and resilience assessment, to construct a complete watershed flood resilience assessment system. First, accurate underlying surface function identification is achieved through multi-source data fusion and multi-level classification. Then, a tolerance assessment system is established based on historical disaster data. Flood processes are simulated using a one- and two-dimensional coupled hydrodynamic model. Finally, resilience characteristics are extracted based on system performance curves and multi-scale evaluation is achieved. This "function identification-disaster bearing capacity-inundation simulation-resilience assessment" technical approach not only takes into account regional functional differences and spatiotemporal variations, but also includes a comprehensive assessment of the system's resistance and recovery capabilities, improving the scientific nature and reliability of the assessment results. In particular, in terms of refined simulation and evaluation at the grid scale, through methods such as adaptive grid optimization, multi-model coupling, and multi-level weight allocation, the shortcomings of traditional evaluation methods in terms of spatial accuracy and temporal dynamics are overcome, providing more accurate technical support for urban flood prevention and disaster reduction decision-making.

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

Claims

1. A grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance, characterized by: The steps include: S1. Obtain high-resolution remote sensing image data and POI data, and obtain detailed classification results through multi-level classification processing; combine with pre-stored socioeconomic data to calculate the weight values ​​of various underlying surface functions; S2. Statistical analysis and grid calculation are performed on the sub-classification results and pre-stored historical disaster data to obtain the grid unit tolerance value; S3, performing simulation calculation based on a preconfigured coupled hydrodynamic model to obtain grid flooding data, including grid flooding depth and flooding time; S4. Based on the grid flooding data and the grid unit tolerance value, the system performance curve is calculated, the toughness characteristic value is extracted, and the toughness characteristic value set of each grid is obtained; S5. Calculate the comprehensive resilience based on the set of resilience characteristic values; perform spatial integration based on the weight values ​​of various underlying surface functions to obtain resilience evaluation results at different time and space scales.

2. The grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance according to claim 1 is characterized in that: Step S1 is further as follows: S11, reading high-resolution remote sensing image data, and sequentially performing radiation correction processing, geometric correction processing, and atmospheric correction processing to obtain a corrected remote sensing image; S12. Based on the corrected remote sensing image, a primary classification model is constructed and classified to obtain a primary classification result; the urban land data therein is read, and block, data enhancement, standardization and deep learning model training are performed to obtain a refined classification model; based on this, a secondary classification is performed on the urban land in the primary classification result to obtain a refined classification result; S13, read POI data, perform coordinate repeatability check, delete duplicate points, functional classification annotation and spatial distribution analysis to obtain POI density data; S14. Based on the refined classification results and POI density data, data superposition, POI point density calculation, dominant function type calculation and urban land use division are performed to obtain refined classification results; S15. Based on the detailed classification results, a hierarchical analysis judgment matrix is ​​constructed and consistency check and standardization are performed to obtain the functional weight values ​​of each type of underlying surface.

3. The grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance according to claim 2 is characterized in that: Step S2 is further as follows: S21. Based on the sub-classification results and pre-stored historical disaster data, statistics on disaster impacts are performed and a corresponding relationship table between land use types and disaster-affected objects is constructed to generate a list of disaster-affected object types; S22. Based on the list of disaster-prone body types and pre-stored historical disaster records, data preprocessing is performed and a survival analysis model is constructed to calculate the probability values ​​of water depth and time impact, and obtain a flood-resistant water depth threshold and a flood-resistant time threshold; S23. Divide the study area into grid units of equal size, perform spatial division on the subdivision results and configure corresponding flooding depth threshold and flooding time threshold, perform comprehensive calculations, and obtain the grid unit tolerance value.

4. The grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance according to claim 3 is characterized in that: Step S3 is further as follows: S31, read the pipe network GIS data, build a one-dimensional pipe network model; read the river data, build a one-dimensional river model; read the surface cover data, build a two-dimensional surface model; based on the terrain data of the study area, couple the one-dimensional pipe network model, the one-dimensional river model and the two-dimensional surface model to generate a one-dimensional and two-dimensional coupled hydrodynamic model; S32, based on the one-two-dimensional coupled hydrodynamic model, calculating the step length data and the interval data; combining the pre-stored design rainfall data and the detailed classification results to perform model calculations to obtain grid flooding depth data and grid flooding time data; S33. Based on the grid flooding depth data and the grid flooding time data, outlier detection, removal of abnormal records, spatiotemporal continuity analysis and data smoothing are performed to obtain corrected grid flooding data.

5. The grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance according to claim 4 is characterized in that: Step S4 is further as follows: S41, based on the grid flooding data and the grid unit tolerance value, calculate the system performance, extract the time when the performance starts to decline, the time when the performance reaches the lowest point, and the time when the performance recovers, and obtain the system performance characteristic time data; The value range of system performance is 0-1; S42. Based on the system performance characteristic time data, the lowest performance value is extracted, and the performance difference and recovery value before and after the standby drainage facilities are activated are calculated, and the results are combined to form a set of resilience characteristic values.

6. The grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance according to claim 5 is characterized in that: Step S5 is further as follows: S51, read the pre-stored historical case data, construct a fuzzy relationship matrix and calculate the membership and subjective weight to obtain the subjective weight value; calculate the entropy value based on the toughness characteristic value set to obtain the objective weight value; Perform arithmetic average calculation on the subjective weight value and the objective weight value to obtain the characteristic weight coefficient; S52, performing weighted calculation on each grid based on the characteristic weight coefficient and the toughness characteristic value set to obtain a grid toughness value; S53. Based on the weight values ​​of various underlying surface functions and the pre-stored time period weight coefficients, the grid resilience values ​​in each sub-basin are weighted averaged to obtain the sub-basin resilience value; the weighted average resilience values ​​of all sub-basins are weighted averaged to obtain the overall resilience value of the basin, that is, the resilience evaluation result.

7. The grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance according to claim 6 is characterized in that: Step S11 is further as follows: S111, read high-resolution remote sensing image data, extract image metadata, perform integrity check and data completion processing, and obtain complete image parameters; S112, calculating the solar zenith angle and azimuth angle based on the complete image parameters; performing noise estimation, selective denoising and angle correction in combination with the pre-stored sensor calibration parameters to obtain a radiation corrected image; S113, based on the radiation correction image and the pre-stored ground control point data, feature matching, geometric transformation parameter calculation and error evaluation are performed to obtain transformation accuracy data; it is determined whether it meets the preset threshold requirement, if not, feature matching is performed again; if so, the radiation correction image is resampled to obtain a geometric correction image; S114, reading the geometrically corrected image and the meteorological observation data acquired in real time, calculating the atmospheric transmittance and the atmospheric scattered radiance, extracting dark pixels, constructing an atmospheric correction model, and obtaining a corrected remote sensing image.

8. The grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance according to claim 6 is characterized in that: Step S31 is further as follows: S311, read the pipe network GIS data, perform topological relationship check, breakpoint detection, connectivity repair, attribute integrity check, parameter supplement and verification, and obtain complete pipe network data; S312, based on the complete pipe network data, construct the Saint-Venant equations, perform discretization processing and solve them, and obtain the pipe network characteristic values; Based on the characteristic values ​​of the pipe network, the discrete format of the finite volume method is constructed, and the solution is verified to obtain the pipe network model; S313, reading pre-stored river channel cross-sectional data and river channel plane data, performing cross-sectional interpolation calculation and curve fitting, and obtaining a river channel cross-sectional curve; reading pre-stored river channel roughness data, performing segmented assignment, and obtaining cross-sectional roughness data; constructing a river channel model based on the river channel cross-sectional curve and cross-sectional roughness data; S314, reading pre-stored DEM data, performing depression filling processing, slope analysis and roughness calculation, and constructing a surface model; S315. Based on the pipe network model, river channel model and surface model, a water exchange equation is constructed, physical constraint processing is performed, and boundary condition data is obtained; S316. Based on the pipe network model, river model, surface model and boundary condition data, a model coupling matrix is ​​constructed, and sparse, interactive iterative format, example verification and optimization solution parameter processing are performed to obtain a one- and two-dimensional coupled hydrodynamic model.

9. The grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance according to claim 6 is characterized in that: Step S41 is further as follows: S411, reading grid inundation data and grid unit tolerance value, calculating inundation depth ratio, constructing a spatial weight matrix, and obtaining weighted depth data; S412, reading grid flooding data, calculating flooding time ratio, performing hierarchical weighting and weighted calculation in combination with pre-stored time weight standards, and obtaining weighted time data; S413, reading weighted water depth data and weighted time data, performing area normalization processing and multi-dimensional weighted calculation, and obtaining a grid performance value; S414, reading grid performance values, performing time series reconstruction, wavelet transform, key time node identification and cluster analysis, and obtaining performance time series data and drop time data; S415, reading performance time series data, building a dynamic time warping model, identifying the lowest performance point and performing confidence interval analysis to obtain the lowest point moment data; S416, reading performance time series data, identifying the performance recovery inflection point, and determining the recovery time data; The descent time data, the lowest point time data and the recovery time data are integrated to obtain the system performance characteristic time data.

10. A grid-scale watershed flood resilience assessment system based on underlying surface function and tolerance, characterized by: include: at least one processor; as well as, a memory communicatively connected to at least one of the processors; wherein, The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the grid-scale watershed flood resilience assessment method based on underlying surface function and tolerance as described in any one of claims 1 to 9.

Images