Grid-scale Watershed Flood Resilience Evaluation Method and System Based on Underlying Surface Function and Tolerance

Abstract

The present invention discloses a grid-scale watershed flood resilience evaluation method and system based on the function and tolerance of underlying surfaces. The method includes obtaining high-resolution remote sensing image data and POI data, obtaining a fine classification result through multi-level classification processing, and calculating the function weight values of various underlying surfaces in combination with socioeconomic data; performing statistical analysis and grid calculation on the fine classification result and historical disaster data to obtain the tolerance capacity values of grid cells; constructing a coupled hydrodynamic model for simulation calculation to obtain grid inundation depth and inundation time data; calculating the system performance curve, extracting the resilience characteristic values, and obtaining a set of grid resilience characteristic values; calculating the comprehensive resilience based on the set of resilience characteristic values, and performing spatial integration in combination with the function weight values of underlying surfaces at different times to obtain the resilience evaluation results at different spatio-temporal scales. The present invention provides accurate decision-making support for urban flood control and disaster reduction.

Description

Technical Field

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

[0002] As one of the global natural disasters, urban flood disasters show a significant upward trend in their occurrence frequency and losses caused with the acceleration of climate change and urbanization. Especially in rapidly urbanizing areas, the significant change in the underlying surface leads to a change in the rainfall-runoff relationship, and coupled with the insufficient capacity of the drainage system, the urban flood risk is continuously increasing. Therefore, establishing a scientific basin flood resilience evaluation method and improving the urban flood prevention and mitigation capabilities are of great significance for ensuring the safe operation of cities and the safety of residents' lives and property. At the same time, the resilience evaluation results can provide a scientific basis for urban planning and infrastructure construction, and have important practical value for enhancing the overall disaster resistance ability of cities.

[0003] At present, scholars at home and abroad have carried out a large number of studies on urban flood resilience evaluation. Traditional evaluation methods mainly rely on a single index system, and quantify the resilience level by constructing an evaluation index and weight system. In terms of data acquisition, they mainly rely on statistical data and field surveys, with a relatively low spatial resolution. In flood simulation, a single one-dimensional or two-dimensional hydrodynamic model is mostly used, which is difficult to accurately describe the flood evolution process under complex terrain conditions. In the identification of disaster-bearing bodies, it mainly relies on manual interpretation and simple remote sensing classification methods, making it difficult to achieve refined function identification. In resilience calculation, static evaluation methods are often used, and the dynamic recovery process of the system is not fully considered.

[0004] However, the existing research methods still have the following technical problems: First, in terms of underlying surface function identification, traditional methods are difficult to accurately identify the subdivision function types of urban land, especially in areas with mixed functions, and the classification accuracy is relatively low; second, in terms of tolerance capacity evaluation, there is a lack of a quantitative method based on historical disaster data, and the determination of thresholds often relies on expert experience; third, in terms of hydrodynamic simulation, there are stability problems in the boundary condition exchange and numerical solution during the model coupling process, and the calculation efficiency is relatively low; fourth, in terms of performance curve analysis, the identification of key time nodes lacks objective mathematical method support, especially when the system shows multiple fluctuations, it is difficult to accurately identify the lowest point of performance and the end time 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 time dynamics, existing methods are difficult to reflect the changes in the importance of underlying surface functions in different time periods (such as weekdays and non-weekdays, crop growth periods), affecting the practicality of the evaluation results. Summary of the Invention

[0005] Objective of the invention: To provide a grid-scale watershed flood resilience evaluation method and system based on the function and tolerance of underlying surfaces, with the expectation of solving at least one technical problem existing in the prior art.

[0006] Technical solution: A grid-scale watershed flood resilience evaluation method based on the function and tolerance of underlying surfaces, comprising the following steps:

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

[0008] S2. Conduct statistical analysis and grid calculation on the fine classification result and the pre-stored historical disaster data to obtain the tolerance capacity value of grid cells.

[0009] S3. Conduct simulation calculation based on a pre-configured coupled hydrodynamic model to obtain grid inundation data, including grid inundation water depth and inundation time.

[0010] S4. Based on the grid inundation data and the tolerance capacity value of grid cells, calculate the system performance curve, extract the resilience characteristic values, and obtain the set of resilience characteristic values for each grid.

[0011] S5. Calculate the comprehensive resilience based on the set of resilience characteristic values; conduct spatial integration in combination with the function weight values of various underlying surfaces to obtain the resilience evaluation results at different spatio-temporal scales.

[0012] A grid-scale watershed flood resilience evaluation system based on the function and tolerance of underlying surfaces, comprising:

[0013] At least one processor; and,

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

[0015] The memory stores instructions executable 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 the function and tolerance of underlying surfaces.

[0016] Beneficial effects: The present invention realizes accurate identification of the function of underlying surfaces through multi-source data fusion and multi-level classification, then establishes a tolerance capacity evaluation system based on historical disaster data, simulates the flood process through a coupled hydrodynamic model, and finally extracts the resilience characteristics and realizes multi-scale evaluation based on the system performance curve; not only considers the differences in regional functions and spatio-temporal variation characteristics, but also includes a comprehensive evaluation of the system's resistance and recovery capabilities, improves the scientificity and reliability of the evaluation results, overcomes the deficiencies of traditional evaluation methods in terms of spatial accuracy and time dynamics, and provides more accurate technical support for urban flood control and disaster reduction decision-making. Description of the Drawings

[0017] Figure 1 This is a flowchart of the method of the present invention.

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

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

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

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

[0022] Figure 6 This is a flowchart of step S5 of the present invention. Detailed Description of the Invention

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

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

[0025] S2. Conduct statistical analysis on the fine classification result and pre-stored historical disaster data to identify the characteristics of disaster-bearing bodies of various land uses and determine the flood tolerance threshold; calculate the tolerance ability value of grid cells based on the flood tolerance threshold through grid calculation;

[0026] S3. Read terrain data, pipe network data, and river data to construct a coupling model; conduct 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. Calculate the system performance curve based on the grid inundation data and the tolerance ability value of grid cells, extract the resilience characteristic value, and obtain a set of resilience characteristic values for each grid;

[0028] S5. Calculate the comprehensive resilience based on the set of resilience characteristic values; combine the comprehensive resilience and the function weight values of various underlying surfaces at different times for spatial integration to obtain the resilience evaluation results at different spatio-temporal scales.

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

[0030] S11. Read the high-resolution remote sensing image data, perform radiometric correction processing in sequence to obtain a radiometrically corrected image; perform geometric correction processing on the radiometrically corrected image to obtain a geometrically corrected image; perform atmospheric correction processing on the geometrically corrected image to obtain a corrected remote sensing image;

[0031] S12. Based on the corrected remote sensing image, combine with the pre-stored support vector machine training sample data to train a support vector machine classifier to obtain a primary classification model; use the primary classification model to classify the corrected remote sensing image to obtain a primary classification result including paddy fields, dry land, forest land, grassland, urban land, and unused land; read the urban land data in the primary classification result, combine with the pre-stored deep learning training sample data, perform block processing on the urban land in the primary classification result to obtain urban image block data; perform data augmentation processing on the urban image block data to obtain enhanced urban image data; perform normalization processing on the enhanced urban image data to obtain normalized urban image data; input the normalized urban image data into a pre-configured deep learning model for training to obtain a fine classification model; use the fine classification model to perform secondary classification processing on the urban land in the primary classification result to obtain a fine classification result;

[0032] S13. Read the original POI data from the map service, perform coordinate repeatability check on the original POI data to obtain duplicate coordinate data; based on the duplicate coordinate data, delete duplicate points to obtain cleaned POI data; read the pre-stored POI function classification standard, perform function classification annotation on the cleaned POI data to generate classified POI data including commercial, residential, industrial, and public facility function types; perform spatial distribution analysis on the classified POI data to obtain POI density data;

[0033] Step S14. Based on the fine classification result, classified POI data, and POI density data, perform data overlay according to spatial location to obtain spatial overlay data; read the pre-stored density threshold standard, calculate the POI point density for each grid in the spatial overlay data to obtain grid density values; based on the grid density values and classified POI data, calculate the dominant function type of each grid to obtain grid function data; according to the grid function data, divide urban land 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 social and economic data, based on the detailed classification result and social and economic data, construct an analytic hierarchy process judgment matrix, calculate the eigenvector to obtain the initial weight; based on the initial weight, perform consistency test to obtain the corrected weight; based on the pre-stored economic contribution values of various types of land, perform normalization processing on the corrected weight to obtain the functional weight values of various underlying surfaces.

[0035] In one embodiment of the present application, the radiation correction method: Radiance L = a·(DN - Lmin) / (Lmax - Lmin) + b·cos(θ); where DN is the gray value of the original image; Lmin and Lmax are radiometric 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 the 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 this type of POI, N0 is the reference number; fi(x, y) is the spatial distribution function, fi(x, y) = exp(-((x - xi) 2 +(y - yi) 2 ) / (2σ 2 ))), (xi, yi) is the POI coordinate, σ 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 realizes accurate identification and weight calculation of the underlying surface function through multi-source data fusion and multi-level classification processing. It improves the accuracy and spatio-temporal adaptability of the underlying surface function identification and provides reliable basic data for subsequent resilience evaluation.

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

[0039] S111. Read high-resolution remote sensing image data, extract image metadata, and obtain image parameter data; perform integrity check on the image parameter data to obtain parameter integrity data; based on the parameter integrity data, perform data complementation processing to obtain complete image parameters;

[0040] S112. Calculate the solar zenith angle and azimuth angle based on the complete image parameters to obtain solar angle data; combine the solar angle data and the pre-stored sensor calibration parameters to calculate the top-of-atmosphere radiance and obtain the initial radiance data; perform noise estimation on the initial radiance data to obtain the noise level data; perform selective denoising processing based on the noise level data to obtain the denoised radiance data; perform angle correction on the denoised radiance data using the solar angle data to obtain the radiometrically corrected image;

[0041] S113. Based on the radiometrically corrected image and the pre-stored ground control point data, perform feature matching on the control points to obtain the matching point pair data; calculate the geometric transformation parameters based on the matching point pair data to obtain the transformation parameter data; perform error evaluation on the transformation parameter data to obtain the transformation accuracy data; determine whether the transformation accuracy data meets the preset threshold requirements. If not, perform feature matching again; if so, perform resampling processing on the radiometrically corrected image using the transformation parameter data to obtain the geometrically corrected image;

[0042] S114. Read the geometrically corrected image and the meteorological observation data obtained in real time, calculate the atmospheric transmittance using the pre-configured 6S radiative transfer model to obtain the transmittance data; calculate the atmospheric scattered radiance based on the transmittance data to obtain the scattered radiance data; extract the dark pixels in the geometrically corrected image based on the scattered radiance data to obtain the dark pixel data; estimate the aerosol optical depth based on the dark pixel data to obtain the optical depth data; combine the transmittance data, the scattered radiance data, and the optical depth data to construct an atmospheric correction model and obtain the correction model parameters; perform correction processing on the geometrically corrected image using the correction model parameters to obtain the corrected remote sensing image.

[0043] In this embodiment, through the systematic remote sensing image preprocessing process, high-quality image data correction is achieved. It not only improves the spectral and spatial accuracy of the image, but also effectively eliminates the atmospheric influence through dark pixel extraction and aerosol optical depth estimation, providing a high-quality data basis for subsequent land use classification. Through strict quality control and precise correction processing, the application value of the remote sensing image is improved, and the reliability of subsequent analysis results is ensured.

[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 the unified projection data; read the attribute information in the unified projection data to obtain the attribute field data; perform integrity check on the attribute field data to obtain the field integrity data; perform data completion processing based on the field integrity data to obtain the complete POI data.

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

[0047] S133. Read the pre-stored POI function classification criteria and the cleaned POI data, construct a function classification dictionary to obtain classification dictionary data; perform word segmentation 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; perform classification summary on the classified POI data to obtain category statistics 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; crop the density interpolation data according to the research area range to obtain POI density data.

[0049] In this embodiment, a method for cleaning and functional analysis of Point of Interest (POI) data is developed. It not only solves the quality problems in the original data, but also realizes the accurate identification of function types through semantic analysis and spatial statistics, providing reliable auxiliary data for the calculation of underlying surface function weights. Especially in density calculation and spatial interpolation, through scientific parameter settings and spatial analysis methods, the rationality and continuity of the density distribution results are ensured.

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

[0051] S21. Based on the fine classification results and the pre-stored historical disaster data, conduct disaster impact statistics on each type of land use, and extract the main disaster-affected object information in the disaster-affected records; based on the main disaster-affected object information, construct a correspondence table between land use types and disaster-affected objects, and generate a list of disaster-bearing body types including the main disaster-bearing entities of each type of land use.

[0052] S22. Based on the inventory of disaster-bearing body types and pre-stored historical disaster records, perform data preprocessing on each type of disaster-bearing body to obtain a disaster impact dataset; based on the disaster impact dataset, construct a survival analysis model to obtain a disaster prediction model; use the disaster prediction model to analyze and process the disaster impact dataset, calculate the survival probability of each type of disaster-bearing body at different inundation depths, and obtain the water depth impact probability value; use the disaster prediction model to analyze and process the disaster impact dataset, calculate the survival probability of each type of disaster-bearing body at different inundation durations, and obtain the time impact probability value; based on the water depth impact probability value, determine the minimum impact water depth and the maximum bearing water depth of each type of disaster-bearing body to obtain the flood tolerance water depth threshold; based on the time impact probability value, determine the shortest impact time and the longest bearing time of each type of disaster-bearing body to obtain the flood tolerance time threshold.

[0053] S23. Based on the fine classification results, the flood tolerance water depth threshold, and the flood tolerance time threshold, divide the study area into grid cells of equal size to obtain basic grid data; divide the fine classification results spatially according to the basic grid data to obtain grid classification data; based on the grid classification data, configure the corresponding flood tolerance water depth threshold for each grid cell to obtain grid water depth threshold data; configure the corresponding flood tolerance time threshold for each grid cell to obtain grid time threshold data; based on the grid water depth threshold data and the grid time threshold data, perform comprehensive calculations to obtain the grid cell tolerance capacity value.

[0054] In an embodiment of the present application, the survival analysis model is specifically: the risk function h(t) = h0(t)·exp(Σβi·Xi(t)); where h0(t) is the baseline risk function, h0(t) =λ·κ·(λt) κ-1 , λ is the scale parameter, κ is the shape parameter; Xi(t) is the time-varying covariate, including the inundation depth D(t) and the 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; the survival function S(t) = exp(-∫h(u)du); the flood tolerance threshold Dc = inf{D: S(t|D) < Sc}, Sc is the critical survival probability.

[0055] Grid tolerance capacity calculation method: The grid tolerance capacity 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, bk is the importance coefficient; Sk is the terrain correction coefficient, Sk = 1 + μ·(▽z / z0) α,▽z is the elevation gradient, z0 is the reference elevation; Tk is the time impact coefficient, Tk = 1 - exp(-t / τk), t is the inundation duration, and τk is the characteristic time; Ek is the environmental adaptability index, Ek = (1 + θ·Dk / D0)·(1 + ω·Pk / P0), Dk is the drainage facility density, Pk is the protection level, D0 and P0 are reference values, and θ and ω are adjustment parameters.

[0056] In this embodiment, historical disaster data is combined with land use types to establish a tolerance assessment method based on survival analysis. Combining the survival analysis theory, not only the inundation depth is considered, but also the influence of the inundation duration is taken into account, making the tolerance assessment more comprehensive and objective. Through grid processing, the assessment results are refined to the grid scale, improving the spatial accuracy and practicability of the assessment, and providing more accurate disaster-bearing capacity 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-bearing body types and the pre-stored historical disaster records, extract the disaster occurrence time to obtain disaster time series data; extract the disaster-affected degree information to obtain loss degree data; extract the inundation characteristics to obtain inundation characteristic data; perform correlation matching on the disaster time series data, loss degree data, and inundation characteristic data to obtain a disaster impact dataset.

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

[0060] S223. Read the disaster prediction model, perform stratified processing on the water depth data in the disaster impact dataset to obtain water depth stratified data; calculate the cumulative risk function for 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 water depth impact probability values.

[0061] S224. Use the disaster prediction model, perform stratified processing on the inundation time data in the disaster impact dataset to obtain time stratified data; calculate the cumulative risk function for each time layer to obtain time risk data; calculate the survival function based on the time risk data to obtain time survival data; perform probability conversion on the time survival data to obtain time impact probability values.

[0062] S225. Read the water depth influence probability value, construct a probability threshold curve to obtain a water depth threshold curve; conduct an inflection point analysis on the water depth threshold curve to obtain water depth inflection point data; determine the critical probability value based on the water depth inflection point data to obtain the water depth critical value; extract the minimum influence water depth and the maximum bearable water depth according to the water depth critical value to obtain the flood tolerance water depth threshold.

[0063] S226. Read the time influence probability value, construct a probability threshold curve to obtain a time threshold curve; conduct an 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 the time critical value; extract the shortest influence time and the longest bearable time according to the time critical value to obtain the flood tolerance time threshold.

[0064] In this embodiment, by introducing the survival analysis theory, a disaster-bearing capacity evaluation model based on historical data is established. It not only considers the comprehensive influence of the inundation depth and duration, but also reflects the vulnerability differences of different disaster-bearing bodies through a probability model, improving the scientificity of the tolerance capacity evaluation. Especially in the determination of the critical value, through hierarchical analysis and inflection point identification, the objective determination of the threshold is realized, avoiding the subjectivity of the traditional empirical method.

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

[0066] S231. Read the fine classification result, extract the boundary of the research area to obtain boundary range data; calculate an appropriate grid size based on the boundary range data to obtain grid size parameters; perform regular meshing on the research area using the grid size parameters 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 the fine classification result, 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; conduct a proportion statistics on the coverage area data to obtain type proportion data; determine the dominant type of the grid according to the type proportion data to obtain grid classification data.

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

[0069] S234. Read the grid classification data and the flood tolerance time threshold, construct the classification mapping relationship to obtain the time mapping data; allocate thresholds 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 the water depth-time tolerance curve to obtain the tolerance curve data; perform normalization processing on the tolerance curve data to obtain the normalized tolerance; read the pre-stored disaster-bearing capacity level standard, perform hierarchical evaluation on the normalized tolerance to obtain the tolerance level data; calculate the comprehensive tolerance ability based on the tolerance level data to obtain the grid unit tolerance ability value.

[0071] In this embodiment, a grid-based tolerance ability evaluation method is developed. It not only improves the spatial accuracy of the evaluation, but also makes the evaluation results more in line with the actual situation by considering local factors such as terrain and drainage. Especially in the calculation of tolerance, by establishing a two-dimensional tolerance curve, the expression of the disaster-bearing capacity under the dual constraints of water depth and time is realized, providing more comprehensive basic data for resilience evaluation.

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

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

[0074] S32. Use the pre-stored design rainfall data as the input boundary condition to obtain the rainfall boundary data; extract the underlying surface parameters based on the fine classification result to obtain the underlying surface parameter data; construct the calculation time step based on the one-two dimensional coupled hydrodynamic model to obtain the calculation step data; construct the output time interval based on the one-two dimensional coupled hydrodynamic model to obtain the output interval data; input the rainfall boundary data, the underlying surface parameter data, the calculation step data and the output interval data into the one-two dimensional coupled hydrodynamic model for model calculation to obtain the calculation result; based on the calculation result, record the waterlogging depth at each calculation grid at each output moment to obtain the grid inundation water depth data; based on the calculation result, count the waterlogging duration for each calculation grid to obtain the grid inundation time data;

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

[0076] In an embodiment of the present application, the rainfall-runoff calculation method is: the unit grid flow rate Q(t) = C·i(t)·A - f(t); where C is the comprehensive runoff coefficient, C = Σ(ci·ai), ci is the runoff coefficient of the underlying surface classification, 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, 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, k is the attenuation coefficient; A is the catchment area.

[0077] In this embodiment, by constructing a two-dimensional coupled hydrodynamic model, the unified simulation of surface runoff, pipe network drainage and river evolution is realized. It not only considers the complexity of the urban drainage system, but also realizes the unified simulation of multi-scale hydrological processes, providing more accurate inundation characteristic data for flood disaster risk assessment.

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

[0079] S311. Read the pre-stored pipe network GIS data, perform topological relationship check to obtain pipe network connectivity data; perform breakpoint detection on the pipe network connectivity data to obtain pipe network breakpoint data; based on the pipe network breakpoint data, perform connectivity repair to obtain corrected pipe network data; perform attribute integrity check on the corrected pipe network data to obtain pipe network attribute data; based on the pre-stored pipe network parameter standard, perform parameter supplementation and verification on the pipe network attribute data 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; perform discretization processing on the initial pipe network equations to obtain discrete equations; use the eigenvalue decomposition method to solve the discrete equations to obtain pipe network eigenvalues; based on the pipe network eigenvalues, construct the discrete format of 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. Read the pre-stored river cross-section data and river plane data, perform cross-section interpolation calculation to obtain dense cross-section data; perform curve fitting on the dense cross-section data to obtain the river cross-section curve; read the pre-stored river roughness data, perform segmented assignment to obtain cross-section roughness data; based on the river cross-section curve and cross-section roughness data, construct a one-dimensional river hydrodynamic equation to obtain a river model;

[0082] S314. Read the pre-stored DEM data, perform depression filling processing to obtain filled DEM data; perform slope analysis on the filled DEM data to obtain surface slope data; based on the surface slope data and the pre-stored land cover data, perform roughness calculation to obtain surface roughness data; based on the filled DEM data and surface roughness data, construct a two-dimensional shallow water equations set to obtain a surface model;

[0083] S315. Based on the pipe network model, river model and surface model, construct the corresponding relationship of grid nodes to obtain node mapping data; based on the node mapping data, construct a water volume exchange equation to obtain exchange flux data; perform physical constraint processing on the exchange flux data to obtain effective flux data; based on the effective flux data, construct a boundary condition equation to obtain boundary condition data;

[0084] S316. Based on the pipe network model, river model, surface model and boundary condition data, construct a model coupling matrix to obtain coupling matrix data; perform sparsification processing on the coupling matrix data to obtain sparse matrix data; based on the sparse matrix data, construct an interactive iteration format to obtain an iterative solution format; apply the iterative solution format to the pipe network model, river model and surface model to obtain a coupled solver; perform case verification on the coupled solver to obtain verification result data; based on the verification result data, optimize the solution parameters to obtain optimized parameter data; configure the optimized parameter data into the coupled solver to obtain a one-two dimensional coupled hydrodynamic model.

[0085] In an embodiment of the present application, the exchange flux Q(t) = Σ(Qij·ηij·φij); where Qij is the flow rate between nodes, Qij = μ·A·(2g·ΔH) 0.5 , μ is the flow coefficient, A is the cross-sectional area of flow, g is the acceleration due to gravity, and Δ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, and 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 spatial step.

[0086] In this embodiment, a complete one - two - dimensional coupled hydrodynamic model is established by integrating pipe network, river channel and surface models. It not only overcomes the limitations of a single model, but also improves the calculation efficiency through sparse matrix technology and interactive iteration format, providing high - precision numerical method support for flood process simulation. Especially in terms of solver optimization, through case verification and parameter optimization, the calculation stability and accuracy of the coupled model are ensured.

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

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

[0089] S332. Based on the water depth threshold data, time threshold data and pre - stored terrain elevation data, conduct physical constraint inspection on each grid to obtain physical over - limit data; calculate the water depth gradient between grids based on the grid inundation water depth data of adjacent grids to obtain water depth gradient data; compare the water depth gradient data with the preset maximum allowable gradient value to obtain gradient anomaly data;

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

[0091] S334. Based on the water depth threshold data, time threshold data, physical over - limit data, gradient anomaly data, time - varying anomaly data and spatial anomaly data, conduct comprehensive analysis to obtain anomaly level data; classify the anomaly level data into normal data, minor anomaly data and serious anomaly data; conduct spatio - temporal interpolation processing on the minor anomaly data to obtain interpolation correction data; mark the serious anomaly data as invalid values to obtain invalid value marking data; merge the interpolation correction data and normal data to obtain the corrected grid inundation data.

[0092] Through multi - dimensional anomaly detection and data correction in this embodiment, the accuracy, spatio - temporal continuity and overall quality of the data are improved, the accuracy and stability of model calculation are enhanced, providing a solid data foundation for flood simulation and risk assessment.

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

[0094] S41. Based on the corrected grid inundation data and the grid cell tolerance values, calculate the system performance for each moment to generate a system performance time series; extract the moment when the performance starts to decline, the moment of the lowest performance point, and the moment when the performance recovery ends from the system performance time series to obtain the system performance characteristic moment 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 moment data, extract the lowest performance value to obtain the robustness value; based on the robustness value, calculate the performance difference before and after the activation of the standby drainage facilities to obtain the redundancy value; based on the robustness value and the system performance characteristic moment data, calculate the recovery value; combine the robustness value, the redundancy value, and the recovery value to form a set of resilience characteristic values. Wherein the recovery value R3 = (1 - Pmin) / (P(te) - Pmin), Pmin is the lowest performance value, and P(te) is the value of the system performance at the moment te when the performance recovery ends.

[0096] In an embodiment of the present application, the system performance calculation method: the system performance value P(t) = 1 - Σ[Wi×Fi(t) ×Si×Ki×Ri(t)]; wherein Wi is the functional weight value of the i-th grid, and the value range is [0, 1]; Fi(t) is the inundation influence factor of the i-th grid at moment t, Fi(t) = α·(Hi(t) / HTi) β + (1 - α)·(Ti(t) / TTi) γ , Hi(t) is the actual inundation depth, HTi is the flood tolerance depth threshold, Ti(t) is the actual inundation time, TTi is the flood tolerance time threshold, and α, β, γ are influence coefficients; Si is the spatial correlation 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, 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 of resilience characteristic value: Comprehensive resilience 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 standby facility is enabled, Pb(t) is the reference performance curve, and T is the evaluation period; Rec is the recoverability value, Rec = ΣVi·exp(-λ·ti), Vi is the recovery rate, Vi = (Pi+1 - Pi) / (ti+1 - ti), ti is the time at the recovery stage, and λ is the time decay 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 improve the overall flood control and disaster reduction capabilities and reduce the losses caused by flood disasters; the generated system performance time series and resilience characteristic value set provide scientific decision-making support data for decision-makers, and can help decision-makers formulate more effective flood control and disaster reduction strategies, and improve the scientificity and effectiveness of dealing with flood disasters.

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

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

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

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

[0103] S414. Read the grid performance value, perform time series reconstruction 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 critical moment data; perform clustering analysis on the critical moment data to obtain moment clustering data; extract the moment when the performance starts to decline from the moment clustering 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 point of performance based on the pattern matching data to obtain lowest point moment data; perform confidence interval analysis on the lowest point moment data to obtain confidence interval data; determine the final lowest point moment of performance based on the confidence interval data to obtain corrected lowest point moment data.

[0105] S416. Read the performance time series data, use a change point detection algorithm to identify the performance recovery inflection point to obtain inflection point moment data; perform trend analysis on the inflection point moment data to obtain recovery trend data; determine the end time of recovery based on the recovery trend data to obtain recovery moment data; integrate the decline moment data, the corrected lowest point moment data, and the recovery moment data to obtain 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 ensures the accuracy of time node identification through various mathematical methods, providing reliable data support for the calculation of resilience characteristic values. Especially in the analysis of performance curves, through the combined application of time warping and change point detection, it effectively avoids the uncertainty in the identification of critical moments in 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, perform smoothing processing on the time series to obtain smoothed sequence data; calculate the change rate of the smoothed sequence data to obtain 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 evaluation 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 facility time series data; perform time alignment on the facility time series data and the system performance time series to obtain aligned sequence data; calculate the system performance before the facility is activated to obtain pre-activation performance value; calculate the system performance after the facility is activated to obtain post-activation performance value; perform difference calculation on the pre-activation performance value and the post-activation performance value to obtain a redundancy value.

[0110] S423. Read the data of the system performance characteristics at a certain moment, extract the moment of the lowest performance to obtain the lowest point time data; extract the moment when the performance recovery ends to obtain the recovery time data; calculate the time difference between the two moments to obtain the recovery period data; read the amount of performance recovery during this period to obtain the recovery amount data; calculate the ratio of the recovery amount data and the recovery period data to obtain the initial recovery rate.

[0111] S424. Read the initial recovery rate, extract the turning points during the recovery process to obtain the turning point data; calculate the recovery rates of each stage to obtain the segmented rate data; perform weighted averaging on the segmented rate data to obtain the 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 the feature vector data; perform normalization processing on the feature vector data to obtain the normalized feature value; read the pre-stored feature importance weights, and perform weighted combination on the normalized feature values to obtain the set of toughness characteristic values.

[0113] In this embodiment, through the eigenvalue extraction method of the system, the quantitative expression of the toughness characteristics is realized. It not only comprehensively reflects the resistance and recovery capabilities of the system, but also realizes the comprehensive evaluation of the toughness characteristics through the introduction of feature importance weights. Especially in the aspect of recovery calculation, by considering the non-linear characteristics of the recovery process, the practical significance of the evaluation results is improved.

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

[0115] S51. Read the pre-stored historical case data, construct a fuzzy relation matrix to obtain a fuzzy evaluation matrix; perform membership degree calculation on the fuzzy evaluation matrix to obtain membership degree data; calculate the subjective weight based on the fuzzy evaluation matrix and the membership degree data to obtain the subjective weight value; calculate the entropy value based on the set of toughness characteristic values to obtain the objective weight value; perform arithmetic mean operation on the subjective weight value and the objective weight value to obtain the characteristic weight coefficient;

[0116] S52. Based on the characteristic weight coefficient and the set of toughness characteristic values, perform weighted calculation on each grid to obtain the grid toughness value;

[0117] S53. Based on the functional weight values of various underlying surfaces and the pre-stored time period weight coefficients, perform weighted average calculation on the grid toughness values within each sub-watershed to obtain the sub-watershed toughness value; perform weighted average calculation on all sub-watershed toughness values to obtain the overall watershed toughness value, that is, the toughness evaluation result.

[0118] In an 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 index, Hi = -k·Σ(pij·lnpij), pij is the standardized value of the j-th sample of the i-th index; Δ 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 embodiment realizes the integration of multi-scale resilience evaluation results by comprehensively applying fuzzy evaluation and entropy weight method. It not only considers the importance differences of regional functions but also reflects the influence of time changes, making the evaluation results more practically guiding. Through spatial integration and time weighting, the reliability and practicality of the evaluation results are improved, providing a scientific basis for urban flood control and disaster reduction decision-making.

[0120] According to one aspect of the present application, the calculation steps of the time period weight coefficient are as follows:

[0121] S53a. Obtain the population heat map data of the research area at different times (such as morning rush hour, evening rush hour, holidays, etc.), and convert it into raster data consistent with the grid cells;

[0122] S53b. Based on the raster data, calculate the population heat value of each grid cell at different time periods as the population activity intensity index of the grid cell, and normalize all the population activity intensity indices to determine the time period weight coefficient.

[0123] This embodiment reflects the human activity intensity at different times in the urban area, which can reflect the needs of humans for different land use types at different times. In contrast, the time scale of the crop growth period is larger, calculated 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 the urban area.

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

[0125] S511. Read the pre-stored historical case data, 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 a correlation analysis on the evaluation index data to obtain index correlation data; perform index screening based on the index correlation data to obtain optimal index data; perform a fuzzy linguistic variable transformation on the optimal index data to obtain fuzzy linguistic data; construct a fuzzy relation matrix based on the fuzzy linguistic data to obtain a fuzzy judgment matrix.

[0127] S513. Read the fuzzy judgment matrix, construct a fuzzy membership function to obtain membership function data; perform parameter optimization on the membership function data to obtain optimized function parameters; perform a membership degree calculation on the fuzzy judgment matrix based on the optimized function parameters to obtain membership degree data.

[0128] S514. Read the set of toughness characteristic values, calculate the sample entropy value to obtain characteristic entropy value data; calculate the information contribution degree based on the characteristic entropy value data to obtain information weight data; perform a normalization process on 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; perform a correction on the weight judgment matrix based on the consistency index to obtain a corrected judgment matrix; perform an eigenvalue calculation on the corrected judgment matrix to obtain a subjective weight value.

[0130] S516. Read the objective weight value and the subjective weight value, perform an arithmetic mean operation to obtain a characteristic weight coefficient.

[0131] This embodiment combines fuzzy evaluation and the entropy weight method to establish an objective-subjective combined weight determination method. It not only avoids the limitations of a single weight method but also improves the scientific nature of the evaluation process through the application of fuzzy theory. Especially in the calculation of subjective weights, through consistency testing and matrix correction, it ensures the rationality of weight determination, improves the reliability and practicality of the toughness evaluation results, and provides 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: Read the corrected remote sensing image, construct an SVM classification model, set the kernel function parameters, train the classifier to obtain the SVM classification model, and use the SVM classification model to perform classification processing on the corrected remote sensing image to obtain a primary classification result including paddy fields, dry land, forest land, grassland, urban land, and unused land.

[0133] Step S13 can also be: Read the original POI data from the map service, delete duplicate coordinate points to obtain cleaned POI data, and classify and label the cleaned POI data according to functional types such as commerce, residence, industry, and public facilities to generate classified POI data.

[0134] Step S15 can also be: reading the fine classification results and the pre-stored socioeconomic data, constructing an analytic hierarchy process judgment matrix to obtain an initial judgment matrix; performing eigenvector calculation on the initial judgment matrix to obtain an initial weight value; combining the pre-stored period basic data, reading the initial weight value and the time information data collected in real time, and calculating the weight values of the urban area for working days and non-working days to obtain the urban period weight value; reading the initial weight value and the pre-stored growth cycle data, and calculating the weight value of the growth period for the agroforestry ecological area to obtain the ecological period weight value; combining the urban period weight value and the ecological period weight value to obtain a period weight coefficient; performing a multiplication operation on the initial weight value and the period weight coefficient to obtain the functional weight values of various underlying surfaces.

[0135] In this embodiment, by constructing an SVM classification model and classifying the corrected remote sensing images, different types of land use, such as paddy fields, dry land, forests, grasslands, urban land, and unused land, can be accurately identified and classified, providing a reliable data basis for subsequent analysis and decision-making. By deleting duplicate coordinate points and classifying and labeling the function types of the cleaned POI data, the accuracy and integrity of the POI data can be effectively improved, helping to more accurately reflect the actual geographical information and providing support for urban planning and management. By constructing an analytic hierarchy process judgment matrix and combining the time information data collected in real time, the weight values of the urban area and the agroforestry ecological area can be dynamically calculated, which can better reflect the actual situation in different time periods and growth periods, and improve the scientificity and rationality of the weight values. By performing a multiplication operation on the initial weight value and the period weight coefficient to obtain the functional weight values of various underlying surfaces, the influence of different factors can be more comprehensively considered, and the accuracy and applicability of the weight values can be improved. This embodiment helps to comprehensively understand and evaluate the actual situation of the research area and provides strong support for scientific decision-making.

[0136] According to one aspect of the present application, step S15 can also be: reading the refined classification result and the pre-stored socioeconomic 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 calculating the fuzzy comprehensive extent value by using the extent analysis method to obtain the fuzzy extent value; performing possibility degree calculation on the fuzzy extent value to obtain a possibility degree matrix; calculating a weight vector based on the possibility degree matrix to obtain a fuzzy weight value; combining the pre-stored period basic data, reading the fuzzy weight value and the real-time collected time information data, and calculating the weight values of the urban area for working days and non-working days to obtain the urban period weight value; reading the fuzzy weight value and the pre-stored growth cycle data, and calculating the weight value of the growth period of the agroforestry ecological area to obtain the ecological period weight value; performing fuzzy arithmetic operation on the urban period weight value and the ecological period weight value to obtain a period weight coefficient; performing fuzzy product operation on the fuzzy weight value and the period weight coefficient to obtain the functional weight values of various underlying surfaces.

[0137] In one embodiment of the present application, the fuzzy weight calculation method: functional weight W = ΣWi·Mi·Ti; where Wi is the basic functional weight, Wi = (Σaij / Σ(Σaij)), aij is the element of the fuzzy judgment matrix; Mi is the 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, τ2 are time scale parameters, β1, β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 non-linear adjustment parameter.

[0138] In this embodiment, by constructing a fuzzy comparison matrix and converting triangular fuzzy numbers, using the extent analysis method to calculate the fuzzy comprehensive expansion value and perform possibility degree calculation, the weight vector can be calculated more accurately. This can effectively handle uncertainty and ambiguity, and improve the accuracy and reliability of weight calculation. Combining the pre-stored period-based data and the real-time collected time information data, dynamically calculating the weight values of urban areas and agroforestry ecological areas can better reflect the actual situation in 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 situation 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-making support data can be provided for decision-makers, helping decision-makers 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-bearing body types and historical disaster records, counting the disaster-affected degrees of various disaster-bearing bodies at different inundation depths, determining the minimum impact water depth and the maximum bearing water depth of each type of disaster-bearing body to obtain the inundation depth threshold, counting the disaster-affected degrees of various disaster-bearing bodies at different inundation durations, and determining the shortest impact time and the longest bearing time of each type of disaster-bearing body to obtain the inundation time threshold.

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

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

[0142] S3a. Reading the pre-stored pipe network GIS data for topological structure processing to obtain pipe network topological data; constructing a one-dimensional pipe network hydraulic model based on the pipe network topological data to obtain a pipe network model; reading the pre-stored river cross-section data and river plane data, constructing a one-dimensional river hydraulic model to obtain a river model; reading the pre-stored DEM data, combining with the real-time collected surface cover data to construct a two-dimensional surface model to obtain a surface model; generating an adaptive calculation grid based on the pipe network model, river model, and surface model to obtain grid data; reading the grid data, performing terrain feature analysis on the grid to obtain terrain feature grids; reading the terrain feature grids to perform water flow feature analysis to obtain water flow feature grids; performing grid dynamic optimization based on the water flow feature grids to obtain optimized grid data; and performing boundary condition exchange on the pipe network model, river model, and surface model based on the optimized grid data to obtain a one-two dimensional coupled hydrodynamic model.

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

[0144] S3c. Read the grid inundation water depth (H(i, t)) and the grid inundation time (T(i, t)), conduct data verification on the calculation results, eliminate outliers to obtain the verified inundation water depth, conduct spatio-temporal continuity analysis on the verified inundation water depth, and perform data smoothing processing to finally obtain the corrected grid inundation data.

[0145] This embodiment adopts the adaptive grid technology to dynamically optimize the grid distribution according to the terrain features and flow features, improving the calculation efficiency and accuracy. This embodiment can support flood analysis and assessment at different scales, helping to comprehensively understand and evaluate the actual situation of the study area and providing strong support for scientific decision-making.

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

[0147] S4a. Read the calibrated grid inundation data and the grid cell tolerance values, calculate the inundation water depth ratio for each grid to obtain the water depth ratio data; calculate the inundation time ratio for each grid to obtain the time ratio data; calculate the area ratio for each grid to obtain the area ratio data; read the functional weight values of various underlying surfaces, and perform weighted calculations with the water depth ratio data, time ratio data, and area ratio data to obtain the grid performance values; perform summary processing on the grid performance value P(t) to obtain the system performance time series; extract the moment when the performance starts to decline from the system performance time series to obtain the decline moment data; extract the moment of the lowest performance point to obtain the lowest point moment data; extract the moment when the performance recovery ends to obtain the recovery moment data; integrate the decline moment data, the lowest point moment data, and the recovery moment data to obtain the system performance characteristic moment 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 inundation depth of the i-th grid at time t; HTi is the flood tolerance water depth threshold of the i-th grid; Ti(t) is the actual inundation time of the i-th grid at time t; TTi is the flood tolerance 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 new parameter Si is the spatial correlation weight and Ki is the key node influence factor.

[0148] S4b. Read the system performance time series and the 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 activation of the standby drainage facilities to obtain the redundancy value; read the lowest point moment and the performance recovery end moment in the system performance characteristic moment data, and calculate the recovery rate to obtain the recoverability value; combine the robustness value, the redundancy value, and the recoverability value to form a set of resilience characteristic values.

[0149] This embodiment proposes a method for extracting resilience characteristic values based on the system performance curve. It not only considers the system's ability to resist disturbances but also includes the system's recovery ability, making the resilience evaluation more comprehensive and scientific. By considering spatial correlation and key node influence, the spatial representativeness of the evaluation results is improved.

[0150] According to another aspect of the present application, step S41 can also be: reading the corrected grid inundation data and grid cell tolerance values, reorganizing the data according to the 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 functional weight values of various underlying surfaces, combining the water depth ratio data, time ratio data and area ratio data to obtain performance feature data; inputting the performance feature data into the trained LSTM model for prediction to obtain a predicted performance sequence; performing time step resampling on the predicted performance sequence 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 of the lowest performance point to obtain lowest point moment data; extracting the moment when the performance recovery ends to obtain recovery moment data; integrating the decline moment data, lowest point moment data and recovery moment data to obtain system performance feature moment data.

[0151] Through accurate system performance prediction and detailed performance feature moment data in this embodiment, the weak links of the system in flood disasters can be identified, and corresponding improvement measures can be proposed, which helps to improve the overall flood control and disaster reduction ability and reduce the losses caused by flood disasters.

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

[0153] S5a. Reading historical case data to construct an expert scoring matrix, calculating the subjective weight using the analytic hierarchy process, reading the set of resilience characteristic values to calculate the entropy value to obtain the objective weight, and performing arithmetic mean on the subjective weight and the objective weight to obtain the characteristic weight coefficients (A1, A2, A3).

[0154] S5b. Reading the set of resilience characteristic values and the characteristic weight coefficients, calculating the weighted sum ResG = A1×R1 + A2×R2 + A3×R3 for each grid, where ResG is the grid comprehensive resilience value, R1 is the robustness value; R2 is the redundancy value; R3 is the recoverability value; A1, A2, A3 are the characteristic weight coefficients; to obtain the grid resilience value.

[0155] S5c. Reading the grid resilience value, the functional weight values of various underlying surfaces and the preset time period weight coefficients, performing weighted average on the grid resilience values within each sub - basin to obtain the sub - basin resilience value, and performing weighted average on all sub - basin resilience values to obtain the overall basin resilience value.

[0156] According to another aspect of the present application, step S51 may 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 the toughness characteristic value set, calculating the membership function to obtain membership data; performing entropy value calculation on the membership data to obtain an objective weight value; reading the fuzzy relationship matrix, calculating the subjective weight value using the fuzzy Delphi method; performing fuzzy weighted average operation on the subjective weight value and the objective weight value to obtain the characteristic weight coefficient.

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

[0158] Step S1: Collect remote sensing image data of the watershed area, conduct a preliminary classification of the watershed underlying surface to obtain paddy fields, dry land, forest land, grassland, urban land, and unused land, and conduct a secondary classification based on POI data to obtain more refined 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: Obtain 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), perform radiometric correction, geometric correction, and atmospheric correction on the remote sensing images to ensure the high precision and consistency of the images;

[0160] Step S12: Use image processing techniques (such as K-means clustering or SVM (support vector machine) classification algorithm) to conduct a preliminary classification of the ground objects in the image, and classify them into major categories such as paddy fields, dry land, forest land, grassland, and urban land;

[0161] Step S13: Conduct a secondary classification of the urban area, use a deep learning model (such as a convolutional neural network (CNN)) to extract more refined spatial features. Identify buildings, roads, green land, etc. through training the model, and generate a preliminary classification result of the urban area, including the general urban land use types;

[0162] Step S14: Combine POI data for secondary classification and functional analysis:

[0163] Step S141: Obtain POI (Point of Interest) data, covering the coordinate and category information of functional areas such as commerce, residence, industry, public facilities, transportation, culture and entertainment, etc. POI data can usually be obtained through channels such as OpenStreetMap (OSM), Baidu Map, and Gaode Map. Clean the POI data, remove duplicate data, invalid data, etc., and classify it according to categories (such as commerce, residence, transportation, etc.).

[0164] Step S142: Conduct spatial overlay analysis on the POI data and remote sensing image data. Use GIS software (such as ArcGIS, QGIS, etc.) to further classify the identified urban areas in the remote sensing image according to the location of the POI data. For example, all POI points marked as "commerce" or "shopping center" may correspond to commercial land in the remote sensing image. For areas in the city that cannot be directly judged by the image (such as the commercial and residential mixed area in the suburbs), the function of the area can be judged through the surrounding POI data to ensure more accurate classification results.

[0165] Step S143: Further subdivide categories such as commercial land, residential land, and industrial land within the initially identified urban areas (such as urban land). Combine the function information of the POI, especially for types such as warehousing and logistics land, culture and entertainment land, and administrative and business land, and determine the specific uses of the corresponding areas through spatial matching and function analysis.

[0166] Step S15: Combine auxiliary data such as urban planning maps and land use planning maps. These data usually mark the detailed land use types of the city, such as commercial areas, residential areas, industrial areas, etc., and have strong authority and accuracy. Use these planning data to verify and optimize the classification results of the remote sensing images.

[0167] Step S16: Use social and economic data, such as population density, traffic flow, economic activity index, etc., as auxiliary variables to optimize the classification model. For example, areas with high population density are usually residential land or commercial land, while low-density areas may be industrial land or green spaces.

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

[0169] Step S171: Use ground truth data, existing land use maps or relevant statistical data as references to conduct classification accuracy evaluation. Verify the accuracy of the classification results by calculating classification accuracy indicators (such as user accuracy, producer accuracy, Kappa coefficient, etc.).

[0170] Step S172: Post-process the classification results to eliminate isolated pixels and misclassified areas, ensuring 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 function to social and economic activities, and calculate the functional scores of each underlying surface type.

[0172] Step S3: Comprehensively consider the main disaster-bearing bodies within each underlying surface type, and determine the tolerance levels of each type of disaster-bearing body in flood disaster events through historical disaster events, including flood-tolerant water depth and flood-tolerant time, and calculate the tolerance capacity of each grid cell in the basin.

[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, conduct flood disaster simulations, and calculate the flood depth and flood time of each grid cell in the basin area.

[0174] Step S5: Construct a method for evaluating the flood resilience of a basin considering the functions and tolerance capacities of the underlying surface. Based on the flood depth and flood time, the system performance curve, and the tolerance capacity of each grid, calculate the robustness, redundancy, and recoverability of each grid, and further determine the weights of each characteristic using subjective and objective weights, and calculate the flood resilience of the grid cell (A1*R1 + A2*R2 + A3*R3); (the robustness R1 is equal to the lowest value of the performance, the redundancy R2 is equal to the difference between the lowest values of the performance before and after the use of the standby drainage capacity, and the recoverability R3 is equal to (1 - the lowest value of the performance) / (the system performance at the end of the simulation - the lowest value of the performance)).

[0175] Step S6: Further, when statistically calculating the flood resilience of sub-basin units or the entire basin, it is necessary to systematically consider the functional weights of each underlying surface and consider the different degrees of attention to different underlying surface functions at different times, and dynamically adjust the weights of each underlying surface to obtain the resilience of the basin at different times, such as weekdays and non-weekdays, day and night, different growth stages of crops, etc.

[0176] In one embodiment of the present application, the system performance curve is a classical theory for evaluating system resilience and has also been widely applied in the field of flood disasters. The system performance curve reflects the change in system performance when the system is impacted by 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, which starts to decline from ts after an 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 flood resilience of the system is quantified using the system performance curve. The change process of flood resilience under one extreme rainfall event is as follows: The initial performance of the system is 1, which starts to decline from ts after the extreme rainfall event, reaches the minimum value of p(t) at tps, and gradually recovers until it fully recovers 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 flood severity, t n is the simulation duration, and p(t) is the system performance value at time t; approximating the flood severity Sev as a rectangular area, the formula for calculating 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 accumulated water time, t n is the total simulation time, and Res0 is the system resilience.

[0178] Existing solutions provide a fixed value for the elasticity of the system and fail to reflect the spatio-temporal variation of elasticity. This embodiment conducts a more refined assessment of the spatio-temporal process of elasticity, 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; Res w Res(i, t) is the flood resilience of the i-th grid point at time t; Sev w Sev(i, t) is the flood severity of the i-th grid point at time t; T is the total simulation time; N is the number of grids.

[0179] When calculating the resilience index of the grid system, many existing studies use the water depth threshold as the judgment 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, the system performance of the grid begins to be affected when encountering a relatively shallow inundation depth. For example, a water depth of 2 cm on an urban road can cause the risk of motor vehicle skidding, and relatively shallow water accumulation (such as 5 cm) on a pedestrian passage can affect pedestrians. At the same time, existing studies do not consider the impact of water accumulation time. However, for different land use types, due to their different functions, the requirements for water accumulation time also vary. For example, in urban areas, the lower threshold of the maximum allowable water recession time is 1 h, and the upper threshold is 24 h; for cultivated land, according to the types of crops, the threshold of the flood tolerance time is about 2 - 6 days, etc. Therefore, in the traditional resilience calculation method, the method considering a single threshold will cause a sudden change in system performance at the water depth threshold, and due to the lack of consideration of the flood tolerance time, it will lead to an overestimation of resilience and cannot reflect the true performance of the system. In view of this, this embodiment proposes two thresholds based on the flood tolerance depth and flood tolerance time to quantify resilience, namely the upper threshold and the lower threshold. The improved formula for calculating the system resilience is as follows:

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

[0181]

[0182] where 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 inundation time of the i-th grid point at time t; H(i, t) is the inundation depth of the i-th grid point at time t; H min (i) is the upper threshold of the inundation depth of the i-th grid point; H max (i) is the lower threshold of the inundation depth of the i-th grid; T(i, T) is the inundation time of the i-th grid point at time T; T min (i) is the lower threshold of the water entry time of the i-th grid; T max (i) is the upper threshold of the inundation time of the i-th grid.

[0183] For the basin scale, the land use types are divided into five categories: urban land, paddy fields, dry land, grassland, and forest land. The flood tolerance depths and flood tolerance times of each land use type are determined respectively, and the results are as follows:

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

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

[0186] This embodiment considers the functions of different underlying surfaces, calculates the function weights, and statistically analyzes the basin resilience. Compared with the conventional mean statistics, it is more reasonable; it considers the tolerance capabilities of different underlying surfaces and introduces the 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, considering different characteristics.

[0187] According to one aspect of the present application, a grid-scale basin flood resilience evaluation system based on the functions and tolerance capabilities of the underlying surface 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 executable by the processor, and the instructions are used to be executed by the processor to implement the grid-scale basin flood resilience evaluation method according to any one of the above embodiments.

[0191] In response to the problem of the recognition accuracy of the underlying surface function, this application adopts two-stage classification and POI-assisted recognition. The support vector machine is used for primary classification to divide land use into large categories such as paddy fields and dry land; then for urban land use, a deep learning model is used for secondary fine classification. In particular, through data augmentation and normalization processing, the generalization ability of the deep learning model is improved. The POI data is introduced to assist recognition. Through spatial overlay analysis and density calculation, combined with the preset density threshold standard, the accurate recognition of urban functions such as commerce and residence is realized, effectively solving the classification problem of functional mixed areas.

[0192] In response to the problem of quantitative assessment of tolerance ability, the survival analysis theory is introduced. By constructing a Cox proportional hazards model, the survival probability under different inundation conditions is calculated based on historical disaster data. Then, through the probability threshold curve and inflection point analysis method, the flood depth and time threshold of different disaster-bearing bodies are objectively determined, overcoming the subjectivity of traditional empirical judgment, and establishing a scientific assessment system based on data-driven.

[0193] In response to the problem of model coupling calculation, the calculation stability problem is solved through a coupling strategy. Specifically, by constructing a node mapping relationship, establishing a water volume exchange equation, and performing physical constraint processing, the rationality of the exchange flux is ensured. At the same time, by constructing a sparse matrix and adopting an interactive iteration format, the calculation efficiency is improved. In particular, through case verification and parameter optimization, the calculation stability of the coupling model is further guaranteed.

[0194] In response to the problem of identifying characteristic points of the performance curve, a time series feature recognition method based on multiple mathematical methods is developed. The time series is reconstructed through wavelet transform, and key nodes are identified using wavelet coefficients; the dynamic time warping model is used for pattern matching, and the lowest point of performance is determined by combining confidence interval analysis; the change point detection algorithm is used to identify the recovery inflection point, improving the accuracy and reliability of key time node identification.

[0195] In response to the problem of multi-scale weight transfer, an objective-subjective combined weight determination method is established. The objective weight is calculated by the entropy method, the subjective weight is determined by combining the fuzzy analytic hierarchy process, and the final characteristic weight coefficient is obtained by fuzzy weighted average. In the scale conversion process, the underlying surface function weight and time period weight coefficient are considered, realizing scientific weight transfer and ensuring the consistency of evaluation results at different scales.

[0196] In response to the problem of expressing time dynamics, the time period weight coefficient is introduced. By combining the time period basic data and time information data, the urban time period weight values on weekdays and non-weekdays, as well as the ecological time period weight values in different growth periods, are calculated respectively. Through the dynamic adjustment of the time period weight coefficient, the time-varying characteristics of the importance of the underlying surface function are expressed, improving the practicality of the evaluation results.

[0197] The present invention integrates multiple technical fields such as remote sensing image processing, POI data analysis, hydrodynamic simulation, and resilience evaluation, and constructs a complete basin flood resilience evaluation system. First, precise identification of underlying surface functions is achieved through multi-source data fusion and multi-level classification. Then, an assessment system for tolerance capacity is established based on historical disaster data, and the flood process is simulated through a one-two dimensional coupled hydrodynamic model. Finally, resilience characteristics are extracted based on the system performance curve and multi-scale evaluation is realized. This technical route of "function identification - disaster-bearing capacity - inundation simulation - resilience evaluation" not only considers the differences and spatio-temporal variation characteristics of regional functions, but also includes a comprehensive assessment of the system's resistance and recovery capabilities, improving the scientificity and reliability of the evaluation results. Especially in the aspect 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 deficiencies of traditional evaluation methods in terms of spatial accuracy and temporal dynamics are overcome, providing more accurate technical support for urban flood control and disaster reduction decision-making.

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

Claims

1. Grid-scale basin flood resilience assessment method based on underlying surface function and tolerance ability, characterized in that, It includes the following steps: S1. Obtain high-resolution remote sensing image data and POI data, and obtain the fine classification results through multi-level classification processing; combine the pre-stored socio-economic data to calculate the functional weight values of various underlying surfaces; S2. Conduct statistical analysis and grid calculation on the fine classification results and the pre-stored historical disaster data to obtain the grid unit tolerance ability values; S3. Perform simulation calculations based on the pre-configured coupled hydrodynamic model to obtain grid inundation data, including grid inundation water depth and inundation time; S4. Calculate the system performance curve based on the grid inundation data and the grid unit tolerance ability values, extract the resilience characteristic values, and obtain the set of resilience characteristic values for each grid; S5. Calculate the comprehensive resilience based on the set of resilience characteristic values; conduct spatial integration in combination with the functional weight values of various underlying surfaces to obtain the resilience evaluation results at different spatio-temporal scales; Step S1 is further as follows: S11. Read the high-resolution remote sensing image data, and successively perform radiometric correction processing, geometric correction processing, and atmospheric correction processing to obtain the corrected remote sensing image; S12. Based on the corrected remote sensing image, construct a primary classification model and conduct classification to obtain the primary classification results; read the urban land data therein, perform block division, data augmentation, standardization, and deep learning model training to obtain a fine classification model; accordingly, conduct secondary classification on the urban land in the primary classification results to obtain the fine classification results; S13. Read the POI data, conduct coordinate repeatability checks, delete duplicate points, functional classification annotation, and spatial distribution analysis to obtain the POI density data; S14. Based on the fine classification results and the POI density data, conduct data overlay, POI point density calculation, dominant function type calculation, and urban land division to obtain the fine classification results; S15. Based on the fine classification results, construct an analytic hierarchy process judgment matrix and conduct consistency test and standardization processing to obtain the functional weight values of various underlying surfaces; Step S3 is further as follows: S31. Read the pipe network GIS data and construct a one-dimensional pipe network model; read the river data and construct a one-dimensional river model; read the surface cover data and construct a two-dimensional surface model; based on the topographic data of the research area, couple the one-dimensional pipe network model, the one-dimensional river model, and the two-dimensional surface model to generate a one-two dimensional coupled hydrodynamic model; S32. Based on the one-two dimensional coupled hydrodynamic model, calculate the step size data and the interval data; combine the pre-stored design rainfall data and the fine classification results to conduct model calculations to obtain the grid inundation water depth data and the grid inundation time data; S33. Based on the grid inundation water depth data and the grid inundation time data, conduct outlier detection, eliminate abnormal records, spatio-temporal continuity analysis, and data smoothing processing to obtain the corrected grid inundation data; Step S4 is further as follows: S41. Calculate the system performance based on the grid inundation data and the grid unit tolerance ability values, extract the moment when the performance starts to decline, the moment of the lowest performance point, and the moment when the performance recovery ends to obtain the system performance characteristic moment data; Wherein the value range of the system performance is 0-1; S42. Based on the moment data of system performance characteristics, extract the minimum performance value, calculate the performance difference and recovery value before and after the activation of standby drainage facilities, and combine them to form a set of resilience characteristic values.

2. The grid-scale watershed flood resilience evaluation method based on underlying surface function and tolerance according to claim 1, characterized in that Step S2 is further as follows: S21. Based on the fine classification results and pre-stored historical disaster data, conduct disaster impact statistics and construct a correspondence table between land use types and disaster-affected objects to generate a list of disaster-bearing body types. S22. Based on the list of disaster-bearing body types and pre-stored historical disaster records, conduct data preprocessing and construct a survival analysis model, calculate the probability values affected by water depth and time, and obtain the threshold of flood tolerance depth and the threshold of flood tolerance time. S23. Divide the study area into grid cells of equal size, spatially divide the fine classification results accordingly, configure the corresponding thresholds of flood tolerance depth and flood tolerance time, and conduct comprehensive calculations to obtain the tolerance capacity values of the grid cells.

3. The grid-scale basin flood resilience evaluation method based on the underlying surface function and tolerance according to claim 2, wherein Step S5 is further as follows: S51. Read the pre-stored historical case data, construct a fuzzy relation matrix, calculate the membership degree and subjective weight, and obtain the subjective weight value; based on the set of resilience characteristic values, calculate the entropy value to obtain the objective weight value. Perform an arithmetic mean operation on the subjective weight value and the objective weight value to obtain the characteristic weight coefficient. S52. Based on the characteristic weight coefficient and the set of resilience characteristic values, perform weighted calculations for each grid to obtain the grid resilience value. S53. Based on the functional weight values of various underlying surfaces and the pre-stored time period weight coefficients, perform weighted average calculations on the grid resilience values within each sub-basin to obtain the sub-basin resilience value; perform weighted average calculations on all sub-basin resilience values to obtain the overall resilience value of the basin, which is the resilience evaluation result.

4. The grid-scale basin flood resilience evaluation method based on the underlying surface function and tolerance according to claim 3, characterized in that Step S11 is further as follows: S111. Read the high-resolution remote sensing image data, extract the image metadata, conduct integrity checks and data complementation processing to obtain complete image parameters. S112. Based on the complete image parameters, calculate the solar zenith angle and azimuth angle; combine the pre-stored sensor calibration parameters to conduct noise estimation, selective denoising and angle correction to obtain the radiometrically corrected image. S113. Based on the radiometrically corrected image and the pre-stored ground control point data, conduct feature matching, geometric transformation parameter calculation and error evaluation to obtain the transformation accuracy data; judge whether it meets the preset threshold requirements. If not, re-conduct feature matching; if it meets, perform resampling processing on the radiometrically corrected image to obtain the geometrically corrected image. S114. Read the geometrically corrected image and the real-time obtained meteorological observation data, calculate the atmospheric transmittance and atmospheric scattered irradiance, conduct dark pixel extraction, construct an atmospheric correction model, and obtain the corrected remote sensing image.

5. The grid-scale basin flood resilience evaluation method based on underlying surface function and tolerance according to claim 3, characterized in that Step S31 is further as follows: S311. Read the pipe network GIS data, conduct topological relationship checks, breakpoint detection, connectivity repair, attribute integrity checks, parameter supplementation and verification to obtain complete pipe network data. S312. Based on the complete pipe network data, construct the Saint-Venant equations, conduct discretization processing and solve them to obtain the pipe network characteristic values. Based on the pipe network characteristic values, construct a discrete format of the finite volume method, conduct solution verification to obtain the pipe network model. S313. Read the pre-stored river cross-section data and river plane data, perform cross-section interpolation calculation and curve fitting to obtain the river cross-section curve; read the pre-stored river roughness data, perform segmented assignment to obtain the cross-section roughness data; construct a river model based on the river cross-section curve and the cross-section roughness data; S314. Read the pre-stored DEM data, perform depression filling, slope analysis and roughness calculation to construct a surface model; S315. Based on the pipe network model, the river model and the surface model, construct a water volume exchange equation, perform physical constraint processing to obtain boundary condition data; S316. Based on the pipe network model, the river model, the surface model and the boundary condition data, construct a model coupling matrix, perform sparsification, interactive iteration format, example verification and optimization of solution parameters processing to obtain a one-two dimensional coupled hydrodynamic model.

6. The grid-scale watershed flood resilience evaluation method based on the underlying surface function and tolerance according to claim 3, wherein Step S41 is further as follows: S411. Read the grid inundation data and the grid cell tolerance values, calculate the inundation depth ratio, construct a spatial weight matrix to obtain the weighted water depth data; S412. Read the grid inundation data, calculate the inundation time ratio, combine the pre-stored time weight standard for hierarchical weighting and weighted calculation to obtain the weighted time data; S413. Read the weighted water depth data and the weighted time data, perform area normalization processing and multi-dimensional weighted calculation to obtain the grid performance value; S414. Read the grid performance value, perform time series reconstruction, wavelet transform, key time node identification and clustering analysis to obtain the performance time series data and the decline time data; S415. Read the performance time series data, construct a dynamic time warping model, identify the lowest point of performance and perform confidence interval analysis to obtain the lowest point time data; S416. Read the performance time series data, identify the performance recovery inflection point and determine 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.

7. Grid-scale basin flood resilience assessment system based on underlying surface function and tolerance, characterized in that Including: At least one processor; And, A memory communicatively connected to the at least one processor; wherein, The memory stores instructions executable by the processor, and the instructions are used to be executed by the processor to implement the grid-scale basin flood resilience evaluation method according to any one of claims 1 to 6 based on the underlying surface function and tolerance.

Images