Urban-scale rainfall infiltration and vegetation configuration collaborative optimization evaluation method

Abstract

The application discloses a kind of urban scale precipitation accommodation and vegetation configuration collaborative optimization evaluation method, belong to urban rain flood management optimization field, including: the basic data of city area are integrated and rainfall event identification;Through building energy consistency, wet surface dissipation and multiple physical gate such as geometric radiation sensitivity, physical mechanism inference is carried out, and the hydrological response characteristics of underlying surface are obtained;Combining the hydrological response characteristics of underlying surface with preset planning target, a robust optimization model is constructed and solved to determine the optimal vegetation configuration scheme;The optimal vegetation configuration scheme is back substitution and closed loop checking and effect evaluation.The application can accurately quantify the hydrological response characteristics of urban underlying surface, and under uncertain conditions, a vegetation configuration scheme with clear probability of achievement and controllable risk is developed, improving the scientificity and reliability of urban water resource resilience planning.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of urban rain flood management optimization, and particularly relates to a method for evaluating and optimizing the cooperation of urban-scale precipitation storage and vegetation configuration. BACKGROUND

[0002] As a region with high concentration of human activities, the hydrological process of city is profoundly changed compared with natural underlying surface. The high-density impervious surface blocks the natural infiltration of precipitation, leading to the rapid collection of surface runoff and increasing the risk of urban waterlogging. At the same time, the urban heat island effect is becoming increasingly significant. Therefore, how to ensure the safety of urban flood control, make full use of rainwater resources, and alleviate the heat island effect through the evapotranspiration of vegetation to achieve the coordinated management of urban water-heat balance has become a key scientific problem for sustainable development of city. As a core component of urban ecosystem, vegetation plays a crucial role in precipitation interception, soil water retention, and evapotranspiration water consumption. Scientific planning and configuration of urban vegetation can not only improve the ability of city to store and absorb precipitation, i.e. to store water during rainfall and restore water storage capacity through evapotranspiration during drought, but also can adjust the local climate through evaporation cooling effect. Therefore, developing a method for accurately evaluating and optimizing the cooperation of urban precipitation storage and vegetation configuration has great theoretical and practical significance for improving the climate resilience of city and building a sponge city.

[0003] Currently, the quantitative research on the overall water storage capacity of city (i.e. the ability to accommodate and store precipitation) is the basis for understanding the hydrological function of city. One of the research paths is based on the evaporation recession theory, which estimates the dynamic total water available for evaporation on the ground, i.e. the water storage capacity of city, by analyzing the decay process of total evapotranspiration (ET) on the ground over time during the dry period after rainfall stops. Specifically, researchers usually use hydro-meteorological models (such as the four-source evapotranspiration model) to simulate the daily ET time series of urban areas. Then, by identifying the dry period events with no rainfall, and assuming that the evapotranspiration process during the entire dry period follows a single exponential recession law, the evapotranspiration value ET0 of the initial day of the dry period and a macroscopic recession time scale λ are fitted. Finally, the two parameters are multiplied (S = λ·ET0), which can estimate the water storage capacity S of city as a whole. In order to ensure the stability of the estimation, a series of empirical constraints are imposed on the selected dry period events, such as limiting the length of dry period, excluding low-temperature snowfall conditions, and statistically averaging the results of multiple events over multiple years. This method overcomes the difficulty of direct measurement caused by the high heterogeneity of city surface, and provides a feasible technical approach for quantifying the water storage capacity of city at regional scale.

[0004] However, the existing technology based on the above macro evaporation recession model still has limitations in the fine characterization of physical processes, dynamic attribution of key hydrological processes, and robustness of decision-making, which is difficult to meet the needs of modern urban fine and forward-looking planning. These limitations mainly manifest in: the black box processing makes it impossible to accurately depict the real physical dynamics of urban dryness after rain; the lack of internal physical process decomposition makes it difficult to make internal mechanism attribution to observed phenomena, making the configuration of vegetation scheme unsatisfactory. SUMMARY

[0005] The present application provides a method for evaluating the optimization of urban-scale rainfall assimilation and vegetation configuration, to solve the above problems existing in the prior art.

[0006] The technical solution is a method for evaluating the optimization of urban-scale rainfall assimilation and vegetation configuration, comprising:

[0007] Obtain the basic data of the urban area, integrate the data and identify the rainfall events to generate a rainfall event data set;

[0008] Based on the rainfall event data set, perform physical mechanism inference to quantitatively obtain the hydrological response characteristics of the underlying surface;

[0009] Combine the hydrological response characteristics of the underlying surface with the preset planning target to construct and solve a robust optimization model, and determine the optimal vegetation configuration scheme;

[0010] Substitute the optimal vegetation configuration scheme back into the rainfall event data set to perform closed-loop checking and effect evaluation, and generate an optimization evaluation report.

[0011] The present application can accurately quantify the hydrological response characteristics of the urban underlying surface, and under uncertain conditions, develop a vegetation configuration scheme with clear probability of achievement and controllable risk, thereby improving the scientificity and reliability of urban water resource resilience planning. BRIEF DESCRIPTION OF DRAWINGS

[0012] Figure 1 A step flowchart of a method for evaluating the optimization of urban-scale rainfall assimilation and vegetation configuration is provided for the embodiments of the present application.

[0013] Figure 2 A step flowchart of quantitatively obtaining the hydrological response characteristics of the underlying surface is provided for the embodiments of the present application.

[0014] Figure 3 A step flowchart of constructing a physical gate and performing joint verification is provided for the embodiments of the present application.

[0015] Figure 4 A step flowchart of determining the optimal vegetation configuration scheme is provided for the embodiments of the present application. DETAILED DESCRIPTION

[0016] In order to better understand the technical scheme of the present application, the technical scheme in the embodiments of the present application will be described clearly and completely below in conjunction with the accompanying drawings of the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative effort should fall within the protection scope of the present application.

[0017] It should be noted that the terms “comprising” and “having” and any variations thereof are intended to cover not exclusively containing, for example, a process, method, system, product or device containing a series of steps or units does not have to be limited to those steps or units clearly listed, but can include other steps or units that are not clearly listed or inherent to these processes, methods, products or devices.

[0018] In the research, it is found that the existing technology simplifies the extremely complex evapotranspiration process after rain into a single exponential decay, which essentially regards the urban underlying surface (including vegetation canopy, soil, impervious surface, etc.) as a linear reservoir with uniform physical properties. This ignores the huge differences in water release mechanisms and rates of different water storage media, especially the inability to distinguish between the extremely fast evaporation process dominated by vegetation canopy interception water and surface thin water film in the early stage after rain, and the slower transpiration and evaporation process limited by soil water diffusion in the later stage. This black box treatment leads to its inability to accurately depict the real physical dynamics of urban drying after rain. In addition, due to the lack of decomposition of internal physical processes, it is difficult to make internal mechanism attribution to the observed phenomena. For example, research has found that increasing urban vegetation sometimes leads to a decrease in estimated total water storage capacity S, which the model itself cannot explain, and must rely on other models that can decompose evapotranspiration components for external and supplementary analysis to find the reason (such as the shading effect of vegetation increase on impervious surface evaporation). The lack of such diagnostic capabilities makes it difficult to guide how to configure vegetation to achieve the goal of optimal water storage. By extension, based on such macro and mechanism-fuzzy models for optimization decision-making, the reliability and robustness of the scheme are questionable. Since the model cannot accurately capture the nonlinear impact of vegetation configuration changes on hydrological physical processes, and does not systematically consider future meteorological conditions and other deep uncertainties, the optimization results may not perform well in some unforeseen scenarios, lack of clear probability guarantee, and are difficult to meet the high reliability requirements of major municipal engineering planning.

[0019] As shown in Figure 1 , a method for co-optimizing and evaluating urban-scale rainfall assimilation and vegetation configuration is proposed, including the following steps:

[0020] Acquire basic data on urban areas, integrate the data and identify rainfall events to generate a rainfall event dataset.

[0021] In this embodiment, the basic data for the urban area includes: meteorological data, rainfall data, remote sensing imagery, and urban morphology data; specific data types include net radiation, surface temperature, leaf area index, etc. That is, the basic data covers, but is not limited to, multiple sources and types. For example, meteorological observation and reanalysis data may include variables such as net radiation, shortwave radiation, air temperature, wind speed, wind direction, and relative humidity; rainfall data, such as minute-level or hourly rainfall records; remote sensing data, which may include surface temperature (LST), leaf area index (LAI), surface albedo, impermeability, or land use type; and urban morphology data, such as digital surface models (DSM) describing building layout and height, street valley orientation, surface roughness length, and building density. To achieve data integration, the basic data from different sources are unified onto a common spatiotemporal reference. Spatially, a target resolution can be set, such as a 30m × 30m grid, and a geometrically overlapping area-weighted projection method is used to resample all spatial data to this target grid network. The calculation process of this method can be expressed as: V target (p) =Σ s∈overlap(p) w(p,s)·V source (s); where V target (p) represents the variable value on the target raster p; Σ is the sign of summation over the set; s is the index of the source raster that has a geometric intersection with the target raster p; overlap(p) is the set of all source rasters that intersect with the target raster p; w(p, s) is the area weighting coefficient, which is calculated as w(p, s) = area intersect (p, s) / area target (p); V source (s) represents the variable value on the source raster s; area intersect (p, s) represents the geometric intersection area of ​​the target raster p and the source raster s; area target (p) represents the area of ​​the target raster p; · represents the real number multiplication symbol. Temporally, it is unified to a fixed time step, such as 1 hour. For cumulative amounts (such as rainfall), resampling is performed using a summation method: A target (t) = Σ k∈cover(t) A source (k); where A target (t) represents the cumulative amount at the target time step t; k is the index of the original time step covered by the target time step t; cover(t) is the set of all original time steps covered by the target time step t; A source (k) represents the cumulative amount at the original time step k. For average quantities (such as temperature and wind speed), a time-weighted average is used: Mtarget (t) = [Σ k∈cover(t) M source (k)·dur(k)] / dur total (t); where M target (t) represents the average value at the target time step t; M source (k) represents the average value at the original time step k; dur(k) represents the duration of the original time step k; dur total (t) represents the total duration covered by the target time step t. This method resolves the inconsistency in the spatiotemporal dimensions of multi-source heterogeneous data, resulting in a preliminarily aligned dataset.

[0022] Furthermore, to ensure the physical consistency and integrity of the data, quality control and missing data repair are performed. In this process, this embodiment provides an adaptive missing data repair process. This process adaptively assigns weight coefficients based on the temporal stability and spatial continuity scores of the variable to be repaired, and uses these weight coefficients to perform a convex combination of the temporal interpolation result and the spatial neighborhood weighted mean to obtain the repaired value. Specifically, for a missing variable V(x, t) at location x and time t, its repaired value V... repaired (x, t) = w t ·V temporal (x, t) + w s ·V spatial (x, t); where V temporal (x, t) represents the result obtained through time interpolation (such as linear interpolation or spline interpolation); V spatial (x, t) is the result obtained by weighted averaging (e.g., inverse distance weighting) of the effective values ​​within its spatial neighborhood; w t with w s These are the time and space weighting coefficients, respectively, and w t +w s = 1, w t w s ≥0; The determination of these two weighting coefficients is based on the analysis of the historical data of the variable: if the variable changes slowly over time (strong time autocorrelation), then w is assigned a weight of ≥0. t A higher value is assigned to w if the variable is spatially continuous (strong spatial autocorrelation). s A relatively high value. This adaptive mechanism allows the repair process to fully utilize the spatiotemporal characteristics of the variables themselves, improving the accuracy of the repair results.

[0023] On this basis, rainfall event identification is performed. Instrument noise is filtered out by setting a minimum distinguishable rainfall threshold (e.g. 0.1 mm / hour), and continuous rainfall time series are cut into a series of independent rainfall-dry period events according to a rolling cumulative rainfall threshold (e.g. cumulative rainfall exceeding 5 mm) and a minimum rain-free interval (e.g. 6 consecutive hours without effective rainfall). The start and end times of each event, the rainfall cessation time, and the evaluation window for subsequent analysis are determined. All processed variables are organized according to grid location, time, and variable name to form an eventized data cube, and metadata are attached to form the rainfall event dataset. This dataset provides a unified, clean, and traceable data entry for all subsequent steps. That is, the rainfall event dataset includes an eventized data cube and an event list. It describes the data collection organized around rainfall events after preliminary processing and is the basis for all subsequent analysis.

[0024] Based on the rainfall event dataset, physical mechanism inference is performed to quantify the underlying surface hydrological response characteristics.

[0025] In this embodiment, after analyzing each rainfall event, the water evaporation process of the urban underlying surface is analyzed to reveal its ability to absorb rainfall. The underlying surface hydrological response characteristics are a collection of quantitative indicators, and their upper concepts can include but are not limited to: the effective water storage capacity of a single rainfall event (in millimeters), the time constant describing the speed of evaporation decay, the switching time from the wet surface dominant phase to the soil limited phase, and the uncertainty information of these characteristic quantities, which are the basis for subsequent optimization configuration. In other words, the underlying surface hydrological response characteristics include the upper generalization of a series of mechanism inference results such as event water storage capacity distribution and uncertainty information, and essentially are a quantitative description of the physical characteristics of how the model responds to rainfall (hydrology) on the ground (underlying surface).

[0026] The robust optimization model is constructed and solved by combining the underlying surface hydrological response characteristics with the preset planning objectives to determine the optimal vegetation configuration scheme.

[0027] Optionally, the planning target can be a specific amount of precipitation to be absorbed (e.g. requiring the urban green space system to absorb an additional 20 mm of precipitation per year), or other indicators related thereto. The robust optimization model is a mathematical model that can seek a reliable solution that can meet the planning target with a high probability under consideration of uncertainty (e.g. randomness of meteorological conditions, estimation error of model parameters). Exemplarily, the robust optimization model is a generalization of the distributionally robust chance-constrained model, seeking a reliable solution (robustness) under uncertainty. By solving the model, an optimal vegetation configuration scheme can be obtained, which is specific to each candidate plot and explicitly indicates the vegetation adjustment to be made, such as increasing the leaf area index of a plot from 1.5 to 2.5, or transforming grassland into shrubs.

[0028] The optimal vegetation configuration scheme is back-substituted into the rainfall event dataset for closed-loop checking and effect evaluation, to generate an optimization evaluation report.

[0029] In the present embodiment, verification and evaluation are performed. The underlying surface parameters (such as leaf area index, surface albedo, etc.) corresponding to the optimal vegetation configuration scheme are updated in the rainfall event dataset to form a new scenario after implementation of the simulation. The physical mechanism inference process is re-run to calculate the underlying surface hydrological response characteristics under the new scenario. By comparing the total amount of water storage capacity before and after implementation, the achievement rate of the scheme to the planning target can be evaluated, and other environmental impacts (such as changes to the urban thermal environment) that can be brought about are analyzed. All analysis results are summarized to generate an optimization evaluation report, providing a scientific basis for actual engineering decision-making.

[0030] In a possible implementation, the urban-scale precipitation assimilation and vegetation configuration co-optimization evaluation method can also be configured to: obtain meteorological data, rainfall data, remote sensing images, and urban form data, generate eventized data cubes and an event list through data integration and event identification; based on the eventized data cubes and the event list, determine the stage switching time from wet surface dominance to soil limitation using a statistical-physical joint identification method, and perform physical correction and parameter estimation to obtain event water storage capacity distribution and uncertainty information; combine the event water storage capacity distribution, uncertainty information, and preset precipitation assimilation targets and engineering constraints to build and solve a distribution robust chance constraint model coupled with physical processes, and output an optimal vegetation configuration scheme; update the eventized data cubes according to the optimal vegetation configuration scheme, and combine the event list to perform back evaluation through mechanism chain recalculation, and generate a technical report containing precipitation assimilation rate and risk assessment. Specifically, meteorological data, rainfall data, remote sensing images, and urban form data are input to complete resampling, missing data repair, and quality review on a unified time scale and spatial grid, so that data from different sources are consistent with each other; according to the thresholds of continuous rain-free period and cumulative rainfall, the rainfall-dry period event is divided, and the rain stop time and evaluation window of each event are determined; the aligned data set and event list are obtained as the only entry for subsequent calculation. With the aligned data set and event list as input, on each grid and each event, the stage switching time from wet surface dominance to soil limitation is first determined through statistical variable point search and joint determination of energy consistency, wet surface dissipation, and street valley radiation geometry; then, before the switching, the intercepted evaporation and impervious surface evaporation are consistently decomposed in the order of energy conservation, time decay, and shadow-sunlight sensitivity, and the time series of each evapotranspiration component is reconstructed; after the switching, a robust joint estimation method is used to obtain key parameters describing the dry period decay rate and reference dissipation intensity, and the effective water storage capacity of a single event scale is accumulated according to the parameters, while the uncertainty range is output; for impervious surface evaporation, the advection enhancement is identified and written back in a short window with clear sky, significant surface temperature gradient, and favorable wind direction, so that the whole process curve is physically self-consistent; finally, the stage switching time distribution, component time series, mechanism parameter set, event water storage capacity distribution, and uncertainty information are obtained.With the event water storage capacity distribution, mechanism parameter set, and uncertainty information as inputs, combined with the target of rainfall to be absorbed, candidate land set, budget, and land and thermal environment constraints, a number of representative leaf area index values are selected in each land, the anchor point data and marginal response are generated by calling the preceding mechanism inference chain one by one, and whether the stage switching occurs is recorded; on this basis, a semi-convex substitution relationship is constructed considering the physical shape constraint and solving efficiency, and a distribution robust chance constraint model coupled with the physical process is established, so that the solution still has a clear achievement probability under uncertain conditions; after solving, the optimal vegetation configuration scheme, expected water storage increment distribution, achievement probability, cost and risk assessment, and implementation list are output, which are used for decision-making and verification before implementation. With the optimal vegetation configuration scheme and the aligned data set, the event list as input, the configuration results are written back to the new leaf area index and geometric parameters, and the mechanism inference process is run again to obtain the evapotranspiration components and event water storage capacity under the implementation scenario; the total water storage capacity before and after the implementation is compared at the scale of the study area, the absorption rate relative to the target rainfall is calculated and the uncertainty range is given, and the changes of component redistribution and surface thermal environment are evaluated to verify whether the planning and safety constraints are met; finally, the water storage capacity distribution after implementation, the absorption rate and its confidence range, the component and thermal environment evaluation index are formed, and are summarized as a technical report and method metadata, providing a basis for subsequent urban or annual expansion and reuse.

[0031] As shown in Figure 2 According to one aspect of the present application, the quantification of the underlying surface hydrological response characteristics includes:

[0032] For the total evapotranspiration process in the rainfall event data set, the division point of the two-section exponential decay model with the optimal fitting is searched, and the candidate switching time is preliminarily identified.

[0033] Specifically, after the rain stops, the total evapotranspiration process usually goes through two main stages: the wet surface dominated stage dominated by interception water and surface thin water film evaporation, whose evapotranspiration rate decays quickly; and the soil limited stage limited by soil moisture supply, whose evapotranspiration rate decays slowly. Both stages can be approximately described by an exponential decay function. Therefore, in the evaluation window after the rain stops, each possible time point can be traversed as a candidate time tc of stage switching. For each tc, the time period t≤tc is fitted with the first exponential decay model, the time period t>tc is fitted with the second exponential decay model, and the total residual sum of squares (SSE) of the two sections is calculated. The mathematical expression of this process is: SSE (tc) =∑ t≤tc [ET obs (t) - E0·exp(-(t-t0) / τ w )] 2 + Σ t>tc [ET obs(t) - ET0·exp(-(t-tc) / τ s )] 2 ; wherein SSE(tc) is the total sum of squared residuals at the candidate switching time tc; ET obs (t) is the observed total evapotranspiration depth at time t in millimeters; E0 is the initial dissipation strength of the wet-surface dominated phase; exp(·) is the natural exponential function; t0 is the rain stop time; τ w is the decay time constant of the wet-surface dominated phase; ET0 is the reference dissipation strength of the soil-restricted phase; τ s is the decay time constant of the soil-restricted phase; 2 is the square operation. By searching for tc that minimizes SSE(tc), one or more statistically optimal candidate switching times can be identified preliminarily. This provides an objective, data-driven starting point for subsequent physical verification, avoiding the bias caused by relying entirely on subjective judgment or fixed thresholds.

[0034] Constructing physical gating and calling multi-source physical quantities in the rainfall event dataset, the candidate switching times are jointly verified; wherein the physical gating includes energy consistency, wet-surface dissipation, and geometric radiation sensitivity.

[0035] Exemplarily, energy consistency, wet-surface dissipation, and geometric radiation sensitivity are constructed as three physical gating, using physical laws to constrain and verify the candidate switching times of the purely statistical results.

[0036] Further, as shown in Figure 3 , constructing physical gating and jointly verifying further include:

[0037] For the energy flux after the candidate switching time, the linear regression relationship between the sum of latent heat and sensible heat flux and the difference between net radiation and surface heat flux is verified, requiring the slope of the linear regression relationship to be close to one and the intercept to be close to zero.

[0038] Specifically, within a short time window (e.g. 24 hours) after the candidate switching time tc, the latent heat flux LE(t) and the sensible heat flux H(t), as well as the net radiation Rn(t) and the ground heat flux G(t) are extracted. The ground energy balance equation is Rn(t) - G(t) = LE(t) + H(t). Therefore, a linear regression model can be constructed: (LE(t) + H(t)) = slope · (Rn(t) - G(t)) + intercept. A physically reasonable switching time means that the ground energy distribution mechanism after the soil limited stage should tend to be stable, i.e. the slope of the regression relationship should be close to 1, and the intercept should be close to 0. In addition, in a preferred implementation, the relative residual variance of the regression model is also required to have a significant drop compared to the period before the switching, indicating that the consistency of the energy closure has been improved. This gating makes the identified stage switching reasonable at the energy balance level.

[0039] A wet surface dissipation index is constructed by weighting the ground temperature rate of change, the ground albedo rate of change and the shortwave radiation recovery ratio, and it is determined whether the wet surface dissipation index exceeds a preset threshold in a continuous period.

[0040] In a possible implementation, the wet surface dissipation index is constructed by: extracting the ground temperature and the ground albedo time series from the rainfall event data set, and calculating the first-order time change rates of the ground temperature and the ground albedo respectively to obtain the ground temperature rate of change and the ground albedo rate of change; calculating the ratio of the current shortwave radiation to the maximum shortwave radiation on the same day as the shortwave radiation recovery ratio; and combining the ground temperature rate of change, the ground albedo rate of change and the shortwave radiation recovery ratio by weighting to construct the wet surface dissipation index. The index is used to determine whether the ground has changed from wet to dry. Its construction formula is: W(t) = a · dLST dt_norm (t) + b · dAlbedo dt_norm (t) + c · Rs ratio (t); where W(t) is the wet surface dissipation index at time t; dLST dt_norm (t) is the normalized first-order time change rate of the ground temperature, which usually accelerates when the ground changes from wet to dry; dAlbedo dt_norm (t) is the normalized first-order time change rate of the ground albedo, which usually rises when the ground changes from wet to dry; and Rs ratio(t) is the shortwave radiation recovery ratio, equal to the shortwave radiation Rs(t) at the current time divided by the possible maximum shortwave radiation of the day (which can be estimated by the clear sky model), reflecting the degree of weather clearing; a, b, c are weight coefficients, which can be calibrated by sensitivity analysis on historical data to reflect the actual influence degree of each component on the total evapotranspiration process. When the index W(t) continuously exceeds the preset threshold (for example, 0.8) for a period of time (for example, 3 hours) after the candidate switching time tc, it is considered that the signal of rapid dissipation of ground water is clear, and the verification is passed. This gating confirms the occurrence of stage switching from the perspective of changes in the physical properties of the ground.

[0041] In the transition period from shadow to sunshine, the response amplitude of interception evaporation and impervious surface evaporation to shortwave radiation is calculated respectively, and the response amplitude of interception evaporation is required to be not higher than that of impervious surface evaporation.

[0042] In this embodiment, this step shows the process of performing physical gating verification of geometric radiation sensitivity. In a possible implementation, according to the urban street valley direction and the solar elevation angle, the transition period of the ground surface from being not in sunshine to being in sunshine is identified; in each transition period, the finite difference response amplitudes of interception evaporation and impervious surface evaporation to shortwave radiation changes are calculated respectively; and the response amplitude of interception evaporation is required to be not higher than that of impervious surface evaporation. Specifically, in the wet surface dominant stage, evaporation mainly comes from interception water (existing on the surface of plant leaves) and thin water film on impervious surface. The interception water has a relatively flat response to light changes due to its three-dimensional structure; while the thin water film on the impervious surface (such as asphalt pavement) will rapidly enhance evaporation after receiving direct light. Therefore, the stage can be judged by comparing the sensitivities of the two to shortwave radiation. Specifically, according to the urban morphology data (street valley direction) and the solar position, the transition period of the ground surface from the shadow area to the sunshine area is calculated. In this period, the response amplitude R sens (q, t) = |ΔE q (t) / ΔRs(t)|; wherein R sens (q, t) is the finite difference sensitivity of component q to shortwave radiation at time t; q can take the value of interception evaporation component Ei or impervious surface evaporation component Eb; ΔE q (t) is the change of component q at the adjacent time step; and ΔRs(t) is the change of shortwave radiation at the adjacent time step. In the real wet surface dominant stage, R sens (Ei, t) ≤ R sens (Eb, t) should be met. If this relationship is no longer stable after the candidate switching time tc, the rationality of the switching time is supported. This gating uses the micro geometry of the urban underlying surface and the radiation transmission mechanism to provide a unique constraint for stage judgment.

[0043] The candidate switching time that passes the joint check is confirmed as the switching time from the wet surface dominant stage to the soil limited stage.

[0044] That is, only when a candidate switching time passes the above three physical gating checks of energy consistency, wet surface dissipation and geometric radiation sensitivity, it is finally confirmed as a reliable stage switching time t star The joint determination method combining statistics and multi-source physical information improves the accuracy and robustness of switching time identification.

[0045] Based on the stage switching time, the subsequent soil limited stage is parameterized and integrated to generate the underlying surface hydrological response characteristics.

[0046] Specifically, after the stage switching time t star is determined, the soil limited stage from t star to the end of the evaluation window is analyzed. Further, based on the stage switching time, parameter estimation and integration are performed, including: in the soil limited stage from the stage switching time, a joint objective function for determining the parameters of the exponential decay model is constructed; the joint objective function includes a robust residual term for measuring the goodness of model fitting, and a physical penalty term for quantifying energy closure consistency; by minimizing the joint objective function, the decay time constant and the reference dissipation intensity are solved simultaneously, and are used for the final generation of the underlying surface hydrological response characteristics. Specifically, the joint objective function can be represented as: J total (τ s , ET0) = Σ t>tstar ρ(ET obs (t) - ET0·exp(-(t-t star ) / τ s )) + β·Σ t>tstar |ε rel (t)|; wherein J total is the joint objective function; τ s and ET0 are the decay time constant and the reference dissipation intensity to be solved; ρ(·) is a robust loss function (such as Huber loss), which is less sensitive to outliers than square loss; β is the weight coefficient of the physical penalty term; ε rel (t) is the relative energy residual, equal to (LE(t)+H(t)) - (Rn(t)-G(t)) divided by the absolute value of (Rn(t)-G(t)), which is used to quantify the degree of energy closure. By minimizing J total , τ swith the optimal estimate of ET0. The data fitting is done in one step with the physical constraints fused in, avoiding physically unreasonable results that might be caused by simple fitting. The effective storage capacity S event for a single event is calculated by time integration of the fitted exponential decay model, with the formula: S event =∫ tstar tstar+T ET0·exp(-(t-t star ) / τ s )dt; where∫is the integral sign over the time interval; t star is the phase switching time; T is the length of the evaluation window; when the evaluation window is much larger than the time constant τ s , an approximation can be used: S event ≈τ s ·ET0; volume conversion is written as V event = S event ·A catch ; V event is the event volume; A catch is the catchment area. The final generated land surface hydrological response characteristics include a series of key parameters such as t star , τ s , ET0, S event , and their uncertainty ranges obtained through regression analysis.

[0047] In addition, as a further refinement of the internal processes of the wet surface dominant phase, after confirming the phase switching time, it also includes: within the time period from the rain stop to the phase switching time, an optimization decomposition problem is established with the objective of minimizing the difference between the sum of the total evapotranspiration and the multi-evapotranspiration component; and at least three types of physical constraints are imposed on the optimization decomposition problem: energy conservation constraint, limiting the total latent heat flux to not exceed the available energy; time sequence constraint, requiring the decay time constant of intercepted evaporation to be less than that of impermeable surface evaporation; in other words, the time sequence constraint specifically includes: requiring the decay time constant representing the rapid dissipation process of intercepted evaporation to be less than the decay time constant representing the slower dissipation process of the thin water film on the impermeable surface; geometric sensitivity constraint, limiting the response amplitude of intercepted evaporation to shortwave radiation to not exceed that of impermeable surface evaporation; solving the constrained optimization decomposition problem to obtain the time series of each evapotranspiration component, which is used for the generation of land surface hydrological response characteristics. Through the mathematical programming method with physical constraints, the total evapotranspiration is decomposed into more detailed physical components, providing more detailed process information for subsequent sensitivity analysis and side effect evaluation (such as changes in the ground thermal environment).

[0048] The embodiment solves the technical problems of unclear physical mechanism and inaccurate dynamic description caused by regarding the evapotranspiration process as a single black box model by decomposing the evapotranspiration process after rain into two stages with different physical mechanisms, namely, a wet surface dominated stage and a soil limited stage, and constructing a joint identification framework composed of statistical optimal search and three physical gates of energy, wet surface and geometric radiation. Specifically, by using a two-section exponential fitting residual minimization algorithm, an hourly total evapotranspiration time series is input, and a statistically optimal candidate switching time tc is output. The candidate time must pass three physical verifications: the energy consistency gate makes the stage switching comply with the thermodynamic law that the surface energy distribution tends to be stable; the wet surface dissipation index gate confirms the physical state change of the surface from wet to dry through remote sensing data such as surface temperature and albedo; and the geometric radiation sensitivity gate uses the unique street valley geometry and light relationship in cities to constrain the response characteristics of different evapotranspiration components. The synergistic effect of statistics and multi-physical domain information (thermodynamics, physical state and radiation geometry) makes the finally identified stage switching time t star has high physical authenticity. The vague and continuous recession process is converted into a structured two-stage problem with an explicit physical inflection point, so that accurate parameter estimation (such as the decay time constant τ s ) of the soil limited stage which really determines the slow water release capacity of the city in the later period is possible, thereby improving the accuracy and mechanism explainability of the evaluation of the water storage capacity of the city.

[0049] In another embodiment of the present application, quantifying the hydrological response characteristics of the underlying surface further includes:

[0050] Based on the shortwave radiation, surface temperature gradient and wind field data in the rainfall event data set, a trigger window of enhanced advection is identified, which meets the conditions of clear sky, temperature difference along the wind direction reaching a preset threshold and wind speed reaching the standard; only in the trigger window of enhanced advection, a non-negative enhancement term related to the product of the horizontal gradient of the surface temperature and the effective wind speed is quantified as the amount of enhanced evaporation of the impervious surface; the amount of enhanced evaporation is written back to the evaporation component of the impervious surface, the evapotranspiration process is corrected, and the hydrological response characteristics of the underlying surface are generated based on the corrected evapotranspiration process.

[0051] In this embodiment, advection refers to the energy transport process in a local area. For example, dry and hot air flowing from a surface (e.g., an asphalt parking lot) to a relatively humid surface (e.g., a thin water film on a stairway after rain) can bring additional energy, causing the evaporation rate of the latter to exceed the level supported by the net radiation it receives. This effect is particularly pronounced in urban environments, which are highly heterogeneous. To accurately identify and quantify this effect, it is necessary to determine the advection enhancement trigger window, which is a set of times, each of which must satisfy the following three conditions: clear sky condition: to ensure that there is enough energy to drive the temperature difference between surfaces. This can be achieved by determining the shortwave radiation value, i.e., Rs(t) > θ Rs ; where Rs(t) is the shortwave radiation at time t, in watts per square meter; θ Rs is the radiation threshold, which can be set to 400 watts per square meter. Significant temperature difference along the wind direction condition: to ensure the potential for energy transport; the ground temperature level gradient indicator GradTs alongwind (t) is calculated, and GradTs alongwind (t) > θ grad is required; where θ grad is the gradient threshold, which can be set to 0.05 Kelvin per meter. GradTs alongwind (t) is calculated by first obtaining the spatial gradient vector of the ground temperature ▽LST(x, y, t), and then taking the dot product of it with the wind direction unit vector v wind (t) at time t. Wind speed meets the condition: to ensure that energy can be effectively transported. The friction wind speed u star (t) can be used as an indicator, and u star (t) > θ ustar is required; where u star (t) is a physical quantity that describes the near-surface turbulence intensity, which better reflects the exchange efficiency of momentum and energy than the conventional wind speed; θ ustar is the friction wind speed threshold, which can be set to 0.15 meters per second. In the identified trigger window, the impervious surface evaporation component Eb(t) is corrected. The corrected impervious surface evaporation Eb corrected (t) is composed of the baseline part Eb base (t) and the non-negative advection enhancement term Eb adv (t). Eb base (t) is the impervious surface evaporation obtained by decomposition without considering the advection effect. The calculation model of the advection enhancement term is: Eb adv (t) = α · max(0, GradTs alongwind (t) · u eff (t) - θ trigger ); where Eb adv(t) is the advection-enhanced evaporation at time t, in millimeters; a is the advection enhancement proportionality coefficient to be estimated, which has the physical meaning of the evaporation increment caused by unit energy transport flux; max(0, ·) is the positive function to ensure the non-negativity of the enhancement term; u eff (t) is the effective wind speed, which can be the friction wind speed u star (t) or the near-surface wind speed after the crown valley resistance correction; θ trigger is the triggering threshold to prevent weak gradients or wind speed fluctuations from being mistakenly amplified as advection effects. After the parameters a and θ trigger are estimated, Eb adv (t) in all triggering windows can be calculated and added back to Eb base (t) to obtain the corrected Eb corrected (t). The correction is updated back to the total evapotranspiration, and the energy balance is rechecked to make the entire process physically self-consistent. Based on the corrected evapotranspiration process curve, subsequent parameter estimation and integration are re-performed to generate more accurate underlying surface hydrological response characteristics.

[0052] Further, the advection enhancement of impermeable surface evaporation is quantified, including: constructing a tool variable that attenuates and weights the temperature difference between the upwind and downwind near-surface temperatures along the wind path, to avoid endogenous bias; and performing a placebo test in the low shortwave radiation period to verify the invalidity of the tool variable under non-triggering conditions; using the verified tool variable to perform regression to robustly estimate the advection enhancement within the advection enhancement triggering window.

[0053] In this embodiment, GradTs alongwind (t) is directly used as the independent variable to regress Eb(t) to estimate the parameter a, which faces a serious endogenous problem. This is because the surface temperature gradient not only leads to evaporation enhancement, but the evaporation enhancement itself (through dissipation of latent heat) also affects the surface temperature in turn, forming a simultaneous equation bias. To solve this problem, the method of instrumental variable (IV) regression is introduced. Specifically, the tool variable IV temp (x, t) is constructed, which is related to the endogenous variable GradTs alongwind (x, t) in space, but not directly related to the random error term of the local evaporation process. A preferred construction is the weighted temperature difference between the upwind and downwind along the wind path. For the target grid x, at time t, determine the upwind and downwind grid sets. Then calculate: IV temp (x, t) = [Σ i∈Upwind(x) w(d i )·LST(i, t)] - [Σ j∈Downwind(x) w(d j) · LST(j, t)]; where Upwind(x) and Downwind(x) are the sets of up- and downwind grid cells of grid x; LST(i, t) is the land surface temperature of upwind grid i; LST(j, t) is the land surface temperature of downwind grid j; w(d) is a distance decay weighting function, e.g., w(d) = exp(-d / L char ), d is the distance of grid i or j from grid x, and L char is a characteristic decay length (e.g., 100 meters). This instrumental variable reflects the temperature field pattern at a larger scale, which drives the local temperature gradient, but whose value is not directly affected by the local evaporation, thus satisfying the instrumental variable requirement. Before performing the two-stage least squares regression with the instrumental variable, a test of the validity of the instrumental variable is also needed. One effective way is to perform a placebo test. Specifically, during non-trigger periods (e.g., nighttime or periods of very low shortwave radiation), the physical mechanism of advection enhancement does not exist. At these times, the instrumental variable regression is repeated. If the instrumental variable is valid, then the regression-derived coefficient a should not be significantly different from zero during these placebo periods. If the test passes, then the instrumental variable can be used to perform the formal regression on the data within the trigger window, resulting in a robust, unbiased estimate of the parameter a.

[0054] In an optional embodiment, further comprising: identifying time breakpoints where wind direction changes abruptly in the rainfall event dataset; and using each time breakpoint as a natural experiment to cross-validate the dynamic changes of the advection enhancement estimated by the instrumental variable regression.

[0055] In this embodiment, to further verify the reliability of the advection enhancement model, a cross-validation method based on natural experiments is introduced. The wind field in a city sometimes changes rapidly and significantly due to the passage of weather systems. This exogenous wind direction mutation can be regarded as an experimental treatment. Specifically, time breakpoints t break where the wind direction changes abruptly are identified in the time series, e.g., the average wind direction changes more than 60 degrees within one hour. For each t break , the dynamic changes of the advection enhancement term Eb adv (t) in a small time window (e.g., 3 hours) before and after t adv are compared. A correct model should predict Eb adv (t) that is consistent with the physical effect of the wind direction mutation. For example, if the wind direction changes from sweeping over a large park (cold source) to sweeping over a large commercial plaza (warm source), then for a wet impervious surface located downwind, Eb breakAfter that, there is a significant jump. By checking whether the dynamic changes predicted by the model are consistent with these physical expectations driven by wind direction mutations, the effectiveness of the model can be cross-verified from multiple angles and dynamically, further enhancing the credibility of the results.

[0056] Optionally, in some embodiments, in addition to wind direction mutations, other exogenous meteorological events, such as the passage of rapidly moving cloud edges, which cause step changes in shortwave radiation, can also be identified and used as natural experiments to test the dynamic response capability of the model. In summary, by introducing the identification and quantification of the advection effect, and using a series of methods such as instrumental variable regression and natural experiment verification, the physical authenticity and accuracy of the simulation of urban evapotranspiration process are improved.

[0057] The present embodiment identifies and solves the physical process of advection enhancement effect, which is generally ignored in existing urban hydrological models but is crucial for energy and water balance calculation, and improves the accuracy of urban evapotranspiration estimation. In the heterogeneous environment of cities, which is composed of hot (such as parking lots) and cold (such as green spaces) surfaces, local thermal advection is an important additional energy source. By setting the physical trigger conditions of the combination of clear sky, wind speed and temperature gradient along the wind direction, the spatio-temporal window in which advection effect is likely to occur is accurately identified. In order to solve the endogenous bias problem in estimating the advection enhancement term due to the mutual influence of temperature and evaporation, a weighted temperature difference along the wind path is constructed and used as an instrumental variable (IV). The instrumental variable is strongly correlated with the local temperature gradient, but is not directly affected by the local evaporation process error, so that an unbiased estimate of the advection enhancement coefficient α can be obtained through IV regression. In addition, by using events such as wind direction mutations as natural experiments for cross-validation, the reliability of the model is further guaranteed. By physically reasonable and statistically robust quantification and correction of this neglected energy term, a key link of urban surface energy balance is completed, making the estimation of total evapotranspiration more accurate, thereby providing a more solid data foundation for upper-layer water storage capacity assessment and optimization decision-making, especially when fine-grained assessment of local small-scale plots is required, the precision improvement brought by this correction is particularly crucial.

[0058] According to an aspect of the present application, before constructing the robust optimization model, it further includes: selecting a representative vegetation configuration value for the candidate plot, and calling the physical mechanism inference process to generate a mechanism anchor point containing the underlying surface hydrological response characteristics corresponding to each value; based on the mechanism anchor point, a semi-convex surrogate function at the plot level is constructed, wherein the semi-convex surrogate function is composed of the local linear approximation at each mechanism anchor point and the upper envelope of the quadratic shape penalty term; the semi-convex surrogate function is used to represent the relationship between vegetation configuration and precipitation assimilation capacity in the robust optimization model.

[0059] In this embodiment, the mechanistic anchor point refers to a set of discrete input-output relationship pairs obtained by running a complete and high-precision physical mechanism inference model. Specifically, for a candidate plot x to be optimized, several representative vegetation configuration values ​​are selected. Vegetation configuration can be quantified using leaf area index (LAI). For example, if the current LAI of plot x is 1.5, the set of representative values ​​that can be selected can be {1.5, 2.0, 2.5, 3.0}. Each representative LAI value is assigned a value of LAI. k As input, the event storage capacity S corresponding to all historical rainfall events under this configuration is calculated. x (LAI k The distribution of ) is determined, and the marginal response at that point, i.e., the gradient g, is calculated. k = dS x / dLAI | LAIk Thus, each pair (LAI) k S x (LAI k ), g k This constitutes a mechanistic anchor point. The semi-convex substitution function S hat_x (LAI) is used to replace the real physical model S in the optimization model. x (LAI) is an approximation function. Its core advantage lies in its ability to maintain high fidelity to the real physical model while possessing good mathematical properties (convexity or near-convexity), thus ensuring that the optimization problem can be solved efficiently. It is constructed by taking the upper envelope of the family of functions formed by all anchor point information, mathematically expressed as: S hat_x (LAI) = max k { S x (LAI k ) + g k ·(LAI -LAI k ) - (κ / 2)·(LAI - LAI k ) 2}; where S hat_x (LAI) represents the alternative event water storage capacity of plot x at any leaf area index (LAI); max k {} represents the operation of taking the maximum value of all anchor index k; S x (LAI k ) represents plot x at anchor point LAI k Event storage capacity at the location (mechanism anchor value); g k The marginal response at the anchor point (mechanistic anchor point gradient); κ is the non-negative upper bound coefficient of the second-order curvature; (κ / 2)·(LAI - LAI) k ) 2i.e. the quadratic shape penalty term. This penalty term makes the surrogate function not too curved, thus keeps good convexity.

[0060] Further, additional physical constraints and nonlinear processing are imposed on the construction of the semi-convex surrogate function. The construction of the semi-convex surrogate function at the plot level includes: imposing a slope upper bound constraint on the semi-convex surrogate function to ensure that the increase in vegetation configuration does not bring more than a physically reasonable upper limit of the increase in precipitation infiltration gain; and imposing a second-order curvature upper bound constraint on the semi-convex surrogate function to ensure that the response surface of the semi-convex surrogate function is smooth and complies with the physical law of diminishing marginal returns.

[0061] Specifically, the slope upper bound constraint refers to, when the mechanism anchor point is generated, limiting the calculated gradient g k , i.e. 0≤g k ≤g max . Wherein g max is a physical upper limit set based on plant physiology and soil physics knowledge, for example, the increase in water storage capacity per unit leaf area index cannot exceed a certain physical extreme value, avoiding unrealistic steep growth of the surrogate function in some regions. The second-order curvature upper bound constraint is achieved by selecting and setting the parameter κ. The setting of the value of κ makes the second-order derivative of S hat_x (LAI) everywhere not greater than κ, which mathematically ensures that the function does not oscillate too much, and physically naturally reflects the law of diminishing marginal returns, i.e. when the vegetation is very dense, the water storage benefit brought by the increase of unit vegetation will gradually decrease.

[0062] Optionally, the construction of the semi-convex surrogate function at the plot level further includes: when the mechanism anchor point indicates that the phase switching time jumps when the vegetation configuration takes different values, identifying the configuration interval where the jump occurs; independently constructing semi-convex surrogate functions on both sides of the configuration interval; and splicing the semi-convex surrogate functions on both sides at the boundary of the configuration interval, so as to retain the nonlinear jump caused by the physical mechanism in the final semi-convex surrogate function.

[0063] In the present embodiment, the jump of the phase switching time refers to the fact that the physical model shows that when LAI is slightly increased from one value to another value, the key feature (the phase switching time t star ) of the evapotranspiration process has a discontinuous, jump-like change. For example, when LAI changes from 2.0 to 2.1, t star may suddenly change from 5 pm to 11 am. This jump represents a qualitative change in system behavior, and if a smooth function is used to approximate it, important physical information will be lost. Illustratively, when the anchor point is generated, the LAI interval where the jump occurs (for example, [2.0, 2.1]) is recorded; semi-convex surrogate functions Shat_x_part1 (LAI), another S hat_x_part2 (LAI) is constructed for anchor points in the interval [2.1, 3.0] independently. hat_x (LAI) = S hat_x_part1 (LAI), when LAI ≤ 2.0; S hat_x (LAI) = S hat_x_part2 (LAI), when LAI > 2.0 (and smoothed at the boundaries). By piecewise construction and stitching, the surrogate function successfully preserves the nonlinear transition behavior driven by physical mechanisms, improving the fidelity of the model.

[0064] According to one aspect of the present application, as Figure 4 shown, determining the optimal vegetation configuration scheme comprises: based on the uncertainty information contained in the underlying hydrological response characteristics, constructing a distribution uncertainty set containing multiple possible probability distributions using a sensitivity-weighted distance metric; establishing a distribution robust chance constraint, requiring that under the worst probability distribution in the distribution uncertainty set, the probability that the total water storage increment meets the planning target is not less than the pre-set confidence level; converting the distribution robust chance constraint into a solvable deterministic constraint form, and combining with the engineering constraints to solve, to determine the optimal vegetation configuration scheme.

[0065] Specifically, the distribution uncertainty set P is a mathematical characterization of the uncertainty of the true probability distribution Q of future meteorological conditions and other random factors. This embodiment does not assume that Q follows any specific distribution (such as a normal distribution), but constructs a distribution sphere with a radius of Δ centered on the historical data empirical distribution P0, P = {Q | D(Q, P0) ≤ Δ}; where D(·,·) is a distance metric function used to measure the distance between two probability distributions. A preferred implementation is to use a sensitivity-weighted Wasserstein distance, which gives greater weight or penalty to the deviation of random variables that are more sensitive to the total water storage capacity (for example, net radiation in the initial period after rain), so that the uncertainty set can more intelligently cover the risk area that has the greatest impact on system performance. The mathematical form of the distribution robust chance constraint is: inf Q∈P Q(Σ x∈Ω [S hat_x (LAI x , ξ) - S x_base (ξ)] ≥ D target ) ≥ 1-ε; wherein inf Q∈P indicates finding the worst (infimum) case among all possible distributions Q in the uncertainty set P; Q(...) represents the probability of the event inside the parentheses occurring under distribution Q; Σ represents the sum of the set Ω of all candidate land parcels x; S hat_x(LAI x ξ) represents the plot x calculated using a substitution function under a random weather scenario ξ, in the configuration of LAI. x Water storage capacity below; S x_base (ξ) represents the baseline (current) water storage capacity; D target The total planned water storage increment target is defined by 1-ε, where 1-ε represents the preset confidence level, such as 95%. This constraint means that the proposed vegetation configuration scheme {LAI} is required to achieve the desired increase. x It is essential to ensure that the probability of the total increase in water storage meeting the target is at least 95%, even under the most unfavorable probability distribution. This is a very robust decision-making criterion.

[0066] Optionally, to solve this model, the probabilistic constraints described above need to be transformed into deterministic, computable constraints. The steps to transform the sub-Brook bar chance constraints into a solvable deterministic constraint form include generating a representative scenario set. This process further includes: using a time-block resampling method that preserves the autocorrelation structure of the dry season to generate the scenario set from historical data; thus, the scenario set can reflect the true meteorological time-series characteristics, allowing for a reliable approximation of the sub-Brook bar chance constraints. Specifically, instead of randomly selecting single-day records from historical meteorological data, continuous time segments (blocks) are extracted, such as 7 or 14 consecutive days. Because the meteorological conditions (such as radiation and humidity) during the dry season after rain have strong temporal persistence, this block sampling method can well preserve this autocorrelation structure, making the generated scenario set statistically closer to the real world, thereby making the constraint transformation results based on this scenario set more reliable. After generating the scenario set, the sub-Brook bar chance constraints can be precisely or approximately transformed into a series of linear constraints or second-order cone constraints, thus transforming the entire optimization problem into a convex optimization problem that can be efficiently solved by standard solvers (such as Gurobi, MOSEK). By combining this robust constraint with engineering constraints such as total investment cost not exceeding the budget and limited redevelopment area of ​​a single plot, and solving the problem with the objective function of minimizing total cost or maximizing benefits, the optimal vegetation configuration scheme {LAI} can be obtained. x} and the corresponding implementation list.

[0067] This embodiment addresses the technical problem of insufficient reliability and uncontrollable risk in optimization decision-making schemes due to unclear model mechanisms and lack of systematic consideration of uncertainties. Specifically, it solves this problem by constructing a semi-convex substitution function S. hat_x(LAI), which solves the problem of embedding complex high-precision physical models into large-scale optimization calculations. The surrogate function approximates the nonlinear and even transition-containing complex relationship between vegetation configuration and precipitation assimilation capacity with high fidelity in the form of a computationally efficient convex function through mechanism anchor points and quadratic shape penalties. On this basis, a distributed robust chance-constrained (DROCC) optimization model is established. This model does not rely on any specific probability distribution assumption for future meteorological scenarios, but instead constructs an uncertainty set containing all possible distributions close enough to historical data, and under the worst (most unfavorable) probability distribution within this set, the probability of meeting the total water storage increment target is not less than the pre-set confidence level (such as 95%), which can effectively resist deep uncertainty. The output optimal vegetation configuration scheme is no longer a simple expected optimal solution, but a robust optimal solution with clear reliability guarantee. The risk-quantified decision-making basis provided for urban planners enables the massive green space system reconstruction project to achieve its predetermined flood control and disaster reduction and rainwater resource utilization goals with a high probability under various possible future meteorological conditions, improving the scientificity and forward-looking nature of urban infrastructure planning.

[0068] In another embodiment of the present application, a continuous phase change model is used to depict the transition of urban land surface from wet to dry, so that the physical reality can be more realistically simulated under certain conditions (e.g., slow recovery of sunshine after rain, diverse underlying surface materials leading to different drying times). Specifically, the quantitatively obtained underlying surface hydrological response characteristics can also be used for:

[0069] The urban heat source separation and multi-scale smoothing processing of the surface temperature data contained in the rainfall event data set is performed to obtain the clean surface temperature time series; a continuous phase change weight function is constructed to represent the gradual transition process from the wet surface dominant stage to the soil restricted stage; the specific form of the phase change weight function is determined by joint optimization of the objective function, wherein the joint optimization objective function couples the phase change rate driven by the dynamic characteristics of the clean surface temperature time series, the energy consistency and the surface resistance continuity; based on the determined phase change weight function, the continuous decomposition and calculation of the evapotranspiration process are performed to generate the underlying surface hydrological response characteristics.

[0070] In this embodiment, the phase transition weight function φ(t) is a scalar function of time t, which takes values in the range [0, 1]. φ(t) = 1 means that the evapotranspiration process is completely dominated by wet surface (intercepted water, impermeable surface water film) at time t; φ(t) = 0 means that it is completely limited by soil moisture supply; and 0 < φ(t) < 1 means that the system is in a transition state of mixed action of the two mechanisms. An exemplary function form can be a Sigmoid function: φ(t) = 1 / (1 + exp((t - tc) / σ)); where tc is the center time of the phase transition process, i.e., the time when φ(tc) = 0.5; σ is the width parameter of the transition zone, representing the speed of the phase transition process. The total evapotranspiration model is correspondingly expressed as a continuous weighted sum: ET model (t) = φ(t)·ET wet (t) + (1-φ(t))·ET soil (t);where ET model (t) is the total evapotranspiration of the model; ET wet (t) is the evapotranspiration dominated by pure wet surface (for example, a fast exponential decay process); ET soil (t) is the evapotranspiration limited by pure soil (for example, a slower exponential decay process). Through a highly physically constrained joint optimization process, the φ(t) function form that best describes the rainfall event is determined (i.e., its parameters tc and σ are determined).

[0071] The key premise for accurate implementation of this joint optimization process is to obtain high-quality ground surface temperature signals that reflect the real evapotranspiration physical process. Therefore, further, to obtain clean ground surface temperature time series, a human heat flux model is established, which includes traffic intensity, building density and daily cycle terms, to represent the influence of human heat sources on ground surface temperature; the to-be-estimated coefficients of the human heat flux model are determined by a variational fitting target that takes into account the minimization of energy balance residual and spatial smoothing constraints; the influence of human heat flux estimated by the human heat flux model with the to-be-estimated coefficients is separated from the ground surface temperature data to obtain the clean ground surface temperature time series.

[0072] In this embodiment, the heat generated by human heat sources (such as traffic, industry, building air conditioning heat dissipation) will be superimposed on the natural solar radiation effect, interfering with the judgment of the latent heat dissipation process of evapotranspiration. Therefore, a semi-empirical human heat flux model H anthro (t, x) = β1·Traffic(t, x) + β2·Building density (x) + β3·sin(2π·t / 24+ φ phase ); where H anthro(t, x) represents the anthropogenic heat flux at location x at time t; Traffic(t, x) represents the traffic intensity index at location x at time t (which can be obtained from traffic models or remote sensing data); Building density (x) represents the building density at location x; the sin(...) term is used to fit the diurnal variation of anthropogenic thermal activity; β1, β2, and β3 are the coefficients to be estimated; φ phase The phase angle is the daily cycle angle. To determine these coefficients to be estimated, this embodiment establishes and solves the variational fitting objective function: J anthro =∫∫|Rn(t,x) - H(t,x) - LE(t,x) - G(t,x) - H anthro (t, x)| 2 dxdt + ω*∫∫|▽H anthro (t, x)| 2 dxdt; where J anthro Let H be the objective function; the first integral term minimizes the residual of energy balance over the entire spatiotemporal domain, where Rn, H, LE, and G are the net radiation, sensible heat, latent heat, and surface heat flux, respectively; the second term is a spatial smoothing regularization term, where ▽ is the spatial gradient operator, and ω* is the smoothing weight. This term makes the estimated H... anthro The space is smooth and continuous, avoiding overfitting. H is obtained by solving this minimization problem. anthro The specific values ​​are obtained, and then the influence is removed from the original surface sensible heat flux or surface temperature to obtain the clean surface temperature time series, which solves the problem of anthropogenic thermal pollution of urban surface temperature (LST) signals.

[0073] In an optional embodiment, based on obtaining the time series of clean surface temperatures, in order to further enhance the joint optimization objective function J phase The physical constraints introduce the following steps: the continuity of coupled surface resistance in the joint optimization objective function is achieved by a penalty term aimed at minimizing the difference between the intercepted surface resistance weighted by the phase transition weight function and the impermeable surface resistance weighted inversely by the same function.

[0074] In this embodiment, based on physical insights, the resistance of the land surface to evapotranspiration is closely related to the surface water content. Moist vegetation or a water film (corresponding to φ(t) ≈ 1) has a surface resistance r... s_wet The surface resistance r of dry soil (corresponding to φ(t) ≈ 0) is very small. s_soil It is very large. Therefore, a physically reasonable phase transition process φ(t) should be related to the phase transition from r... s_wet to r s_soil Coupled with the smooth transition process. Specifically, by jointly optimizing the objective function J phase Add a penalty item Φ resistTo achieve this coupling: Φ resist = ∫|φ(t)·r s_wet - (1-φ(t))·r s_soil - r s_eff (t) 2 dt (a more complete form); where r s_eff (t) is the instantaneous effective surface resistance estimated by other means (e.g. Penman-Monteith inverse algorithm). Alternatively, in a simplified implementation, the penalty term aims to minimize the difference between the weighted resistance model and its expected behavior. With this constraint, the solution of phase change weight no longer depends only on the evapotranspiration itself, but also on the evolution of the surface energy exchange efficiency (characterized by the surface resistance), enhancing the physical realism of the model. Assembling the above elements into a complete joint optimization objective function J phase , to determine the optimal parameters of the phase change weight function φ(t). The schematic structure of the objective function J phase (tc, σ) is: J phase = (data fitting error term) + λ1·(phase change rate constraint term) + λ2·(surface resistance continuity penalty term) + λ3·(smoothing regularization term). The phase change rate constraint term utilizes the dynamic characteristics of clean land surface temperature time series. Specifically, in the process of the surface changing from wet to dry, its temperature rising rate will experience a process of acceleration and then deceleration, and this inflection point information can be captured by the first and second time derivatives of the land surface temperature. An indicator Ψ(t) = (d 2 LST clean / dt 2 ) / (dLST clean / dt +ε) can be constructed, and the peak position of the indicator Ψ(t) should theoretically correspond to the peak position of the phase change rate dφ / dt. Therefore, the phase change rate constraint term ∫|dφ / dt - f(Ψ(t))| 2 dt is to require that the change rate of φ(t) should be consistent with the physical dynamics revealed by the clean land surface temperature. By minimizing J phase , the optimal tc and σ can be obtained, and thus the function φ(t) is determined. Based on this function and the weighted total evapotranspiration model ET model (t), the subsequent integral calculation can be carried out to obtain the underlying surface hydrological response characteristics such as event water storage capacity.

[0075] Optionally, in areas where data quality is extremely high or urban morphology is relatively simple, the anthropogenic heat separation model can be simplified, for example, only considering the daily cycle term. In some embodiments, the phase transition weight function φ(t) can also take other functional forms, such as an error function or a piecewise polynomial, as long as it can meet the boundary condition of smooth transition from 1 to 0.

[0076] The embodiment provides a more refined and physically more realistic evapotranspiration process simulation path, solves the problems of anthropogenic heat pollution of ground surface temperature signals in complex urban environments and gradualness of evapotranspiration stage transition. Specifically, the anthropogenic heat flux model H anthro is obtained by separating the anthropogenic heat influence from the original land surface temperature (LST) data, and using a variational fitting algorithm that takes into account energy closure and spatial smoothing, to obtain clean LST time series that only reflect natural heat processes, and to provide unbiased data input for subsequent analysis that relies on temperature dynamics. On this basis, instead of finding a discrete switching point, a continuous phase transition weight function φ(t) is defined to describe the limited gradual transition from a wet surface to soil. The determination process of φ(t) is a highly coupled joint optimization, which integrates three major physical constraints of evapotranspiration data fitting, phase transition rate driven by the kinetic characteristics (first and second derivatives) of clean LST time series, and surface resistance continuity that should be followed by water and heat state transition of the ground surface into a unified objective function J phase . The final phase transition process description not only fits the water quantity data, but also keeps consistent with the internal rhythm of energy change and resistance change. It can capture real physical processes that are not instantaneous and gradual changes with high fidelity, and due to the built-in urban heat source separation module, it can still maintain high accuracy and reliability in the core area of the city where human activities are intense.

[0077] Overall, the application discloses a city scale precipitation assimilation and vegetation configuration collaborative optimization evaluation method, comprising: integrating city multi-source heterogeneous data and identifying rainfall events; identifying and quantifying the phase switching process from wet surface dominated to soil limited evapotranspiration by building multiple physical gates such as energy consistency, wet surface dissipation and geometric radiation sensitivity, or using a continuous phase change model coupled with city heat source separation and surface resistance continuity, so as to accurately obtain the underlying surface hydrological response characteristics; on this basis, a distribution robust chance constraint optimization model coupled with mechanism anchor and semi-convex surrogate function is built, and the optimal vegetation configuration scheme is solved under the premise of considering deep uncertainty; and the optimal scheme is back substituted for closed loop checking and effect evaluation. The application can accurately quantify the hydrological response characteristics of the city underlying surface, and under uncertain conditions, the vegetation configuration scheme with clear probability and controllable risk is formulated, and the scientificity and reliability of the city water resource resilience planning are improved.

[0078] Exemplarily, to back substitute the optimal vegetation configuration scheme to the rainfall event data set, perform closed loop checking and effect evaluation, and generate an optimization evaluation report, a possible implementation manner is provided. Through the feedback verification cycle, the theoretically optimal scheme can also achieve the expected effect in the simulated real world, so as to provide a solid decision basis for the final engineering implementation. Specifically, the closed loop checking and effect evaluation process can include the following steps:

[0079] An optimal vegetation configuration scheme is obtained as the final output, which is usually represented as a list detailing the target leaf area index LAI x_optimal or the vegetation type to be selected for each to-be-transformed plot x. These optimization results are written back to the benchmark rainfall event data set. Specifically, the writing back operation not only includes updating the leaf area index parameter of each plot x, but also can include the linkage updating of other physical parameters related thereto. For example, if the scheme requires transforming a grassland (lower LAI) into a shrubbery (higher LAI), in addition to updating the LAI value, the ground albedo, zero plane displacement height, ground roughness and a series of parameters describing the physical characteristics of the underlying surface of the plot need to be adjusted accordingly according to the change of the vegetation type. Through this step, a new data set representing the post-implementation scenario is generated, which is physically consistent and is a bridge connecting optimization decision and physical verification.

[0080] On the basis of generating the post-implementation scenario data set, the complete physical mechanism inference chain is called for re-computation. Preferably, the physical mechanism inference method and the advection enhancement correction method can be called. For each plot and each historical rainfall event, the phase switching time (or continuous phase change process), the time series of each evapotranspiration component, the key hydrological parameters (such as decay time constant) and the final event water storage capacity S under the optimal vegetation configuration are re-calculatedevent_after This process is equivalent to a full virtual implementation of the optimization scheme in the computer, and thus the simulation of the entire urban hydrological response after the implementation of the scheme.

[0081] The target achievement rate and side effect assessment are carried out, mainly from two dimensions: the degree of achievement of the expected target and the potential unintended impact. On the one hand, the target achievement rate is quantitatively evaluated. The event water storage capacity S event_after after implementation is compared with the baseline water storage capacity S event_before before implementation, and the total water storage increment in the entire study area Ω and all historical events E is calculated. The target absorption rate η can be calculated as: η = [∑ x∈Ω，e∈E S event_after (x, e) - ∑ x∈Ω，e∈ E S event_before (x, e)] / D targettotal ; where D targettotal is the total amount of precipitation to be absorbed set by the plan. In addition, in order to evaluate the reliability of the results, the uncertainty propagation also needs to be carried out. The uncertainty (e.g. confidence interval) of each S event value obtained in the mechanism inference process is transferred to the final absorption rate η through Monte Carlo simulation or error propagation formula, giving a confidence range. For example, the evaluation result can be expressed as that the scheme is expected to complete 98% of the absorption target, with a 90% confidence interval of [91%, 105%]. This evaluation with uncertainty provides more complete risk information for decision makers. On the other hand, the potential side effects, especially the thermal environmental impact, are evaluated. The change of vegetation configuration will affect the distribution of surface energy. For example, increasing vegetation coverage can increase evapotranspiration (latent heat, LE), but may also change the surface sensible heat flux (H). Sensible heat flux is one of the main driving forces of urban heat island effect. Therefore, it is necessary to compare the change of average sensible heat flux of the reconstructed land in the afternoon (such as 13:00-15:00) under typical summer sunny and hot weather conditions before and after implementation. If the sensible heat flux increases significantly, it may indicate that the scheme may exacerbate the deterioration of the local thermal environment while enhancing the absorption of precipitation, which needs to be weighed by the decision makers. Alternatively, other aspects of side effects, such as the impact on water quality (different vegetation has different interception capacity for pollutants) or ecology (species singularity), can also be evaluated.

[0082] All the analysis results, including the optimal vegetation configuration scheme, the spatial distribution of expected water storage increment, the target achievement rate and its confidence range, the investment and maintenance cost estimation, and the thermal environment risk assessment indicators, are integrated and written into a structured technical report. This report can be used for project delivery or public communication. At the same time, the key thresholds used in the entire evaluation process (such as event identification thresholds), shape constraint parameters (such as kappa value), trigger conditions, model parameters, solver logs, etc. are archived into a method metadata package. This metadata package ensures the traceability and reusability of the entire research process, and when future evaluations of other cities or in the next year are needed, model migration and updating can be quickly carried out based on this metadata package.

[0083] In a specific embodiment, assume the following scenario: Candidate plot X: located in the city center, area 10000 square meters, current status is grassland, its average leaf area index (LAI) is 1.0. Rainfall event Y: a 3-hour-long rainfall with a total rainfall of 20 mm stops at t = 0. After the rain stops, the total evapotranspiration water depth time series ET obs (t) (unit: mm / hour) of the plot is obtained through sensor network and remote sensing data. The ET obs (t) sequence after the rain stops is fitted using a two-section exponential decay model. It is assumed that a traversal search is performed within the window of t = 1 to t = 24 hours, and the residual sum of squares SSE(tc) is calculated for each potential partition point tc. The calculation results show that SSE(12) reaches a minimum value when tc = 12 hours. Therefore, the candidate switching time is preliminarily determined as tc = 12 hours. Physical gating joint verification is performed: energy consistency gating: extract the energy flux data from t = 13 to t = 24 hours, and perform linear regression on (LE(t) + H(t)) and (Rn(t) - G(t)). The slope of the regression equation is 0.97, and the intercept is 6.5 W / m2. It is assumed that the preset passing criteria are that the slope is between [0.9, 1.1] and the intercept is between [-10, 10] W / m2. This candidate time passes this check. Wet surface dissipation gating: calculate the wet surface dissipation index W(t) from t = 13 to t = 15 hours, and the obtained values are {0.82, 0.88, 0.91}. It is assumed that the preset threshold is 0.8, and since the index exceeds the threshold for 3 consecutive hours, the candidate time passes this check. Geometric radiation sensitivity gating: (this is a simplified example, assume that the check passes). Since the candidate time tc = 12 hours passes all three physical gating, it is confirmed as the final phase switching time, i.e. t star = 12 hours. For the ET obs (t) data of t > 12 hours, parameter estimation is performed using a joint objective function. By minimizing the objective function, the soil-restricted phase decay time constant τ s= 72 hours (i.e., 3 days), the baseline dissipation strength ET0 = 0.0417 mm / hour (i.e., 1.0 mm / day). The effective water storage capacity S event of the plot for this rainfall event, under the current status (LAI = 1.0), can be approximated as S eventbase ≈ τ s · ET0 = 3 days · 1.0 mm / day = 3.0 mm. In volume, this is 3.0 mm · 10,000 m2= 30 m3.

[0084] Now consider upgrading the vegetation of the plot, e.g., raising the LAI to 2.0. Repeating the entire process of physical mechanism deduction above, with LAI = 2.0 as input, we obtain the new event water storage capacity S event_upgraded = 4.5 mm. Thus, we have two mechanism anchors: Anchor 1 (k = 1): LAI1= 1.0, S x (LAI1) = 3.0 mm; and Anchor 2 (k = 2): LAI2= 2.0, S x (LAI2) = 4.5 mm. Between Anchor 1 and Anchor 2, we can estimate the average gradient g 1_2 = (4.5 - 3.0) / (2.0 - 1.0) = 1.5 mm / LAI unit. Based on these two anchors, and assuming a second-order curvature upper bound coefficient K = 0.5, we can construct a semi-convex surrogate function S hat_x (LAI) for plot X. This function is an upper envelope (i.e., taking the maximum) of two functions f1(LAI) and f2(LAI): f1(LAI) = S x (LAI1) + g1 · (LAI - LAI1) - (K / 2) · (LAI - LAI1) 2 ; and f2(LAI) = S x (LAI2) + g2 · (LAI - LAI2) - (K / 2) · (LAI - LAI2) 2 ; substituting the numerical values (and assuming g1and g2are both approximated as g 1_2 ), we obtain: S hat_x (LAI) = max{ 3.0 + 1.5 · (LAI - 1.0) - 0.25 · (LAI - 1.0) 2 , 4.5 + 1.5 · (LAI - 2.0) - 0.25 · (LAI - 2.0) 2}. This specific function expression S hat_x(LAI), i.e. the input of plot X to the whole city-scale robust optimization model. It describes the non-linear relationship between the precipitation absorption capacity of the plot and the vegetation configuration (LAI) in a computationally efficient and physically reasonable way, making large-scale, multi-plot collaborative optimization possible.

[0085] In an embodiment of the present application, a more detailed and preferred implementation is provided for data integration and quality control. In the complex urban environment, the quality of multi-source heterogeneous data is the key to determine the success of all subsequent analysis. Therefore, this embodiment describes a set of refined data preprocessing and quality assurance processes to solve the problems of spatial and temporal mismatch, physical inconsistency, noise interference, etc., thereby providing a high-quality and reliable data base. Specifically, a preferred data preprocessing and quality control process includes:

[0086] High-fidelity spatial and temporal resampling. In order to fuse data of different sources and different resolutions into a unified grid and time step, a method with higher fidelity is adopted. In spatial resampling, unlike simple nearest neighbor or bilinear interpolation, the method of geometric overlap area weighted projection is preferred. This method ensures the conservation of area type variables such as land cover type proportion when converting between different resolutions. In temporal resampling, different physical quantities are treated differently. For cumulative quantities such as rainfall, a simple summation is used; for average quantities such as temperature, a time-weighted average is used. This differentiated treatment ensures that the physical meaning of the variables is not distorted during the time scale aggregation process.

[0087] Multi-stage, physical and statistical quality control. After the data is spatially and temporally aligned, the following quality control steps are performed in sequence: Physical feasibility check: used to eliminate data points that obviously violate physical laws. A core check is the surface energy balance constraint. That is, at any time t, the sum of latent heat flux LE(t) and sensible heat flux H(t) should not exceed the available energy, i.e. the difference between net radiation Rn(t) and ground heat flux G(t). Considering sensor errors, the constraint can be written as (LE(t) + H(t)) ≤ (Rn(t) - G(t)) + ε tol ; where ε tol is a small energy residual tolerance upper limit, for example 10 watts per square meter. Data points that do not meet this constraint will be marked as suspicious. Robust outlier suppression: sensor data in urban environments is often subject to local, transient strong interference, forming outliers. To effectively identify these outliers, a robust statistical method based on a sliding window is preferred. Specifically, the dimensionless outlier score S out (x, t) = [X(x, t) - median window (X)] / MAD window(X); where X(x, t) is the current observation of the variable; median window (X) is the median of the variable within the same sliding time window; MAD window (X) is the median of the absolute deviation within the window. This score is not sensitive to the overall shape of the distribution of the data, and is more robust than the traditional Z-score method based on mean and standard deviation. When S out (x, t) exceeds a pre-defined threshold (e.g. 3.5), the data point is identified as an outlier and repaired. Cross-variable correlation harmonization: There can be systematic biases between different sensors. To harmonize the consistency between different components of the energy balance, a dynamic harmonization method is preferred. Specifically, within a sliding time window, a linear regression is performed on (LE(t) + H(t)) vs. (Rn(t) - G(t)). In an ideal, unbiased system, the regression slope should be 1 and the intercept should be 0. If the regression result deviates significantly, it indicates that some fluxes are systematically overestimated or underestimated. At this time, the values of LE(t) and H(t) can be adjusted in proportion according to the regression coefficients, so that the relationship between LE(t) and H(t) and the net available energy Rn(t) - G(t) is more stable and consistent over the entire time series.

[0088] Through the above preferred preprocessing steps, a high-quality data set with self-consistent internal physical logic, low noise level, and unified spatiotemporal reference can be obtained, providing a solid guarantee for the accuracy of subsequent physical mechanism inference.

[0089] In another embodiment of the present application, a more robust and refined implementation is provided for the physical mechanism inference process, particularly the physical gating check and parameter estimation involved therein. By introducing stronger physical constraints and better statistical methods, the reliability of stage switching time identification and the accuracy of subsequent hydrological parameter estimation are improved. Specifically, the preferred improvements to the physical mechanism inference process include:

[0090] On the basis of the three physical gates, the judgment criteria are deepened. The enhanced test of energy consistency gate includes: in addition to testing the slope and intercept of linear regression, the second-order statistic test is added, that is, it is required that the variance of relative energy residual within the time window after the candidate switching time should be significantly reduced compared with that before the switching. The relative energy residual is defined as (LE+H)-(Rn-G) divided by (Rn-G). The decrease of variance shows that the model after switching not only conforms to the energy balance in mean, but also has smaller fluctuations around the equilibrium state, that is, the energy distribution mechanism is more stable, which provides stronger evidence for the occurrence of stage switching. The adaptive weight of the wet surface dissipation index is increased: the weights a, b, c of the three components of the wet surface dissipation index W(t) are preferably not fixed values, but are determined by dimensionless sensitivity scaling. Specifically, during the model training phase, the sensitivities ΞET / ΞLST, ΞET / ΞAlbedo, ΞET / ΞRs of the total evapotranspiration ET to the land surface temperature LST, the land surface albedo Albedo and the shortwave radiation Rs can be calculated respectively by numerical perturbation experiment on historical data, where Ξ is the partial derivative. After dimensionless processing and normalization of these sensitivities, the weight coefficients a, b, c with clear physical meaning can be obtained. This adaptive scaling makes the construction of the index W(t) reflect the dominant degree of different physical processes on evapotranspiration in a specific research area.

[0091] In order to determine the decay time constant τ of the soil limited stage s A joint objective function is constructed with the reference dissipation intensity ET0. An optimal implementation is provided for the construction and solution of the function. The complete form of the joint objective function can be expressed as: J total (t star , τ s , ET0) = Σ t ρ(ET obs (t) - ET model (t)) + β·Σ t |ε rel (t)| + μ·Penalty Gates (t star );where the first term is the data fitting term, and the loss function ρ(·) preferably adopts a robust loss function that is not sensitive to outliers, such as the Huber loss function. Compared with the conventional square loss, the Huber loss behaves as a square loss when the residual is small, and behaves as a linear loss when the residual is large, thereby reducing the excessive influence of a single abnormal data point on the fitting result. The second term is the energy consistency penalty term, and ε rel(t) is the relative energy residual, and β is its penalty weight. This term enforces energy conservation as a soft constraint, guiding the parameters to converge towards physically more reasonable values during optimization. The third term is the physical gating penalty term, Penalty Gates (t star ) is the weighted sum of three physical gating violations. This term enforces that the phase switching time t star itself can also be fine-tuned within its candidate neighborhood during the final joint minimization process, so as to find a globally optimal solution in terms of statistical fitting and physical reasonability. By simultaneously optimizing the phase switching time t star and the key physical parameters τ s and ET0 of the subsequent phase in a unified objective function, and coupling robust statistical methods with multiple physical constraints, more accurate, reliable and physically consistent characteristics of the underlying surface hydrological response can be obtained than traditional step-by-step solution methods.

[0092] In another embodiment of the present application, a preferred implementation is provided for the mathematical construction and solution strategy of the distributionally robust optimization model. The distributional uncertainty set P is used to describe the uncertainty of the real probability distribution of future weather scenarios. The preferred construction method provided in this embodiment not only considers historical data, but also incorporates sensitivity information of the system to different uncertainty factors. Specifically, the Wasserstein distance (also known as the optimal transport distance) is preferably used to construct the uncertainty set. Compared with other statistical distances, the Wasserstein distance can measure the geometric structural difference of the distribution and performs better when dealing with high-dimensional random variables. The uncertainty set P is defined as a distribution sphere with a radius of Δ centered on the empirical distribution P0 of the historical data: P = {Q | D W (Q, P0) ≤ Δ}; where D W (Q, P0) is the Wasserstein distance between Q and P0. D W here is the sensitivity-weighted Wasserstein distance. The specific implementation is to introduce a sensitivity weight matrix W sens when defining the movement cost function c(ξ1, ξ2) relied on by the Wasserstein distance. The cost function can be defined as: c(ξ1, ξ2) = ||W sens ·(ξ1 -ξ2) ||; where ξ1, ξ2 are two different weather scenario vectors; W sens is a diagonal matrix, and the elements w iThe sensitivity of the total water storage increment to the i-th component of the meteorological vector (e.g., the net radiation at the 24th hour after rain). This sensitivity can be obtained by running the physical model and performing a perturbation analysis. By introducing this sensitivity weight, the constructed uncertainty set P is asymmetric in shape. It extends further in the directions of uncertainties that have the most impact on the system performance (i.e., the directions with high sensitivity), i.e., more extreme changes are considered. In other words, this uncertainty set can intelligently and purposefully cover those most critical risks, making the robust optimization decision based on it more efficient and targeted.

[0093] To enable numerical solution, the general form of distributionally robust chance constraint (DROCC) needs to be transformed into a deterministic mathematical programming problem. This embodiment provides a preferred transformation path. A time block resampling method that preserves the autocorrelation structure of dry periods is employed to generate N high-quality representative scenarios {ξ1, ξ2,..., ξN} from historical data. N It ensures that the discrete scenario set used for approximation can reflect the temporal persistence characteristics of the real meteorological process. Based on the powerful optimization duality theory, it can be proved that the distributionally robust chance constraint under the uncertainty set constructed based on the sensitivity-weighted Wasserstein distance can be accurately or closely approximated as a deterministic large-scale convex optimization problem. In a preferred implementation, this problem can be transformed into a second-order cone program (SOCP) problem. Although SOCP is nonlinear, it belongs to the category of convex optimization, and there are mature and efficient interior point method solvers (such as MOSEK, Gurobi) that can find its global optimal solution in polynomial time.

[0094] In another optional implementation, by employing a simpler distance metric (such as variational distance under a specific norm) or linearization approximation of the chance constraint (such as Conditional Value-at-Risk, CVaR approximation), this problem can also be transformed into a large-scale linear programming (LP) problem. Although the accuracy is slightly lost, the solving speed of LP is faster, especially suitable for extremely large-scale urban planning problems.

[0095] In summary, this embodiment provides a complete and computationally feasible path for solving the deep uncertainty problem in urban vegetation planning by employing sensitivity-weighted Wasserstein distance to refine the construction of the uncertainty set and combining the solving techniques based on block sampling and convex optimization reconstruction.

[0096] In another embodiment of the present application, the urban unique anthropogenic heat physical process, the thermodynamic process of surface energy exchange, and the hydrological process of surface evapotranspiration are deeply coupled, so as to realize a highly realistic simulation of the urban post-rain drying process. Specifically, the specific model and solving process of the physical separation and inversion of the urban anthropogenic heat source are unfolded. The specific form and parameters of the anthropogenic heat flux model H anthro (t, x) need to be adapted according to the data availability of the research area. For example, in a data-rich smart city environment, the traffic intensity index Traffic(t, x) can be provided by real-time traffic flow data; in areas with less data, a static road network density can be used as a proxy. The building density Building density (x) can be extracted from a high-precision building model. The solving of the variational fitting objective function J anthro is a large-scale optimization problem, and the core is to balance the two goals of energy closure and spatial smoothing. The selection of the weight coefficient ω* is crucial, which controls this balance. If the value of ω* is too small, the estimated H anthro will have a lot of noise, although it can well close the energy at each grid, but the spatial distribution is unreasonable; if the value of ω* is too large, an overly smoothed H anthro distribution will be generated, which will mask the real heat source differences in the city. Preferably, the value of ω* can be determined by the cross-validation method, that is, fitting the model on part of the data, and then verifying the prediction error on another part of the data, and selecting the ω* value that minimizes the prediction error.

[0097] The phase change weight function φ(t) is finally determined by solving the joint optimization objective function J phase that couples multiple physical processes. A preferred J phase (tc, σ) has the complete form: J phase =∑[ET obs (t) - ET model (t)] 2 +λ1·∫|dφ / dt(t) - f(Ψ(t))| 2 dt+λ2·Φ resist (t)+λ3·∫|d 2 φ / dt 2 (t)| 2 dt; where λ1, λ2, λ3 are the penalty weights of each term, which can also be determined by cross-validation and other methods. The four components of the objective function each carry different physical constraints: the data fitting term ∑[ET obs (t) - ET model (t)] 2 is the basic term, which requires the model to predict the total evapotranspiration ETmodel (t) = φ(t) · ET wet (t) + (1 - φ(t)) · ET soil (t) should be as close as possible to the observed value ET obs (t). The phase rate constraint term λ1 · ∫|dφ / dt(t) - f(Ψ(t))| 2 dt links the rate of hydrological state transition dφ / dt to the dynamic characteristic of the surface thermal state Ψ(t). Ψ(t) is a dynamic indicator calculated from the time series of clean surface temperature, reflecting the changes in the rate of surface heating. This constraint requires the rhythm of hydrological phase transition to be consistent with the rhythm deduced from thermodynamics, which is a deep data-driven physical coupling. The surface resistance continuity penalty term λ2 · Φ resist (t) requires the evolution of the effective surface resistance implied by the phase transition weight φ(t) to be consistent with the physical model of the transition from the wet surface resistance r s_wet to the dry surface resistance r s_soil . This constraint couples hydrology and micrometeorology theory. The smoothness regularization term λ3 · ∫|d 2 φ / dt 2 (t)| 2 dt is a mathematical regularization term to penalize excessive bending or oscillation of the phase transition weight function φ(t), ensuring that it evolves smoothly as a physical process. By minimizing this highly coupled joint objective function containing data fitting, thermodynamic dynamics, micrometeorological resistance model, and mathematical regularization, the embodiment can solve a physically plausible continuous phase transition process φ(t), thereby achieving a fine and dynamic simulation of the urban evapotranspiration process.

[0098] In one exemplary embodiment, a macro-scale estimation method of urban water storage capacity based on a single exponential decay model can include the following steps: obtain daily evapotranspiration (ET) time series data of a study area over a period of time (e.g., 2009-2018). The daily ET data can be simulated by a mechanistic model, such as a four-source evapotranspiration model, driven by atmospheric forcing data (e.g., CMFD dataset), to provide a basis of continuous daily-scale ET data for subsequent decay analysis. Identify and select all eligible dry spell events. A dry spell event is defined as a continuous period of no precipitation. To ensure the effectiveness of the decay characteristics, a series of selection criteria can be applied to the dry spell events. In one specific implementation, the selection criteria can include: dry spell length criterion: requires the dry spell to last at least three consecutive days without effective precipitation; irrigation influence avoidance criterion: to reduce the disturbance of urban green space irrigation on the natural evapotranspiration decay process, the upper limit of the dry spell length for analysis can be limited to within 10 days; non-snow condition criterion: requires the daily average temperature to be above 0°C throughout the dry spell to exclude the influence of snowfall or frozen soil conditions on evapotranspiration. For each selected dry spell event, fit a single exponential decay model to estimate the initial evapotranspiration ET0 and the decay time scale λ. The core assumption is that the consumption process of total dynamic water storage (including soil water, interception water, etc.) during the entire dry spell can be approximated by a single exponential decay process. The initial evapotranspiration ET0 is defined as the evapotranspiration value on the first day of the dry spell, and the decay time scale λ is the time constant of the exponential decay process. According to the estimated ET0 and λ, calculate the urban water storage capacity S. The urban water storage capacity S is defined as the total amount of water that can be removed from the system through evapotranspiration during a complete drying process. Its calculation formula is the simple product of the two parameters: S = λ·ET0; where S is the total dynamic water storage consumed by evapotranspiration during complete drying, with a unit of millimeters (mm); λ is the time scale, with a unit of days (day); ET0 is the evapotranspiration on the initial day of the dry spell, with a unit of millimeters per day (mm / day). Perform statistics and filtering on the calculation results of all dry spell events to obtain the multi-year average water storage capacity of the region. To obtain a stable and reliable regional water storage capacity representation, the S values calculated for all dry spell events need to be statistically averaged. Before statistics, reasonable filtering of the estimation results of single events can also be performed, for example, extreme or abnormal results with λ greater than 35 days or ET0 greater than 10 millimeters per day are removed. The multi-year average S value obtained is the estimated urban water storage capacity.

[0099] It can be understood that this embodiment has the advantages of simplicity and convenience of application as a macro-estimation model. However, compared with the above-mentioned embodiments of the present application, it has several inherent technical limitations, and these limitations are exactly the technical problems to be solved by the present application: regarding the complex urban water storage system (including interception, depression, soil, etc.) as a single, undifferentiated water tank, and the evaporation process is controlled by a uniform decay constant λ. This black box treatment cannot reveal the contribution and dynamic evolution of different internal physical processes. In contrast, by identifying phase switching or simulating continuous phase change, the present application can clearly distinguish and quantify the two stages of wet surface dominated and soil limited, which have completely different physical mechanisms, and provide a more profound physical insight. This embodiment is based on daily scale data for calculation, and the time scale λ obtained is also at the level of days. This makes it unable to capture the more intense hydrological dynamics that occur in the early stage after rain (usually within a few hours to a day), which are dominated by the rapid evaporation of interception water and surface water film. The present application can handle higher time resolution data (such as hourly), and accurately depict this rapid response process through component decomposition and continuous phase change weight, which is crucial for evaluating the rapid storage capacity of cities in response to short-term heavy rainfall. Using this embodiment, some abnormal phenomena may be observed, such as in some cities, increasing vegetation (increasing LAI) leads to a decrease in total water storage capacity S. The model of this embodiment itself cannot explain the internal reason of this phenomenon. In order to explain it, additional external attribution analysis is needed to decompose the ET component (such as vegetation transpiration Et, impermeable surface evaporation Eb). The present application can directly diagnose the reason of this phenomenon within the model, for example, the shading inhibition effect of increased vegetation on impermeable surface evaporation Eb exceeds the increase of its own transpiration Et, thereby providing direct, mechanism-level explanation and decision support. The application of this embodiment relies on the key assumption that the water storage is completely filled after each rainfall. In reality, for small rain events that occur continuously or have insufficient rainfall, this assumption often does not hold, which will introduce estimation bias. The present application does not rely on this strong assumption, but simulates the actual hydrological response process of each event by strictly adhering to a series of physical constraints such as energy balance and surface resistance continuity, starting from the first principle of physics, and thus has stronger universality and accuracy.

[0100] In summary, through comparison with this comparative embodiment, it can be seen that the urban-scale rainfall storage and vegetation configuration co-optimization evaluation method proposed by the present application overcomes the limitations of existing macro-estimation methods by introducing refined mechanism inference based on physical process decomposition and multiple physical constraints, and can provide more accurate, more physically deep, and more supporting fine-grained planning decisions for urban hydrological response characteristics.

[0101] The preferred embodiments of the present application are described in detail above, but the present application is not limited to the specific details of the above-described embodiments, and various equivalent transformations of the technical solutions of the present application can be made within the technical concept of the present application, and these equivalent transformations all belong to the protection scope of the present application.

Claims

1. A method for evaluating the synergistic optimization of urban-scale precipitation absorption and vegetation configuration, characterized in that, include: Acquire basic data of urban areas, perform data integration and rainfall event identification, and generate a rainfall event dataset; Based on rainfall event datasets, physical mechanisms are inferred and the hydrological response characteristics of the underlying surface are quantified. By combining the hydrological response characteristics of the underlying surface with the preset planning objectives, a robust optimization model is constructed and solved to determine the optimal vegetation configuration scheme. The optimal vegetation configuration scheme is back-substituted into the rainfall event dataset for closed-loop verification and effect evaluation, and an optimization evaluation report is generated. Quantitatively obtain the hydrological response characteristics of the underlying surface, including: For the total evapotranspiration process in the rainfall event dataset, the split point of the best-fitting two-stage exponential decay model is searched to preliminarily identify candidate switching moments; A physical gating mechanism is constructed, and multi-source physical quantities from the rainfall event dataset are invoked to jointly verify candidate switching moments; the physical gating mechanism includes energy consistency, wet surface dissipation, and geometric radiative sensitivity. The candidate switching moments that pass the joint verification will be confirmed as the stage switching moments from wet surface dominance to soil constraint. Based on the stage transition time, parameters are estimated and integrated for the subsequent soil-constrained stages to generate hydrological response characteristics of the underlying surface.

2. The method according to claim 1, characterized in that, Constructing physical gating and performing joint verification further includes: For the energy flux after the candidate switching time, examine the linear regression relationship between the sum of latent heat and sensible heat flux and the difference between net radiation and surface heat flux, requiring that the slope of the linear regression relationship is close to one and the intercept is close to zero. A wet surface dissipation index is constructed, which is a weighted average of the rate of change of surface temperature, the rate of change of surface albedo, and the shortwave radiation recovery ratio. The wet surface dissipation index is then used to determine whether it exceeds a preset threshold within a continuous period. During the transition period from shadow to sunlight, the response amplitudes of intercepted evaporation and impermeable surface evaporation to shortwave radiation are calculated separately, and the response amplitude of intercepted evaporation is required to be no higher than that of impermeable surface evaporation.

3. The method according to claim 1, characterized in that, Quantitatively obtaining the hydrological response characteristics of the underlying surface also includes: Based on shortwave radiation, surface temperature gradient and wind field data in the rainfall event dataset, the advection enhancement trigger window is identified that meets the conditions of clear sky, temperature difference along wind direction reaching a preset threshold and wind speed meeting the standard. Within the advection enhancement trigger window only, the non-negative enhancement term related to the product of the surface temperature horizontal gradient and the effective wind speed is quantified as the advection enhancement amount of impermeable surface evaporation. The advection enhancement is written back to the evaporation component of the impermeable surface to correct the evapotranspiration process, and the hydrological response characteristics of the underlying surface are generated based on the corrected evapotranspiration process.

4. The method according to claim 1, characterized in that, After confirming the phase transition time, the following is also included: During the time period from the end of the rain to the transition point, an optimization decomposition problem is established with the objective of minimizing the difference between total evapotranspiration and the sum of multiple evapotranspiration components. At least three types of physical constraints are imposed on the optimization decomposition problem: Energy conservation constraints limit the total latent heat flux to no more than the available energy. The time sequence constraint requires that the decay time constant of intercepted evaporation be less than the decay time constant of evaporation on the impermeable surface; Geometric sensitivity constraints limit the response amplitude of intercepted evaporation to shortwave radiation to not exceed that of impermeable surfaces. Solve the optimization decomposition problem with physical constraints to obtain the time series of each evapotranspiration component, and use it to generate the hydrological response characteristics of the underlying surface.

5. The method according to claim 1, characterized in that, Determining the optimal vegetation configuration includes: Based on the uncertainty information contained in the hydrological response characteristics of the underlying surface, a distribution uncertainty set containing multiple possible probability distributions is constructed using a sensitivity-weighted distance metric. Establish a partial bluff chance constraint, requiring that under the worst probability distribution within the uncertainty set, the probability that the total water storage increment meets the planning target is not lower than the preset confidence level; The opportunistic constraints of the pluripotent rods are transformed into solvable deterministic constraints, and then combined with engineering constraints for solving to determine the optimal vegetation configuration scheme.

6. The method according to claim 1, characterized in that, Before constructing the robust optimization model, the following is also included: Representative vegetation configuration values ​​are selected for candidate plots, and the physical mechanism inference process is invoked to generate mechanism anchor points containing the underlying surface hydrological response characteristics corresponding to each value. Based on the mechanism anchor points, a semi-convex substitution function is constructed at the plot level. The semi-convex substitution function is composed of the local linear approximation at each mechanism anchor point and the upper envelope of the quadratic shape penalty term. Semi-convex substitution functions are used to characterize the relationship between vegetation configuration and precipitation absorption capacity in robust optimization models.

7. The method according to claim 6, characterized in that, Constructing a semi-convex substitution function at the parcel level further includes: An upper bound constraint on the slope is imposed on the semi-convex substitution function to ensure that the precipitation absorption gain brought about by the increase in vegetation configuration does not exceed the physical reasonable upper limit. At the same time, a second-order upper bound constraint on the curvature is applied to the semi-convex substitution function to ensure that the response surface of the semi-convex substitution function is smooth and conforms to the physical law of diminishing marginal utility.

8. The method according to claim 1, characterized in that, Quantitative analysis of the underlying hydrological response characteristics can also be used to: Urban heat source separation and multi-scale smoothing were performed on the surface temperature data contained in the rainfall event dataset to obtain the time series of clean surface temperature. A continuous phase transition weighting function is constructed to characterize the gradual transition process from the wet surface-dominated stage to the soil-constrained stage. The specific form of the phase transition weight function is determined by a joint optimization objective function, which couples the phase transition rate, energy consistency, and surface resistance continuity driven by the dynamic characteristics of the clean surface temperature time series. Based on a defined phase transition weighting function, the evapotranspiration process is continuously decomposed and calculated to generate the hydrological response characteristics of the underlying surface.

9. The method according to claim 8, characterized in that, Obtaining the time series of clean surface temperatures includes: Establish an anthropogenic heat flux model that includes traffic intensity, building density, and diurnal periodic terms; By taking into account both the minimization of energy balance residuals and the spatial smoothness constraint of variational fitting objectives, the coefficients to be estimated for the artificial heat flux model are determined. The anthropogenic heat flux effects estimated by the anthropogenic heat flux model with determined coefficients are separated from the surface temperature data to obtain the clean surface temperature time series.

Images