Complex Slope Runoff Simulation Method Based on Multi-Source Data Fusion

Abstract

This invention discloses a method for simulating complex slope runoff based on multi-source data fusion. It includes: fusing UAV LiDAR, ground-penetrating radar, and remote sensing images, and considering physical coupling relationships such as vegetation root systems, to generate a high-precision and physically consistent initial surface parameter field; identifying micro-topographic depressions using graph cutting algorithms and constructing a dynamically evolving slope overflow network based on a depression-filling-overflow mechanism; establishing a dynamic feedback control mechanism of overflow-scour-infiltration: the scour intensity calculated based on overflow data is used to correct the soil saturated hydraulic conductivity in real time, and the corrected parameters are then fed back into an improved Green-Ampt infiltration model; and outputting the spatiotemporal distribution of slope runoff through multi-timescale algorithm coupling, which can trigger early warning and invert critical rainfall conditions. This invention can simulate the dynamic feedback process of the underlying surface, improve simulation reliability, and provide technical support for flash flood early warning.

Description

Technical Field

[0001] This invention relates to the field of slope runoff simulation methods, and in particular to a complex slope runoff simulation method based on multi-source data fusion. Background Technology

[0002] Slope runoff, a crucial link in the surface hydrological cycle, serves as a conduit between rainfall and watershed runoff. Its accurate simulation is vital for the rational use of water resources, early warning of flash floods, assessment of soil erosion, and prevention of non-point source pollution. Particularly in ecologically fragile areas such as the Loess Plateau in my country, soil erosion driven by slope runoff is one of the major geological hazards. Despite various mitigation measures, accurately predicting the generation and evolution of runoff under extreme rainfall conditions remains a significant challenge in disaster prevention and mitigation. Therefore, developing simulation methods capable of finely characterizing slope runoff generation and confluence processes under complex underlying surface conditions is a cutting-edge research topic in hydrology and disaster science.

[0003] Currently, slope runoff simulation research has developed multiple technical approaches. Among them, lumped or semi-distributed hydrological models, such as SWAT (a watershed-scale continuous distributed hydrological model) and HEC-HMS, have been widely used in large-scale watershed hydrological simulations by generalizing complex slopes into homogeneous units and using empirical or semi-empirical formulas to describe runoff generation and confluence processes. Distributed physical models based on the Saint-Venant equations or their simplified forms (such as the equations for diffused waves and kinematic waves) have also been developed, aiming to describe the physical details of water flow at a grid scale. In terms of data acquisition, emerging technologies such as UAV remote sensing and LiDAR are increasingly being used to acquire high-resolution topographic data. In the simulation of infiltration processes, classic infiltration theories such as the Green-Ampt model and the Horton model remain commonly used modules in various current models.

[0004] However, existing simulation methods still face several deep-seated technical bottlenecks when dealing with complex slopes exhibiting high spatial heterogeneity. These bottlenecks limit the physical realism and prediction accuracy of the simulation results. These problems mainly manifest in the static treatment of the underlying surface physical state, the oversimplification of micro-hydrological processes, and the fragmented application of the physical implications of multi-source data. Existing models generally treat key hydraulic parameters such as soil saturated hydraulic conductivity as static constants throughout the entire rainfall process, neglecting a physical feedback mechanism: the scouring effect of slope runoff (especially concentrated flow) alters the surface soil structure in real time, such as forming surface closures or washing away fine particles, leading to the exposure of large pores and dynamically changing its infiltration performance. The lack of dynamic coupling between overflow, scouring, and infiltration processes results in a continuously accumulating bias in the model's estimation of cumulative infiltration and total runoff during long-duration, high-intensity rainfall events. This problem is exacerbated by the simplification of the micro-topographical control effect. Traditional models, relying on coarse-resolution DEM data, cannot identify micro-depressions at the centimeter level and neglect the nonlinear runoff generation mechanism of depression filling-connection-overflow. This not only affects the accurate determination of runoff timing but also prevents the model from capturing the formation location and timing of the priority flow paths driving the aforementioned scour feedback mechanism. The inadequacy of physical process simulation is also related to the way front-end data is fused. Even with multi-source data, traditional methods often treat them as independent layers and simply overlay them, failing to consider the physical coupling relationships between data, such as the intrinsic influence of vegetation roots on ground-penetrating radar soil moisture signals. This results in physical inconsistencies in the generated initial surface parameter field, limiting the reliability of the simulation. Summary of the Invention

[0005] The purpose of this invention is to provide a method for simulating complex slope runoff based on multi-source data fusion, so as to solve the above-mentioned problems existing in the prior art.

[0006] According to one aspect of this application, a method for simulating complex slope runoff based on multi-source data fusion includes:

[0007] Acquire multi-source spatial data to generate an initial surface parameter field containing a digital elevation model and soil parameter distribution;

[0008] Based on the initial surface parameter field and rainfall data, a dynamic overflow network and infiltration model are constructed, and dynamic overflow data and dynamic infiltration data are calculated.

[0009] Based on dynamic overflow data and initial surface parameter field, a feedback control mechanism between overflow, scour and infiltration is established, and the slope runoff results are calculated by combining dynamic infiltration data.

[0010] Analyze slope runoff results and output runoff early warning information.

[0011] Beneficial effects: Through the above technical solutions, the present invention can simulate the dynamic feedback process of the underlying surface, improve the simulation reliability, and provide technical support for flash flood early warning. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the overall process of a complex slope runoff simulation method based on multi-source data fusion;

[0013] Figure 2 A schematic diagram illustrating the feedback control mechanism between overflow, scouring, and infiltration;

[0014] Figure 3 A schematic diagram illustrating the process of dynamically adjusting the soil hydraulic conductivity in the initial surface parameter field;

[0015] Figure 4 To perform subsequent infiltration and confluence calculations within this adaptive time step, a pre-defined execution sequence flowchart is shown. Detailed Implementation

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

[0017] The terms "first," "second," etc., used in the specification and accompanying drawings of this invention are used to distinguish different objects, not to describe a specific order. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, apparatus, or product comprising a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or products. The reference to "embodiment" herein means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of the invention. The appearance of this phrase in various places in the specification does not necessarily refer to the same embodiment, nor is it a separate or alternative embodiment mutually exclusive with other embodiments. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0018] This specification adopts a unified convention for symbols and data items. The main characters used in the embodiments of this invention are defined below:

[0019] ICP stands for IterativeClosestPoint, the iterative closest point algorithm.

[0020] Ks represents the saturated hydraulic conductivity of soil (m / s), which is the water conduction capacity of soil when it is fully saturated.

[0021] S represents soil suction (m), soil matrix potential, which reflects the soil's ability to retain water.

[0022] θi represents the initial soil moisture content (m 3 / m 3 (Soil volumetric moisture content before rainfall begins).

[0023] θs represents the saturated soil moisture content (m³). 3 / m 3 ), the volumetric water content when the soil is fully saturated.

[0024] Δθ represents the change in soil moisture content, Δθ = θs - θi.

[0025] NDVI stands for Normalized Difference Vegetation Index.

[0026] NDVI _max This represents the maximum NDVI value within the study area.

[0027] NDVI _min This represents the minimum NDVI value within the study area.

[0028] L _root This indicates the characteristic depth (m) of the vegetation root system.

[0029] ρ _root Root density (kg / m³) 3 ).

[0030] ρ _canopy This represents the canopy point cloud density (points / m²). 3 ).

[0031] σ _vegetation This represents the vegetation cover influence coefficient, which is dimensionless.

[0032] Hmax represents the maximum water storage depth (m) of the micro-depression.

[0033] V _storage Indicates the water storage volume of the micro-depression (m³) 3 ).

[0034] H _spill This indicates the overflow elevation (m), which is the water level elevation at which the depression begins to overflow.

[0035] H_spill_ij The elevation threshold (m) represents the overflow from depression i to depression j.

[0036] A _dry Indicates the area of ​​the non-waterlogged area (m²) 2 ).

[0037] A _pond Indicates the area of ​​the waterlogged region (m²) 2 ).

[0038] A _eff Indicates the effective infiltration area (m²) 2 ).

[0039] A _total Represents the total area of ​​the calculation unit (m²) 2 ).

[0040] h _pond Indicates the depth of the water (m).

[0041] h _ref This represents the reference depth (m), with a value of 0.5 × S.

[0042] h _i This represents the water level height (m) of the i-th depression.

[0043] i(t) represents the infiltration rate (m / s) at time t.

[0044] I(t) represents the cumulative infiltration volume (m).

[0045] p represents rainfall intensity (m / s).

[0046] θ represents the slope angle (degrees or radians).

[0047] t _pond Indicates the time (s) at which water begins to accumulate.

[0048] G(V,E) represents a graph structure, where V is the set of nodes (concave areas) and E is the set of edges (overflow paths).

[0049] C _ij The value represents the connectivity between depressions i and j, where 1 indicates connectivity and 0 indicates disconnection.

[0050] C(t) represents the dynamic connectivity matrix.

[0051] L _ij This represents the distance (m) between depressions i and j.

[0052] Q _overflow_ij The overflow flow rate (m) from depression i to j represents the overflow flow rate from depression i to j. 3 / s).

[0053] K_weir This represents the weir flow coefficient, with a value of 0.4, and is dimensionless.

[0054] W _ij This represents the overflow width (m) between depressions i and j.

[0055] φ(θ _ij ) represents the flow direction angle correction factor, φ(θ) _ij )=cos(θ _ij / 2).

[0056] θ _ij This represents the flow angle (in radians) from the depression i to j.

[0057] v _path This represents the average flow velocity (m / s) along the preferred flow path.

[0058] Q _total Represents the cumulative path flow (m) 3 ).

[0059] t _duration Indicates the duration of the flow (s).

[0060] E _path Indicates the flow path scour intensity index (Pa·m) 2 / s).

[0061] E _critical Indicates the critical erosion strength (Pa·m) 2 / s).

[0062] τ represents the actual bed shear stress (Pa).

[0063] τ _c This represents the critical shear stress (Pa).

[0064] τ _bed This represents the bed surface shear stress (Pa).

[0065] τ _critical This represents the critical shear stress (Pa) used for potential field calculations.

[0066] ρ represents the density of water, 1000 kg / m³. 3 .

[0067] g represents the acceleration due to gravity, 9.8 m / s². 2 .

[0068] R represents the hydraulic radius (m).

[0069] S _f This represents the friction slope, which is dimensionless.

[0070] Ks _initialThis represents the initial saturated hydraulic conductivity of the soil (m / s).

[0071] Ks _modified This represents the corrected saturated hydraulic conductivity of the soil (m / s).

[0072] Fine _loss This represents the fine particle loss rate, dimensionless, ranging from 0 to 0.3.

[0073] β _soil This represents the soil erosion sensitivity coefficient, which is dimensionless.

[0074] α represents the scouring sensitivity coefficient (in Method 3), which is dimensionless.

[0075] γ represents the porosity increase coefficient (in Method 3), which is dimensionless.

[0076] Δt _overflow This indicates the time step of the overflow process, 0.1 seconds.

[0077] Δt _infiltration This indicates the time step of the infiltration process, in 1 second.

[0078] Δt _erosion This indicates the time step of the flushing process, in 10-second increments.

[0079] Δt _adaptive This represents the adaptive time step (s).

[0080] h represents the water depth (m).

[0081] q represents the unit width flow rate (m) 2 / s).

[0082] Q _outlet (t) represents the flow rate process curve at the outlet section (m). 3 / s).

[0083] u _max This indicates the maximum flow velocity (m / s).

[0084] h _max This indicates the maximum water depth (m).

[0085] Q _peak Peak flow rate (m 3 / s).

[0086] u _threshold This indicates the flow velocity warning threshold, 2.0 m / s.

[0087] h _threshold This indicates the water depth warning threshold, 0.3m.

[0088] Q _threshold This indicates the traffic warning threshold, 5.0m.3 / s.

[0089] p _critical This indicates the critical rainfall intensity (m / s) required to trigger a warning.

[0090] t _critical This indicates the critical rainfall duration (in seconds) that triggered the warning.

[0091] p _min This represents the minimum critical rainfall intensity (m / s).

[0092] t _min This represents the minimum critical rainfall duration (s).

[0093] ε _raw This represents the original dielectric constant, which is dimensionless.

[0094] ε _corrected This represents the corrected dielectric constant, which is dimensionless.

[0095] z represents soil depth (m).

[0096] T _root This represents the root-soil transfer function, which is dimensionless.

[0097] M _validity This represents the cross-validation matrix, which is dimensionless.

[0098] σ _local This represents the local standard deviation.

[0099] r represents the correlation coefficient, which is dimensionless.

[0100] Φ _depression It represents a concave potential field and is dimensionless.

[0101] Φ represents the abbreviation for the depression potential value, which is dimensionless.

[0102] ▽ 2 h represents the terrain curvature (1 / m).

[0103] κ1, κ2, and κ3 represent the potential field weighting coefficients, which are 0.4, 0.3, and 0.3, respectively.

[0104] Example 1: This example describes the overall process of a complex slope runoff simulation method based on multi-source data fusion disclosed in this invention. Through a series of steps, a complete simulation chain is realized from multi-source data input to runoff early warning output, such as... Figure 1 As shown.

[0105] Step 1: Acquire multi-source spatial data to generate an initial surface parameter field containing a digital elevation model and soil parameter distribution.

[0106] Multi-source spatial data refers to various types of data describing surface characteristics acquired from different platforms or sensors. In this embodiment, this may include, but is not limited to, high-density point cloud data acquired using LiDAR on a UAV, slope detail data supplemented by a terrestrial 3D laser scanner, underground dielectric constant data obtained using ground-penetrating radar (GPR), and surface vegetation information acquired through multispectral or hyperspectral remote sensing imagery. The initial surface parameter field is one or more spatialized raster datasets that digitally represent the topography, soil physical properties, and vegetation cover of the simulated area at the initial moment.

[0107] Specifically, the process of generating the initial surface parameter field may include: registering and fusing UAV data with ground point cloud data to generate a digital elevation model (DEM) with centimeter-level resolution; processing ground-penetrating radar echo signals to invert and obtain spatial distribution maps of parameters such as soil saturated hydraulic conductivity and soil suction; and analyzing remote sensing images to extract the normalized vegetation index (NDVI) and obtain a vegetation cover distribution map.

[0108] In this embodiment, multi-source data fusion is employed to address the challenge of a single data source comprehensively and accurately characterizing the high heterogeneity of complex slopes. For example, using only UAV data may fail to capture ground details obscured by vegetation, while combining it with ground laser scanning can compensate for this deficiency. Fusion of multiple data sources can improve the accuracy and reliability of the initial surface parameter field.

[0109] Step 2: Based on the initial surface parameter field and rainfall data, construct a dynamic overflow network and infiltration model, and calculate the dynamic overflow data and dynamic infiltration data.

[0110] In this embodiment, this step identifies micro-topographic depressions on the slope and uses them as nodes for water collection and redistribution, simulating the infiltration and runoff generation process of rainfall at these nodes (depressions) and on the slope. The dynamic overflow network refers to a network composed of micro-depressions on the slope and the temporary hydraulic connections formed during rainfall, whose topology dynamically evolves with the change in water depth in each depression.

[0111] For example, this step identifies all micro-depressions and their geometric features (such as maximum volume and overflow elevation) based on a high-precision DEM. After the simulation begins, the water storage of each depression is tracked in real time based on the input rainfall data. When the water depth of a depression reaches its overflow elevation, it forms a hydraulic connection with another depression downstream, and an edge in the network is activated. Dynamic overflow data describes the overflow flow rate over time along each connection path in the network. The model calculates the infiltration process on the slope (including waterlogged and non-waterlogged areas) in parallel to obtain dynamic infiltration data, i.e., the spatial distribution of infiltration rate over time.

[0112] This dynamic modeling approach abandons the simplified treatment of uniform slopes in traditional models, and can accurately describe the physical process of filling depressions, connecting and overflowing dominated by micro-topography, more realistically reflecting the nonlinear characteristics and spatiotemporal evolution of runoff generation.

[0113] Step 3: Based on dynamic overflow data and the initial surface parameter field, establish a feedback control mechanism between overflow, scour and infiltration, and calculate the slope runoff results by combining dynamic infiltration data.

[0114] This step addresses the interaction between physical processes, a problem generally overlooked in existing models. The feedback regulation mechanism refers to simulating the scouring effect of surface runoff (overflow) on the soil surface. This scouring effect, in turn, changes the physical properties of the soil (such as hydraulic conductivity), affecting the subsequent infiltration process and influencing the total runoff on the slope, thus forming a closed-loop feedback.

[0115] Specifically, the model utilizes the dynamic overflow data obtained in the previous step (especially the preferential flow path data for high velocities) and combines it with the soil erosion resistance parameters in the initial surface parameter field to calculate the scour intensity distribution field characterizing the water erosion capacity. Based on this scour intensity, the soil hydraulic conductivity parameters at corresponding locations are dynamically adjusted. For example, in areas of intense scour, surface fine particles may be washed away, leading to increased macropore connectivity and thus increased hydraulic conductivity; or, under specific conditions, water compaction may also cause a decrease in hydraulic conductivity. The corrected hydraulic conductivity parameter field is then applied to the infiltration calculation in the next time step, realizing a feedback process.

[0116] By establishing the above feedback mechanism, this invention enables the model to simulate the self-evolution of the land surface, and can capture the changes in underlying surface characteristics caused by the interaction between water and soil during rainfall. It is particularly important for simulating runoff and soil erosion under long-term, high-intensity rainfall events, and the calculation results are more consistent with reality.

[0117] Step 4: Analyze the slope runoff results and output runoff early warning information.

[0118] This step transforms high-precision simulation results into disaster early warning information with practical application value. Slope runoff results typically include data on the changes in hydraulic elements such as velocity, water depth, and flow rate over time at any location on the slope and at the outlet section.

[0119] Specifically, this step presets a series of safety thresholds in the model, such as maximum allowable flow velocity, maximum safe water depth, or hazardous flow level. During the simulation, the system compares the calculated slope runoff results with these thresholds in real time. When the hydraulic elements at any monitoring point exceed their corresponding preset threshold, the system automatically triggers an early warning mechanism, generating and outputting warning information. Optionally, this warning information may include the warning level, the specific location exceeding the threshold, the time of exceeding the threshold, and the predicted impact range.

[0120] Example 2: This example describes in detail the feedback control mechanism described in Example 1. This mechanism improves the physical realism of the simulation by quantifying the interaction between overflow, flushing and infiltration.

[0121] The establishment of this feedback control mechanism, such as Figure 2 As shown, specifically, it includes the following steps:

[0122] Step 2.1: Analyze the dynamic overflow data and combine it with the initial surface parameter field to evaluate the scour intensity distribution field that characterizes the erosive capacity of the water flow.

[0123] Specifically, the scour intensity distribution field is assessed by: determining the bed shear stress based on the flow velocity and duration in the dynamic overflow data and in combination with the initial surface parameter field, and quantifying the scour intensity distribution field.

[0124] In this embodiment, the dynamic overflow data is preferably high-flow-rate paths and their hydraulic characteristic parameters, such as average flow velocity and duration, extracted by a priority flow path identification algorithm. The scour intensity distribution field, denoted by E(x,y), is a two-dimensional spatial field, where the value of each grid cell quantifies the intensity of water erosion at that location.

[0125] To obtain the scour intensity distribution field, it is necessary to calculate the actual bed shear stress τ within the channel. This calculation can be based on τ = ρ × g × R × S. _f Where τ is the actual bed shear stress (Pa); ρ is the density of water (kg / m³). 3 ), usually taken as 1000 kg / m 3 g is the acceleration due to gravity (m / s²) 2 (, usually taken as 9.8 m / s) 2 R is the hydraulic radius (m), which can be calculated based on water depth and wetted perimeter; S _f The friction slope can be approximated as the slope of a waterway.

[0126] After obtaining the actual shear stress τ, the critical soil shear stress τ is obtained from the initial surface parameter field. _c (Pa), and the average velocity v of the preferred flow path obtained from the dynamic overflow data. _path(m / s) and duration t _duration (s) can be used to calculate the scour intensity index E for each path. _path The preferred calculation formula is:

[0127] E _path =(τ-τ _c )×v _path 2 ×t _duration ; where, (τ-τ _c The term represents the net water flow force exceeding the soil's erosion resistance, v _path 2 The term reflects the kinetic energy of the water flow, t _duration The term reflects the cumulative effect of erosion. The E values ​​for each path are... _path The values ​​are assigned to the grid through which the flow passes, forming the scour intensity distribution field E(x,y).

[0128] Step 2.2: Based on the scour intensity distribution field, the soil hydraulic conductivity in the initial surface parameter field is used as the initial value to generate the corrected soil hydraulic conductivity parameter field.

[0129] Specifically, generating the corrected soil hydraulic conductivity parameter field includes: determining the fine particle loss rate based on the scour intensity distribution field, and dynamically adjusting the soil hydraulic conductivity in the initial surface parameter field by integrating the scour intensity distribution field and the fine particle loss rate to obtain the corrected soil hydraulic conductivity parameter field.

[0130] This step transforms the force calculated in the previous step (scour intensity) into a change in the underlying surface state (soil hydraulic conductivity). This process simulates two main physical effects: first, water scour causes surface soil compaction or fine particles to clog pores, reducing hydraulic conductivity; second, scour carries away fine particles, potentially exposing large pores or enhancing their connectivity, thus increasing hydraulic conductivity. In this embodiment, dynamic adjustment is achieved through a comprehensive correction formula, which includes a hydraulic conductivity reduction coefficient and a hydraulic conductivity increase coefficient: the hydraulic conductivity increase coefficient defines the impact of fine particle loss on pore connectivity. This can be calculated by determining the fine particle loss rate. _loss To achieve, Fine _loss =min(0.3,β _soil ×E _path 0.5 ), where β _soil E is the soil erosion sensitivity coefficient, a dimensionless parameter that reflects the sensitivity of different soil textures to erosion. _path This refers to the scouring intensity index calculated in the previous step; min(0.3,...) indicates that there is a physical upper limit to the loss rate, for example, not exceeding 30%. Therefore, the coefficient of increase in hydraulic conductivity can be expressed as (1+γ×Fine) _loss), where γ is the weighting coefficient for the increase in porosity, which is preferably 1.5 in this embodiment.

[0131] Furthermore, a coefficient representing the reduction in hydraulic conductivity, characterizing the effect of erosion on soil compaction, is defined. This coefficient can be expressed as (1-α×E) _path / E _duration ), where α is the weighting coefficient for scouring and compaction, preferably 0.2; E _duration This is the critical scouring strength threshold. When the actual scouring strength exceeds this value, the compaction effect becomes more pronounced.

[0132] Based on this, the initial soil saturated hydraulic conductivity Ks _initial Multiplying these two coefficients yields the corrected soil saturated hydraulic conductivity Ks. _modified The formula is:

[0133] Ks _modified =Ks _initial ×(1-0.2×E _path / E _duration )×(1+1.5×Fine _loss ).

[0134] It is understandable that, such as Figure 3 As shown, the steps for dynamically adjusting the soil hydraulic conductivity in the initial surface parameter field can also be as follows: based on the scour intensity distribution field, define a hydraulic conductivity reduction coefficient that characterizes the effect of scour on soil compaction; based on the fine particle loss rate, define a hydraulic conductivity increase coefficient that characterizes the effect of fine particle loss on pore connectivity; multiply the initial soil saturated hydraulic conductivity by the hydraulic conductivity reduction coefficient and the hydraulic conductivity increase coefficient to obtain the corrected soil saturated hydraulic conductivity.

[0135] Step 2.3: Apply the corrected soil hydraulic conductivity parameter field to the subsequent infiltration and runoff calculations.

[0136] The corrected soil hydraulic conductivity parameter field Ks is obtained _modified After (x,y,t), this parameter field will replace the initial Ks in the next computation time step of the model. _initial The infiltration rate is used as an input parameter for infiltration models (such as the Green-Ampt model) to calculate the infiltration rate at various locations on the slope. Correspondingly, changes in the infiltration rate directly affect surface runoff, influencing the subsequent slope runoff calculation results. Through this method, the alteration effect of scouring on the underlying surface is fed back into the hydrological process in real time, forming a dynamic, closed-loop simulation system.

[0137] Example 3: This example details the subsequent infiltration and confluence calculations based on the feedback mechanism. This example proposes a multi-timescale separation algorithm to address the contradiction between computational efficiency and stability caused by the significant differences in the characteristic timescales of different physical processes (overflow, infiltration, scour) in the model.

[0138] In this embodiment, the modified soil hydraulic conductivity parameter field is applied to subsequent infiltration and runoff calculations, executed through a multi-timescale separation algorithm. The physical processes in a slope runoff system exhibit varying response speeds: surface water flow (overflow) responds very rapidly, typically in the sub-second range; water infiltration into the soil is next, with a response in the second range; while soil property changes caused by water erosion are a cumulative effect, with the slowest response, in the tens of seconds range or even longer. Using a uniform, extremely small time step (determined by the fastest overflow process) to simultaneously calculate all processes would result in large, unnecessary computations. This algorithm achieves a balance between accuracy and efficiency by decoupling the computation frequencies of different processes.

[0139] Specifically, the algorithm includes:

[0140] Step 3.1: Differentiate and obtain the characteristic time scales corresponding to the overflow process, infiltration process, and the formation process of the scour intensity distribution field. The characteristic time scale refers to the typical time required to characterize a significant change in a physical process. In this embodiment, by analyzing the physical processes, characteristic time scales can be set for different processes:

[0141] characteristic time scale of overflow process (Δt) _overflow Because it is dominated by gravity and micro-topography, the response is extremely fast, and it is preferred to set it to sub-second level, such as 0.1 seconds.

[0142] Characteristic time scale of infiltration process (Δt) _infiltration This process involves the movement of water in a porous medium, with a moderate response speed, preferably set to the second level, such as 1 second.

[0143] Characteristic time scale of scouring process (Δt) _erosion This process involves the stripping and transport of soil particles, and is a relatively slow cumulative process, preferably set in the order of ten seconds, such as 10 seconds.

[0144] Step 3.2: Determine an adaptive time step based on the obtained multiple characteristic time scales, within which subsequent infiltration and confluence calculations will be performed. Adaptive time step (Δt) _adaptive Δt is the global time step used by the model for iterative calculations. To ensure the stability of the numerical calculations, this step size must be less than or equal to the smallest characteristic time scale in all the solutions. Therefore, in this embodiment, Δt _adaptive Set as Δt_overflow That is, 0.1 seconds. Furthermore, the algorithm does not apply to every Δt... _adaptive Instead of updating all processes internally, different update frequencies are used based on their respective characteristic time scales. For example, the overflow process updates at each Δt... _adaptive Calculations are performed every 0.1 seconds; the infiltration process can be calculated every 10 Δt seconds. _adaptive (That is, it updates once every 1 second); the flushing process can be done every 100 Δt. _adaptive It updates every 10 seconds, thus ensuring accurate capture of fast processes and avoiding redundant calculations for slow processes.

[0145] In some alternative implementations, Δt _adaptive It can also be dynamically variable. For example, the step size can be adjusted based on the rate of change of water depth or flow velocity in the model. When the water flow changes drastically, the step size can be automatically shortened, while when the water flow is stable, the step size can be appropriately extended to further improve computational efficiency.

[0146] Step 3.3: Within this adaptive time step, perform subsequent infiltration and confluence calculations, following a preset execution order. In Δt _adaptive Within this framework, the order of calculations for each process is crucial for ensuring physical causality. For example... Figure 4 As shown, in this embodiment, a fast-process-first serial computation order is preferably adopted, which is as follows:

[0147] The calculations corresponding to the overflow process are performed, specifically, based on the slope water depth distribution and infiltration loss at the previous moment, the hydrodynamic equations are solved, and the slope water depth, flow velocity and overflow between depressions at the current moment are updated.

[0148] Perform calculations corresponding to the infiltration process, specifically, based on the updated water depth distribution in the first step (especially in areas with water accumulation), calculate the infiltration rate and cumulative infiltration volume at the current moment.

[0149] The calculations corresponding to the formation and correction of the scour intensity distribution field and the soil hydraulic conductivity parameter field are performed. Specifically, based on the updated flow velocity and water depth in the first step, the scour intensity is calculated, and the soil hydraulic conductivity parameter field is updated according to the method in Example 2. The updated parameter field will be used for the next infiltration calculation cycle.

[0150] The above algorithm executes the physical process response in descending order of speed, ensuring that the calculation of the slow process is based on the latest physical state updated by the fast process, thus guaranteeing the integrity of the computational logic of the entire coupled system.

[0151] For example, the algorithm may further include integrated computation and result output. Specifically, at each time step, the updated parameters, particularly the dynamically changing infiltration rate i and the corrected soil hydraulic conductivity Ks, are integrated. _modified The affected hydraulic parameters are substituted into the governing equations of the slope confluence flow for solution. In this embodiment, the governing equation is preferably the diffusion wave equation dh / dt + dq / dx = pi. Where h is the water depth (m); t is time (s); and q is the unit width flow rate (m³ / s). 2 / s); p is the rainfall intensity (m / s); i is the infiltration rate (m / s). This equation can be solved discretely using numerical methods such as the finite difference method or the finite element method, and the output is the outlet cross-sectional flow process line Q as the result of slope runoff. _outlet (t) and the flow field distribution (u,v,h) on the slope space.

[0152] Example 4: This example describes in detail the specific process of generating the initial surface parameter field in Example 1.

[0153] Step 4.1: Fusing UAV and terrestrial laser point cloud data to generate a digital elevation model. Specifically, this step utilizes a UAV equipped with a LiDAR sensor to perform a flight scan of the entire study area, acquiring a large-scale airborne laser point cloud. For local areas with dense vegetation cover or complex terrain, a terrestrial 3D laser scanner is used for supplementary scanning to obtain a more accurate terrestrial point cloud. Using the Iterative Closest Point (ICP) algorithm or its variants, the two sets of point cloud data from different sources are registered to a unified coordinate system, forming a fused point cloud dataset. To ensure the accuracy of subsequent micro-topography analysis, the density of this fused point cloud is preferably up to 1000 points per square meter.

[0154] Building upon this, a progressive morphological filtering algorithm is employed to process the fused point cloud, separating ground points from non-ground points such as vegetation and buildings. Then, using only ground points, the discrete point cloud data is converted into a continuous raster surface through irregular triangular mesh interpolation, generating an original digital elevation model (DEM) with a spatial resolution of, for example, 5 cm. To eliminate high-frequency noise that may be generated during interpolation, median filtering is used to smooth the original DEM, yielding a centimeter-level DEM for simulation.

[0155] By integrating drone and ground scanning data, both acquisition efficiency and data accuracy are balanced, and the generated DEM can accurately depict the micro-topographic features that have a decisive impact on slope runoff paths.

[0156] Step 4.2: Process ground-penetrating radar echo signal data to obtain soil parameter distribution. This step non-invasively acquires the hydraulic properties of the topsoil on the slope. Specifically, a ground-penetrating radar (GPR) device is used to probe the slope along a pre-set survey line, collecting echo signal data of radar waves propagating and reflecting in the underground medium. The collected echo signals are processed (e.g., filtering, gain adjustment) to calculate the propagation speed of electromagnetic waves and the dielectric constant of the soil. Since there is a clear empirical relationship between the dielectric constant of the soil and its water content, the dielectric constant distribution map can be inverted into the distribution map of the initial soil water content θi based on this relationship. The correlation between soil texture data from field borehole sampling and indoor geotechnical test results can be combined to establish the relationship between soil texture and other hydraulic parameters (such as saturated hydraulic conductivity Ks and soil suction S). Using geostatistical methods such as Kriging interpolation, discrete survey line data and sampling point data are spatially interpolated to generate distribution maps of soil saturated hydraulic conductivity Ks, soil suction S, and initial water content θi covering the entire study area.

[0157] Step 4.3 involves analyzing multispectral remote sensing images to extract vegetation cover information as a component of the initial surface parameter field. This step quantifies the impact of vegetation on runoff processes. Specifically, multispectral remote sensing images of the study area (from UAVs or satellite platforms) are acquired and radiometrically calibrated and atmospherically corrected to eliminate errors introduced by the sensor itself and the atmosphere. After correction, the Normalized Difference Vegetation Index (NDVI) is calculated, which is an indicator that effectively reflects surface vegetation cover and growth vitality. By thresholding the NDVI images, a vegetation cover distribution map can be quickly extracted.

[0158] To further differentiate the varying impacts of different vegetation types on hydrological processes (such as root depth and canopy interception capacity), supervised classification methods can be employed. For example, a Support Vector Machine (SVM) classifier is preferred, using a small sample of vegetation types from a ground survey as training data to classify the entire image and generate a vegetation type distribution map. This distribution map can then be linked to a parameter library, assigning corresponding root density, canopy interception parameters, etc., to different vegetation types, serving as an integral part of the initial land surface parameter field.

[0159] Example 5: This example provides a preferred alternative, describing a method for generating surface parameter fields that considers the physical coupling relationships between multiple data sources. Example 4 treats topography, soil, and vegetation as independent layers for data processing, even though there are actually close physical connections between them (e.g., vegetation roots affect soil structure and dielectric properties). This example generates a more physically consistent and accurate initial parameter field by explicitly modeling the physical coupling relationships.

[0160] Step 5.1: Extract vegetation feature information from multi-source spatial data; construct the physical coupling relationship between vegetation feature information and data signals from different sources in the multi-source spatial data. Specifically, constructing the physical coupling relationship involves building a three-dimensional root-soil transfer function characterizing the influence of roots on soil electromagnetic properties based on vegetation feature information. Specifically, from the fused point cloud dataset, point clouds above a certain threshold (e.g., 0.2 meters) above the ground are separated as vegetation canopy point clouds. NDVI distribution maps are extracted from multispectral images. The three-dimensional root-soil transfer function T is constructed. _root This function quantifies the influence of vegetation roots on physical detection signals (in this case, ground-penetrating radar signals) at any point (x, y, z) in three-dimensional space. The function can be constructed as: T _root (x,y,z)=1+ρ _root ×(NDVI / NDVI _max ) 0.6 ×exp(-z / L _root ).in,

[0161] T _root It is a dimensionless root-soil transfer function;

[0162] (x,y) represents the spatial plane coordinates, and z represents the soil depth (m);

[0163] ρ _root Root density (kg / m²) 3 It can be estimated based on the NDVI value, for example:

[0164] ρ _root =0.3×(NDVI-NDVI _min ) / (NDVI _max -NDVI _min );

[0165] NDVI, NDVI _max NDVI _min These are the maximum and minimum normalized vegetation indices at the current point and in the study area, respectively.

[0166] L _root The characteristic root depth (m) for the vegetation type can be obtained from the vegetation type map.

[0167] The function shows that the effect of the root system increases with increasing root density and NDVI, and decreases exponentially with increasing depth z.

[0168] Step 5.2: Based on the physical coupling relationship, mutual correction and fusion of data signals from different sources are performed to generate an initial surface parameter field. Specifically, the mutual correction and fusion of data signals from different sources involves applying a three-dimensional root-soil transfer function to perform root attenuation correction on the ground-penetrating radar echo signal data from the multi-source spatial data, and then inverting the soil parameter distribution based on the corrected ground-penetrating radar echo signal data.

[0169] The original dielectric constant ε measured by ground penetrating radar _raw It is simultaneously affected by soil moisture and root systems. To isolate the influence of root systems and obtain more accurate soil moisture information, the T model constructed in the previous step is used. _root Function with respect to ε _raw Perform correction. The correction formula can be: ε _corrected =ε _raw / (1+0.2×T _root ×log(1+ρ _root )). Among them, ε _corrected This is the corrected dielectric constant.

[0170] Specifically, T _root and ρ _root The larger the area, the stronger the attenuation and scattering effect of the root system on the radar signal, thus requiring a larger correction value. The obtained ε _corrected By substituting the empirical relationship between dielectric constant and water content mentioned in Example 4, a more accurate soil water content profile, corrected for root influence, can be obtained through inversion, and other soil parameter fields can be generated accordingly.

[0171] Step 5.3 verifies the physical consistency of the fusion results. This includes: extracting topographic morphology indicators and soil moisture distribution indicators from the initial surface parameter field; calculating the correlation coefficient between the topographic morphology indicators and the soil moisture distribution indicators to quantify the physical consistency between surface morphology and soil moisture distribution; and adaptively determining the final values ​​of the initial surface parameter field based on the correlation coefficient.

[0172] To ensure the physical consistency of the fused parameter field, this step introduces a verification and adaptive decision-making mechanism. In nature, topographic features and soil moisture distribution are often correlated (e.g., soil moisture content in depressions is usually higher than that at the top of slopes). This step utilizes this prior knowledge to verify the rationality of the data fusion. Specifically, topographic features (such as surface curvature ▽ calculated using a DEM) are extracted from the final generated parameter field. 2 (DEM) and soil moisture distribution indices (such as root-corrected topsoil moisture content). Calculate the spatial correlation coefficient r between these two indices across the entire study area. Make decisions based on this correlation coefficient:

[0173] If r is greater than the preset physical consistency threshold (e.g., r > 0.6), the coupling relationship between topographic and soil moisture data is considered reasonable, the data reliability is high, and further weighted fusion can be performed to optimize parameters, such as Ks. _fused =Ks _gpr ×(1+0.3×T _root ×σ _vegetation ).

[0174] Otherwise, if there is a significant physical inconsistency between the two, it may mean that one of the data sources has a large error. In this case, the system will retain the result from the single data source with higher credibility, thus avoiding erroneous fusion and ensuring the robustness of the final output parameter field.

[0175] Optionally, as a more refined quality control method, an anomaly feature identification step can also be introduced. For example, by calculating the cross-validation matrix M... _validity This matrix can assess the fusion confidence of local regions based on the consistency of parameters derived from different data sources (such as LiDAR and GPR) in spatial gradients. When M _validity When the data is below the threshold, the area can be marked on the data anomaly map and its original anomaly features can be preserved without fusion processing. This is of great significance for protecting some special but real geological or geomorphological anomalies.

[0176] According to one aspect of this application, the steps of vegetation cover information extraction and classification can also be a vegetation-mediated multi-source data adaptive fusion process, specifically as follows:

[0177] Points 0.2 meters above the ground were extracted from the fused point cloud dataset as vegetation canopy point clouds, and the local point density ρ of each canopy point was calculated. _canopy Extract the NDVI values ​​at corresponding locations from multispectral remote sensing images and establish a correlation function f between point density and NDVI. _correlation =ρ _canopy ×(NDVI / NDVI _max ) 0.8 Generate a canopy-spectral correlation coefficient field.

[0178] Based on the vegetation type distribution map, query the standard root depth L of different vegetation types. _root Calculate root density ρ using NDVI distribution map _root =0.3×(NDVI-NDVI _min ) / (NDVI _max -NDVI _min Construct a three-dimensional root-soil transfer function:

[0179] T _root (x,y,z)=1+ρ _root ×(NDVI / NDVI_max ) 0.6 ×exp(-z / L _root Output the root soil transfer function field T _root .

[0180] Read the ground-penetrating radar echo signal data and use the root-soil transfer function field T _root Differential correction is applied to radar signals at different depths, and the corrected dielectric constant ε _corrected =ε _raw / (1+0.2×T_root×log(1+ρ_root)), the root-corrected soil moisture profile is obtained by inverting the corrected dielectric constant.

[0181] The canopy-spectral correlation coefficient field was used as a weight to adjust the confidence of LiDAR data, resulting in a vegetation-corrected DEM, and the surface curvature ▽ was calculated. 2 The correlation between the DEM and the root-corrected soil moisture profile was analyzed, and a weighted fusion Ks algorithm was performed when the correlation coefficient r > 0.6. _fused =Ks _gpr ×(1+0.3×T _root ×σ _vegetation Otherwise, retain a single, more reliable data source and output a physically consistent soil parameter field.

[0182] Calculate the local standard deviation σ at each location in a physically consistent soil parameter field. _local When |P _current -P _mean |>2σ _local Time markers are used to identify outliers, and a cross-validation matrix is ​​constructed:

[0183] M _validity [i,j]=correlation(dP _lidar / dx,dP _gpr / dx)×T _root When M _validity When the value is less than 0.3, the region is marked on the data anomaly marking map, and the original anomaly features are preserved without fusion processing.

[0184] Example 6: This example describes in detail the process of identifying micro-topographic depressions as network nodes from a high-precision digital elevation model (DEM) before constructing a dynamic overflow network.

[0185] The runoff simulation method of the present invention includes a step of identifying micro-topographic depressions, as follows:

[0186] Step 6.1: Based on the digital elevation model (DEM) in the initial surface parameter field, identify local confluence points as depression centers. Local confluence points, often referred to as depressions or sinks in hydrology, are raster cells in the DEM raster data whose elevation values ​​are lower than all their neighboring raster cells. In this embodiment, this identification process is implemented using a standard DEM hydrological analysis algorithm. Specifically, the flow direction matrix of the DEM can be calculated, for example using the D8 algorithm, which determines the flow direction for each raster cell, pointing towards the raster with the steepest slope in its eight neighborhoods. Based on this, by analyzing the flow direction matrix, all raster cells with no outflow (i.e., no flow towards neighboring raster cells), or where the flow from all neighboring raster cells points towards the raster cell, are identified as local confluence points. These points constitute the seeds for subsequent regional growth.

[0187] Step 6.2: Starting from the center of the depression, apply the seed point region growing algorithm to determine the catchment area of ​​each depression. The seed point region growing algorithm is an image segmentation technique used in this embodiment to determine the precise spatial extent of each depression unit, i.e., the boundary of its catchment area. The principle of this algorithm is similar to the process of filling a depression with water in the physical world. Specifically, the algorithm starts from a seed point (local catchment point) and iteratively merges grid cells in its neighborhood that meet specific conditions into the current region. The condition here is that the elevation of the neighboring grid cells must be lower than the final overflow elevation of the depression. In other words, the algorithm continuously searches for the point with the lowest elevation on the boundary of the current region and merges all adjacent grid cells with an elevation lower than that point into the region until the elevation of all points on the boundary of the region is higher than the lowest point on that boundary. At this point, the region stops growing, and the set of grid cells it contains constitutes the catchment area of ​​the depression, while the lowest point on the boundary is the overflow outlet of the depression. Repeating this process for all identified seed points allows for the segmentation of all micro-depressions on the slope.

[0188] Step 6.3: Calculate and extract the geometric feature parameters of each depression based on its catchment area to form a micro-depression feature database for subsequent steps. After determining the catchment area of ​​each depression, its geometric and hydrological characteristics need to be quantified for subsequent water storage and overflow calculations. The micro-depression feature database is a structured dataset that stores the identifier of each depression and its corresponding feature parameters. The main parameters include:

[0189] Overflow elevation H _spill (m): Calculate and determine the elevation value of the lowest point on the boundary of the catchment area.

[0190] Maximum water storage depth Hmax (m): Calculates the overflow elevation H _spill The difference between the elevation of the lowest point inside the depression (i.e., the seed point) and the elevation of the seed point.

[0191] Water storage volume V _storage (m3 ): This is calculated by volume integration over all grids within the catchment area. Specifically, the overflow elevation H of each grid is calculated. _spill Multiply the difference between the grid's elevation and its own elevation (if positive) by the grid area, and then sum the calculation results for all grids.

[0192] After storing these parameters in the database, the necessary basic data was provided for constructing the dynamic overflow network and simulating the depression filling process in Example 8.

[0193] Example 7: This example, as a preferred alternative to Example 6, provides a micro-depression identification method based on physical meaning and with higher computational efficiency. Traditional depression filling or region growing algorithms involve enormous computational costs when processing high-resolution (e.g., centimeter-level) DEMs, and their segmentation is based purely on geometric elevation, failing to incorporate other physical factors affecting depression formation and development. This example solves the above problems by introducing a depression potential field and a graph cutting algorithm.

[0194] In a preferred embodiment, the step of identifying micro-topographic depressions can also be implemented using a graph-based segmentation method, including:

[0195] Step 7.1: Calculate the depression potential field based on the initial surface parameter field. Depression potential field Φ _depression This is a two-dimensional spatial field of the same size as a DEM, where the value of each grid cell represents the probability or potential of the location to become or belong to a water-bearing depression. The construction of this potential field integrates factors of topography, soil, and hydrodynamics, not just elevation. Specifically, this step determines topographic factors characterizing topographic geometry, soil property factors characterizing soil infiltration capacity, and hydrodynamic factors characterizing water erosion potential from the initial surface parameter field. In this embodiment, these factors are determined as follows:

[0196] Topographic factors: Calculation of topographic curvature based on digital elevation model ▽ 2 h, and take its square (▽) 2 h) 2 The concave region typically has a large positive curvature (bowl-shaped in shape), so the larger the value, the higher the concave potential.

[0197] Soil property factors: Infiltration resistance 1 / Ks is determined based on the distribution of soil parameters, where Ks is the saturated hydraulic conductivity of the soil. Areas with poor infiltration capacity (i.e., large 1 / Ks values) are more likely to accumulate water on the surface after rainfall, thus having a higher potential for depression.

[0198] Hydrodynamic factors: used to estimate the bed shear stress τ _bed and compare it with the critical shear stress τ _durationComparison. Areas with low bed shear stress have weak water flow energy, are less prone to erosion, and are more likely to become areas of siltation and depression development. This factor can be expressed as exp(-τ _bed / τ _duration ).

[0199] By weighting and combining the three factors, a concave potential field is generated:

[0200] Φ _depression =κ1·(▽ 2 h) 2 +κ2·(1 / Ks)+κ3·exp(-τ _bed / τ _duration ), where κ1, κ2, and κ3 are weighting coefficients, for example, 0.4, 0.3, and 0.3 respectively.

[0201] Step 7.2: Apply a graph cutting algorithm to segment the digital elevation model to identify concave regions. Graph cutting is an image segmentation technique that finds the minimum cut in a graph, dividing the graph's nodes into two disjoint sets. In this embodiment, the entire DEM is constructed as a graph G(V,E), where each raster is a node V. The algorithm first defines the source point S (e.g., all Φ...). _depression Nodes with values ​​greater than a certain high threshold) and sink T (e.g., all Φ _depression (Nodes whose values ​​are less than a certain low threshold). By solving the maximum flow minimum cut, the algorithm can find the optimal boundary, separating all nodes connected to the source node S (determined as concave regions) from the nodes connected to the sink node T (determined as non-concave regions).

[0202] Furthermore, this step preferably employs a multi-resolution strategy to improve computational efficiency: the centimeter-level DEM is downsampled to a lower resolution representation (e.g., 20cm resolution), and preliminary graph cutting segmentation is performed on this coarse grid to quickly obtain the preliminary location and approximate outline of the concave region. Based on the preliminary location of the concave region, a boundary refinement region is defined along its outline. Only within the boundary refinement region, the graph cutting algorithm is performed again on the high-resolution representation of the digital elevation model (i.e., the original 5cm resolution) to determine the precise boundary of the concave region at a lower computational cost.

[0203] Optionally, to capture the dynamic evolution of the land surface during rainfall, this step can also introduce a dynamic update mechanism. For example, a 5-minute update cycle can be set to periodically recalculate the potential field and perform map cutting during the simulation to identify micro-depressions that are newly formed or disappear due to erosion or siltation.

[0204] Step 7.3: Calculate and extract the geometric feature parameters of each depression based on the depression region to form a micro-depression feature database for subsequent steps. This step is similar to step 6.3 in Example 6; after determining the precise depression boundary through graph cutting, the overflow elevation H of each depression is also calculated. _spill Maximum water storage depth Hmax and water storage volume V _storage And store it in the database.

[0205] Example 8: This example describes in detail the two main destinations of rainwater that are processed in parallel: water that seeps into the ground (infiltration) and water that flows to the surface (overflow).

[0206] The first part involves calculating dynamic infiltration data. This calculation includes the following steps:

[0207] Step 8.1: Dynamically determine the effective infiltration area. Traditional infiltration models often assume uniform infiltration throughout the entire area, which does not match reality when surface water is present. This embodiment introduces an effective infiltration area A. _eff The concept of infiltration is that the infiltration capacity per unit area is enhanced in waterlogged areas due to the presence of water head pressure. Specifically, for areas that have accumulated water, the effective infiltration area is determined based on the actual waterlogged area and an enhancement calculation based on the water depth. In this embodiment, the effective infiltration area A... _eff Calculate according to the following formula:

[0208] A _eff =A _dry +A _pond ×(1+h _pond / h _ref ); where A _dry To calculate the area of ​​the unit without water accumulation (m²) 2 A _pond To calculate the area of ​​water accumulated within the unit (m²) 2 );h _pond h represents the average water depth (m) of the flooded area. _ref The reference depth (m) related to soil suction is used to calibrate the enhancing effect of water depth on infiltration. In this embodiment, h _ref The preferred value is 0.5 times the soil suction force S.

[0209] Step 8.2: Generate dynamic infiltration data based on the effective infiltration area. The A calculated in the previous step... _eff Integrating this into the Green-Ampt infiltration model, an improved formula for calculating the infiltration rate i(t) is formed:

[0210] i(t) = Ks × (A _eff (t) / A _total)×[1+(S×Δθ) / (I(t)+S×Δθ×h _pond / h _ref )];

[0211] Where i(t) is the infiltration rate at time t (m / s); Ks is the saturated hydraulic conductivity of the soil (m / s); A _total To calculate the total area of ​​the unit (m²) 2 S represents soil suction (m); Δθ represents the difference between the initial soil moisture content and the saturated soil moisture content; I(t) represents the cumulative infiltration (m). To accurately simulate the transition from no water accumulation to water accumulation, the model also determines the starting time t of water accumulation by comparing the effective rainfall intensity p×cosθ with the instantaneous infiltration capacity i(t). _pond Using this as a boundary, different infiltration calculation modes are adopted to generate an infiltration rate field that covers the entire rainfall duration and is dynamically changing in space and time, i.e., dynamic infiltration data.

[0212] The second part is the calculation of dynamic overflow data, which includes the following steps:

[0213] Step 8.3: Construct a dynamic connectivity matrix that evolves over time. This step, based on the micro-depression database obtained in Example 6 or 7, analyzes the DEM using the steepest descent method to determine all potential overflow paths between depressions and construct a static topology network.

[0214] Based on this, the model tracks the change in water depth h of each depression caused by rainfall data in real time at each time step t. _i (t). By comparing the real-time water depth h of the depression. _i (t) and its overflow elevation H _spill_ij To determine its connectivity with the downstream depression j. Based on this, a dynamic connectivity matrix C(t) is constructed. If h _i (t)>H _spill_ij Then matrix element C _ij (t) = 1 (indicating connectivity), otherwise 0.

[0215] Step 8.4 uses the dynamic connectivity matrix as the basis for calculating dynamic overflow data. For any connected pair of depressions (i,j) (i.e., C) determined by the dynamic connectivity matrix at any time, _ij (t)=1), and the overflow flow rate Q is calculated using the broad-crested weir formula. _overflow_ij The application of this formula is based on the real-time water depth h of the overflow depression. _i (t) and its overflow elevation H _spill_ij The determined head difference, and based on the overflow width W provided by the micro-indentation feature database. _ij In this embodiment, the formula is preferably:

[0216] Q _overflow_ij =C _ij ×K _weir ×W _ij ×(h _i -H _spill_ij ) 1.5 ×φ(θ _ij ); where K _weir The weir flow coefficient is preferably 0.4; W _ij The overflow width (m) between depressions i and j.

[0217] Optionally, the application of this formula further includes: based on the flow direction angle θ between connected pairs of depressions. _ij Calculate the flow direction angle correction factor φ(θ) used to correct for the influence of overflow direction. _ij This correction factor is then applied to the broad-crested weir formula. In this embodiment, the preferred correction factor is φ(θ). _ij )=cos(θ _ij / 2), used to simulate energy loss caused by the inconsistency between the overflow direction and the downstream slope.

[0218] The overflow flow rate of all connected depression pairs is calculated to form the dynamic overflow data at that moment.

[0219] Step 8.5: Identify the preferred flow path and extract its features. Specifically, this embodiment also includes identifying, based on dynamic overflow data, situations where the flow rate continuously exceeds a preset threshold (e.g., 0.001m). 3 Overflow paths with a velocity of v / s are identified as priority flow paths. For each identified priority flow path, its hydraulic characteristic parameters, including the average velocity v, are calculated. _path Cumulative traffic Q _total and duration t _duration The algorithm generates a priority flow path distribution map. This map visually displays the main collection channels of slope runoff, and its extracted feature parameters will serve as direct inputs for calculating the scour intensity distribution field in Example 2, closely linking the water collection process with the soil erosion process.

[0220] Example 9: This example describes how the slope runoff simulation results obtained in Example 3 are transformed into early warning information with practical disaster prevention and mitigation significance, and further reverse-engineering the critical rainfall conditions that trigger the early warning.

[0221] Step 9.1 involves comparing hydraulic element thresholds and triggering early warnings. This step enables real-time assessment of runoff risk. Specifically, the system uses the simulation results of slope runoff, namely the spatial flow field distribution (u, v, h) and the outlet cross-section flow process line Q, to determine the risk. _outlet In (t), key hydraulic risk indicators are extracted in real time. In this embodiment, these indicators are preferably: the maximum flow velocity u in the study area. _maxMaximum water depth h _max And the peak flow rate Q at the outlet section _peak The system pre-sets safety thresholds corresponding to the disaster-bearing capacity of the study area. These thresholds can be determined based on local historical disaster data, geographical information, or engineering design standards. For example, a flow velocity warning threshold u can be set. _threshold =2.0 m / s (This flow velocity may pose a threat to personnel safety or begin to scour small structures); water depth warning threshold h _threshold =0.3m (this water depth may cause road traffic disruptions or flooding of rural houses); flow warning threshold Q _threshold =5.0m 3 / s (This flow rate may exceed the safe discharge capacity of the downstream river channel).

[0222] At each time step in the simulation process, the system will extract u in real time. _max h _max Q _peak Compare with the preset thresholds mentioned above. When the condition is met:

[0223] (u _max >u _threshold OR(h) _max h _threshold OR(Q) _peak Q _threshold If any of the following conditions are met, the system will determine that the risk level has exceeded the safe range and generate a warning trigger flag.

[0224] In some alternative implementations, multiple threshold levels (e.g., attention level, warning level, danger level) can be set to achieve fine-grained classification of warning levels. Warning thresholds can also vary spatially, for example, setting stricter thresholds in areas with important protected objects (such as villages or roads).

[0225] Step 9.2 involves inverting the critical rainfall conditions. Specifically, simply knowing that current rainfall has caused a hazard is insufficient; for weather forecasting and emergency management, a more important question is how much rainfall is needed to trigger a hazard. This step addresses this issue through model inversion. When the warning in Step 9.1 is triggered, the system records the rainfall intensity p of the rainfall event that triggered the warning. _duration and rainfall duration t _duration The system automatically initiates the inversion calculation program to find the minimum rainfall combination (minimum critical rainfall intensity p) that can just trigger an early warning. _min and minimum critical rainfall duration t _min ).

[0226] In this embodiment, the inversion process is preferably implemented using the bisection method to solve for p under a fixed duration. _minFor example: the algorithm will be in [0, p _duration Perform an iterative search within the interval. In the first iteration, take the midpoint of the interval, 0.5*p. _duration As a new input for rainfall intensity, a full slope runoff simulation is run. If this simulation triggers an alert, it indicates that p _min Located in [0, 0.5*p _duration [Interval]; if not triggered, it means p _min Located at (0.5*p _duration ,p _duration The search interval will converge rapidly by repeatedly performing bisections and simulations on the given interval until a sufficiently accurate p is found. _min Value. For rainfall duration t _min A similar method can be used to solve this problem.

[0227] Alternatively, for more complex rainfall events (such as rainfall pattern changes), the inversion calculation can employ more advanced optimization algorithms, such as simulated annealing or genetic algorithms, to find a series of rainfall event curves that can trigger an early warning.

[0228] Step 9.3 involves generating and outputting early warning information. This step integrates all analysis results into a clear and usable early warning product. Specifically, the system will integrate the early warning trigger identifier sequence and the critical rainfall intensity p obtained in step 9.2. _min and critical rainfall duration t _min The system also collects the spatial location information of the hydraulic elements exceeding the limit when the warning is triggered. Based on this, the system generates a comprehensive warning report. This report may include the following:

[0229] Warning level: Determined based on the threshold level exceeded.

[0230] Warning trigger time and location: Clearly indicate the location and time when the risk first appears.

[0231] Expected impact area: Plotted based on the area covered by the excess water depth and flow velocity.

[0232] Critical rainfall conditions: providing forecasting departments with clear thresholds for disaster-causing rainfall.

[0233] All key parameters of early warning events, including rainfall events, simulation results, and inversion thresholds, will be stored in a historical database for subsequent model validation, calibration, and training of AI-assisted forecast models, forming a continuously optimized closed-loop early warning system.

[0234] Example 10: This example provides a specific numerical calculation case to describe the specific application of the overflow-scour-infiltration feedback control mechanism and multi-timescale coupled calculation.

[0235] Assume that at a certain adaptive time step Δt _adaptive Within 0.1 seconds, the model identified the preferred flow path using the method described in Example 8. This case demonstrates the calculation of scour intensity and feedback correction of soil hydraulic conductivity for a single computational unit along this path. The input parameters for this computational unit are set as follows:

[0236] Water flow and soil properties: average velocity v along the preferred flow path _path =1.5m / s; water flow duration t _duration (This is the cumulative time) = 120s; hydraulic radius R = 0.05m; friction slope S _f =0.03; initial soil saturated hydraulic conductivity Ks _initial =5.0×10 -5 m / s; Critical soil shear stress τ _c =1.2Pa; Soil erosion sensitivity coefficient β _soil =0.1. Model and physical constants: Critical scour strength E _duration =500Pa·m 2 / s; the density of water ρ = 1000 kg / m³ 3 The acceleration due to gravity is g = 9.8 m / s². 2 .

[0237] Step A: Calculate the actual bed shear stress τ. Specifically, according to Example 2, calculate the actual force exerted by the water flow on the riverbed. The formula used is: τ = ρ × g × R × S _f Substituting the parameter: τ = 1000 kg / m 3 ×9.8m / s 2 ×0.05m×0.03=14.7Pa.

[0238] Step B: Evaluate the flow path scour intensity E _path Specifically, the overall erosive capacity of water flow is quantified using the formula:

[0239] E _path =(τ-τ _c )×v _path 2 ×t _duration ;

[0240] Substitute parameter: E _path =(14.7Pa-1.2Pa)×(1.5m / s) 2 ×120s=13.5×2.25×120=3645Pa·m 2 / s.

[0241] Step C: Calculate the fine particle loss rate. _lossSpecifically, the loss of fine soil particles is estimated based on the erosion intensity using the following formula:

[0242] Fine _loss =min(0.3,β _soil ×E _path 0.5 );

[0243] Substitute parameters: Fine _loss =min(0.3,0.1×(3645) 0.5 =min(0.3,0.1×60.37)=min(0.3,6.037)=0.3. It can be seen that due to the high scouring intensity, the calculated fine particle loss rate has reached the set physical upper limit of 0.3.

[0244] Step D: Calculate the corrected soil saturated hydraulic conductivity Ks _modified Specifically, the initial hydraulic conductivity is dynamically adjusted based on the scouring intensity and the fine particle loss rate. The formula used is:

[0245] Ks _modified =Ks _initial ×(1-0.2×E _path / E _duration )×(1+1.5×Fine _loss );

[0246] Substitute parameter: Ks _modified =(5.0×10 -5 )×(1-0.2×3645 / 500)×(1+1.5×0.3); Calculate the coefficients within the parentheses: Hydraulic conductivity reduction coefficient = 1-0.2×3645 / 500 = 1-1.458 = -0.458; Hydraulic conductivity increase coefficient = 1+1.5×0.3 = 1.45. It is understandable that the calculated result of the hydraulic conductivity reduction coefficient is negative, which is physically impossible. Therefore, in the specific implementation of this invention, a lower limit of 0 is set for this coefficient, that is, when (1-0.2×E)×(1-0.2×3645 / 500)×(1+1.5×0.3), the coefficient is reduced to -0.458. _path / E _duration When Ks < 0, the value is 0. This represents extreme erosion that has caused the surface to become completely sealed, reducing its water conductivity to zero. Therefore, the corrected calculation is: Ks _modified =(5.0×10 -5 )×0×1.45=0m / s.

[0247] Based on the above calculations, the soil saturated hydraulic conductivity of the computational unit located on this preferred flow path, due to the intense water flow erosion, decreases from the initial 5.0 × 10⁻⁶ at the current moment. -5 m / s was dynamically corrected to 0 m / s. The newly calculated Ks _modifiedThe initial value is replaced by a new value, which serves as the input parameter for the unit in the next infiltration calculation cycle (e.g., after 1 second). That is, in the following time, all rainwater falling into the unit will be unable to infiltrate and will be 100% converted into surface runoff. This process demonstrates the overflow-scouring-infiltration feedback control mechanism of this invention: high-intensity overflow leads to severe scouring, severe scouring alters the soil's infiltration properties, and the change in infiltration properties affects future overflow processes, achieving a complete and physically consistent dynamic coupling.

[0248] This invention relates to several computational processes that differ from existing technologies, specifically including:

[0249] Calculate the effective infiltration area enhancement: A _eff =A _dry +A _pond ×(1+h _pond / h _ref This is used to dynamically quantify the enhancing effect of water depth on infiltration area.

[0250] Calculate the improved Green-Ampt infiltration rate:

[0251] i(t) = Ks × A _eff (t) / A _total ×[1+(S×Δθ) / (I(t)+S×Δθ×h _pond / h _ref This couples the dynamic effective infiltration area with the classical infiltration model.

[0252] Calculate the potential field of the depression: Φ _depression =0.4×(▽ 2 h) 2 +0.3×(1 / Ks)+0.3×exp(-τ _bed / τ _duration This method integrates topographic, soil, and hydrodynamic factors to provide a basis for map segmentation.

[0253] Calculate the cascade overflow flow rate: Q _overflow_ij =C _ij ×K _weir ×W _ij ×(h _i -H _spill_ij ) 1.5 ×φ(θ _ij Based on the broad-crested weir formula and with the introduction of angle correction, the overflow between depressions is accurately calculated.

[0254] Evaluate flow path scour intensity: _path =(τ-τ _c )×v _path 2 ×t _durationTaking into account shear force, flow velocity, and duration, the erosive power of water flow is quantified.

[0255] Corrected soil hydraulic conductivity dynamics: Ks _modified =Ks _initial ×(1-0.2×E _path / E _duration )×(1+1.5×Fine _loss Establish a physical feedback mechanism between scour intensity and soil hydraulic conductivity.

[0256] The slope runoff integration is calculated as: dh / dt+dq / dx=pi. The diffusion wave equation is used to integrate all sub-processes and solve for the final slope runoff result.

[0257] This invention details a dynamic feedback control mechanism for overflow-scour-infiltration. This mechanism calculates the scour intensity of runoff in real time and dynamically updates the soil saturated hydraulic conductivity parameter accordingly, enabling the model to capture the evolution of the underlying surface physical state due to soil-water interaction during rainfall. This allows the simulation process to reflect physical reality rather than being based on a static, unchanging surface, improving the accuracy of infiltration and total runoff estimation during long-duration heavy rainfall events and solving the problem of existing models treating soil hydraulic parameters as static constants.

[0258] This invention utilizes multi-source data to generate centimeter-level DEMs and employs advanced algorithms such as graph cutting based on physical potential fields to accurately identify micro-depressions. By constructing a dynamic connectivity matrix to explicitly simulate the nonlinear process of depression filling-connection-overflow, it reproduces the real physical mechanism of surface runoff generation. This allows for accurate prediction of runoff initiation time, spatial location, and the formation of preferential flow paths driving erosion, providing physical input for establishing the aforementioned feedback mechanism. It also addresses the problem of oversimplification in the regulation of micro-topography.

[0259] This invention proposes a vegetation-mediated optimization scheme. By constructing physical models such as the three-dimensional root-soil transfer function, it corrects the mutual influence between signals from different sources (such as GPR and remote sensing imagery) at the data fusion level and verifies physical consistency. This ensures the physical and logical integrity and reliability of the initial parameter field of the model and solves the problem of lack of physical coupling in multi-source data fusion.

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

Claims

1. A complex slope runoff simulation method based on multi-source data fusion, characterized in that, The method comprises: acquiring multi-source spatial data to generate an initial land surface parameter field containing a digital elevation model and soil parameter distribution; based on the initial land surface parameter field and rainfall data, a dynamic overflow network and infiltration model are constructed to calculate dynamic overflow data and dynamic infiltration data; a feedback regulation mechanism between overflow, erosion and infiltration is established according to the dynamic overflow data and the initial land surface parameter field, and the slope runoff result is calculated in combination with the dynamic infiltration data; analyzing the slope runoff result to output runoff early warning information; wherein the feedback regulation mechanism between overflow, erosion and infiltration comprises: analyzing the dynamic overflow data, and combining the initial land surface parameter field to evaluate a scouring intensity distribution field representing the erosion ability of water flow; based on the scouring intensity distribution field, the soil water conductivity in the initial land surface parameter field is used as an initial value to generate a corrected soil water conductivity parameter field; the corrected soil water conductivity parameter field is applied to subsequent infiltration and confluence calculation; wherein: The assessment yielded the scour intensity distribution field, including: determining the bed shear stress based on the flow velocity and duration in the dynamic overflow data and in conjunction with the initial surface parameter field; quantifying the scour intensity distribution field, specifically including: calculating the actual bed shear stress τ within the channel, and combining it with the critical soil shear stress τ obtained from the initial surface parameter field. _c (Pa), and the average velocity v of the preferred flow path obtained from the dynamic overflow data. _path (m / s) and duration t _duration (s) can be used to calculate the scour intensity index E for each path. _path ; E of each path _path The values ​​are assigned to the grid through which the scour flow passes, forming the scour intensity distribution field E(x,y); where E _path =(τ-τ _c )×v _path 2 ×t _duration ;(τ-τ _c The term represents the net water flow force exceeding the soil's erosion resistance, v _path 2 The term reflects the kinetic energy of the water flow, t _duration The term reflects the cumulative effect of erosion; τ = ρ × g × R × S _f ρ is the density of water, g is the acceleration due to gravity, R is the hydraulic radius, and S is the velocity of water. _f For friction slope; generating the corrected soil water conductivity parameter field comprises: determining a fine particle loss rate according to the scouring intensity distribution field, and fusing the scouring intensity distribution field and the fine particle loss rate to dynamically adjust the soil water conductivity in the initial land surface parameter field to obtain the corrected soil water conductivity parameter field; dynamically adjusting the soil water conductivity in the initial land surface parameter field comprises: According to the scouring intensity distribution field, a hydraulic conductivity reduction coefficient (1-α×E _path / E _duration ) is defined to represent the influence of scouring on soil compaction, where α is a weight coefficient of scouring compaction, and E _duration is a critical scouring intensity threshold value; Based on the fine particle loss rate, the increase in hydraulic conductivity (1+γ×Fine) is defined as the factor characterizing the effect of fine particle loss on pore connectivity. _loss ), where γ is the weighting coefficient for the increase in porosity, Fine _loss For fine particle loss rate; The initial soil saturated hydraulic conductivity is multiplied by the hydraulic conductivity reduction factor and the hydraulic conductivity increase factor to obtain the modified soil saturated hydraulic conductivity Ks _modified = Ks _initial × (1 - 0.2 x E _path / E _duration ) x (1 + 1.5 x Fine _loss ), Ks _initial is the initial soil saturated hydraulic conductivity.

2. The method of claim 1, wherein, applying the corrected soil water conductivity parameter field to subsequent infiltration and confluence calculation is performed by a multi-time scale separation algorithm, which comprises: distinguishing and obtaining the characteristic time scales corresponding to the overflow process, the infiltration process and the formation process of the scouring intensity distribution field; determining an adaptive time step according to the obtained multiple characteristic time scales to perform subsequent infiltration and confluence calculation within the time step.

3. The method of claim 2, wherein, performing subsequent infiltration and confluence calculation within the adaptive time step follows a predetermined execution order, which is: performing calculation corresponding to the overflow process; performing calculation corresponding to the infiltration process; performing calculation corresponding to the formation and correction of the scouring intensity distribution field and the corrected soil water conductivity parameter field.

4. The method of claim 1, wherein, generating the initial land surface parameter field comprises: fusing unmanned aerial vehicle and ground laser point cloud data to generate a digital elevation model; processing ground penetrating radar echo signal data to obtain soil parameter distribution by inversion; analyzing multi-spectral remote sensing images to extract vegetation cover information as part of the initial land surface parameter field.

5. The method of claim 1, wherein, Generating the initial land surface parameter field can also include: extracting vegetation feature information from multi-source spatial data; building a physical coupling relationship between vegetation feature information and different source data signals in multi-source spatial data; based on the physical coupling relationship, mutually correcting and fusing different source data signals to generate the initial land surface parameter field; wherein: building a physical coupling relationship comprises: based on vegetation feature information, building a three-dimensional root-soil transfer function representing the influence of root system on soil electromagnetic properties; mutually correcting and fusing different source data signals comprises: applying the three-dimensional root-soil transfer function to root attenuation correction of ground penetrating radar echo signal data in multi-source spatial data, and obtaining soil parameter distribution based on the corrected ground penetrating radar echo signal data.

6. The method of claim 5, wherein, The method further comprises physically consistent verification of the fusion result, comprising: Extracting a topographic index and a soil moisture distribution index from the initial land surface parameter field; Calculating a correlation coefficient between the topographic index and the soil moisture distribution index to quantify the physical consistency between the land surface topography and the soil moisture distribution; Adaptively determining a final value of the initial land surface parameter field based on the correlation coefficient.

Images