A dynamic DEM spatial interpolation method

Abstract

The application discloses a dynamic DEM spatial interpolation method, fills the DEM data for depression, flow direction analysis, flow accumulation calculation, and adopts a minimum catchment area threshold method to extract a natural sub-basin unit; a terrain feature matrix is constructed, a first-order basin is finely classified through an improved self-organizing mapping network, a hydrological response unit is generated through boundary processing, and a hydrological attribute database is established; multi-source data is fused, attribute interpolation of a lower surface and rainfall is supplemented, an interpolation auxiliary parameter system is constructed; a basin is divided into regular grids as neurons, a dynamic neural network containing dynamic states and static attributes is constructed, and nonlinear mapping of DEM correction parameters is realized through dynamic activation functions and loss functions optimization. The method overcomes defects of traditional interpolation algorithms, such as not considering hydrological boundary constraints and fixed topological structure of a neural network model, and improves DEM interpolation precision and hydrological simulation effect in a complex terrain area.

Description

Technical Field

[0001] This invention relates to the field of interdisciplinary technology of geographic information systems and hydrological models, specifically to a dynamic DEM spatial interpolation method, which is applicable to high-precision hydrological simulation, geological disaster early warning and water resource management. Background Technology

[0002] Digital elevation models (DEMs), as fundamental data for depicting the topographic relief of the Earth's surface, are core support for fields such as hydrological cycle simulation, watershed water resource management, geological disaster early warning, and ecological pattern analysis. Their accuracy directly determines the reliability of key topographic parameters such as river network extraction, slope and aspect calculation, and confluence path simulation, thus affecting the scientific validity of flood control engineering design, water resource allocation schemes, and ecological protection planning. For example, in flood forecasting of small and medium-sized watersheds, a DEM interpolation error of only 0.5 meters can lead to a confluence time deviation of more than 15%, significantly increasing the risk of inaccurate flood peak warnings. In permafrost ecological research, distortion in the depiction of topographic relief directly affects the simulation of permafrost thawing rates, weakening the accuracy of hydrological process assessment in cold regions.

[0003] The limitations of traditional spatial interpolation methods are particularly pronounced in complex geographical environments. Taking Kriging as an example, its assumption of optimal unbiased estimation based on spatial autocorrelation is difficult to apply in fragmented karst landforms or plateau canyon areas—these regions exhibit dramatic spatial heterogeneity due to abrupt topographic changes caused by lithological differences and fault cutting. Kriging's semi-variogram model is prone to oversmoothing, obscuring key topographic features such as steep cliffs and gullies. The inverse distance weighting method, influenced by the subjective setting of the distance attenuation coefficient, often results in an imbalanced weight distribution in data-sparse mountainous areas (such as the marginal watersheds of the Qinghai-Tibet Plateau), leading to the "isolated point domination" phenomenon. This causes the elevation values ​​of adjacent sub-basins to be incorrectly flattened, disrupting the hierarchical structure of the natural drainage system.

[0004] While existing neural network models demonstrate advantages in nonlinear fitting, they suffer from a disconnect between data and physics. Most models use randomly sampled discrete elevation points as training samples, failing to incorporate watershed hydrological topology (such as the confluence levels of tributaries and main streams, and the hydrological isolation effect of watersheds), resulting in spatial patterns lacking physical plausibility. For example, in rainfall-induced soil erosion zones, existing models cannot dynamically adjust topographic weights using real-time rainfall data, making it difficult to capture the effects of raindrop splash erosion and slope runoff on micro-topography. This leaves the updated DEM at a static topographic representation level, with insufficient coupling to actual hydrological processes.

[0005] At the spatial partitioning level, conventional regular grid partitioning (e.g., 1km×1km) in the Hengduan Mountains with its dramatic topographic relief can divide a complete slope of the same aspect into multiple grids, leading to fragmented calculations of slope length and gradient. Furthermore, partitioning based on administrative boundaries (e.g., county-level units) may split the entire watershed into different management units, causing "cross-boundary breaks" in the simulation of confluence processes and failing to reflect the natural continuity of hydrological processes. This partitioning method contradicts the principle that Hydrological Response Units (HRUs) must adhere to the synergy of "topography-soil-vegetation," directly increasing the uncertainty of parameters in distributed hydrological models.

[0006] In recent years, although some studies have attempted to introduce hydrological constraints (such as interpolation results based on river network corrections), these have mostly been limited to single-scale watersheds, failing to form a multi-level nested mechanism of "whole watershed - sub-watershed - micro-watershed," making it difficult to simultaneously consider macroscopic confluence trends and microscopic topographic details. Furthermore, some dynamic interpolation models rely excessively on high-resolution observation data, limiting their applicability in underdeveloped areas where data is scarce. Therefore, constructing a DEM spatial interpolation method that combines physical constraints, dynamic response, and multi-scale adaptability has become a key issue in overcoming the accuracy bottleneck of hydrological simulation in complex topographic regions. Summary of the Invention

[0007] The purpose of this invention is to provide a dynamic DEM spatial interpolation method that addresses the shortcomings of existing DEM spatial interpolation methods in areas such as watershed classification processing, attribute supplementation, and neural network construction, and provides a dynamic DEM spatial interpolation method with higher accuracy and better dynamic response.

[0008] To achieve the above functions, this invention designs a dynamic DEM spatial interpolation method, including the following steps S1-S4, to complete the spatial interpolation of DEM data and the hydrological simulation of the target watershed:

[0009] Step S1: For the target watershed, collect raw DEM data, perform primary segmentation of the raw DEM data based on ArcGIS hydrological analysis toolkit, and extract natural sub-watershed units through depression filling, flow direction analysis, runoff accumulation calculation and minimum catchment area threshold method.

[0010] Step S2: For natural sub-basin units, construct a topographic feature matrix, use an improved self-organizing mapping network for secondary classification, optimize clustering through dynamic grid mechanism and flow direction constraint, and generate hydrological response units through Voronoi segmentation with hydrological constraints and boundary morphology optimization. Establish a hydrological attribute database for hydrological response units.

[0011] Step S3: Integrate multi-source data, including meteorological elements, geological data, land use type data, and vegetation parameters; supplement attribute interpolation for each multi-source data to construct an interpolation auxiliary parameter system;

[0012] Step S4: Construct a dynamic neural network, divide the target watershed into a neuron grid, build topological connections based on the D8 flow direction, simulate the hydrological response through a dynamic activation function, and train and optimize the parameters of the dynamic neural network using the root mean square error (RMSE) as the loss function.

[0013] Beneficial effects: Compared with the prior art, the advantages of the present invention include:

[0014] This invention designs a dynamic DEM spatial interpolation method. This method is not limited to a single-scale watershed, but forms a multi-level nested mechanism of "whole watershed - sub-watershed - micro-watershed", taking into account both macro-level confluence trends and micro-level topographic details. It constructs a dynamic DEM spatial interpolation method that combines physical constraints, dynamic response, and multi-scale adaptability. This method overcomes the shortcomings of traditional interpolation algorithms, such as not considering hydrological boundary constraints and fixed topology of neural network models. It improves the DEM interpolation accuracy and hydrological simulation effect in complex terrain areas, and is applicable to the fields of high-precision hydrological simulation, geological disaster early warning, and water resource management. Attached Figure Description

[0015] Figure 1 This is a flowchart of a dynamic DEM spatial interpolation method provided by an embodiment of the present invention;

[0016] Figure 2 This is a schematic diagram of the watershed DEM of the study area provided according to an embodiment of the present invention;

[0017] Figure 3 This is a schematic diagram of a primary watershed and hydrological response unit provided according to an embodiment of the present invention;

[0018] Figure 4 This is a conceptual diagram of a dynamic neural network provided according to an embodiment of the present invention. Detailed Implementation

[0019] The present invention will be further described below with reference to the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and should not be used to limit the scope of protection of the present invention.

[0020] This invention provides a dynamic DEM spatial interpolation method, referring to... Figure 1 The process includes the following steps S1-S4, which complete the spatial interpolation of the DEM data and the hydrological simulation of the target watershed:

[0021] Step S1: For the target watershed, collect raw DEM (Digital Elevation Model) data with a resolution of 30m. Based on the ArcGIS hydrological analysis toolkit, perform first-level segmentation of the raw DEM data, and extract natural sub-watershed units through depression filling, flow direction analysis, runoff accumulation calculation and minimum catchment area threshold method.

[0022] The specific steps of step S1 are as follows:

[0023] Step S1.1: Filling to eliminate depressions in the original DEM data;

[0024] Step S1.2: Perform flow direction analysis based on the D8 algorithm to generate the D8 flow direction matrix;

[0025] Step S1.3: Flow Accumulation generates the runoff path network;

[0026] Step S1.4: Based on the minimum catchment area threshold method (threshold range: 1-5km) 2 Automatically extract natural sub-basin unit boundaries;

[0027] Step S1.5: Output the vectorized set of natural sub-basin units U = {u1, u2, ..., u} n}, where u1, u2, ..., u n This represents each natural sub-basin unit, where n is the total number of natural sub-basin units.

[0028] Step S2: For natural sub-basin units, construct a topographic feature matrix, use an improved self-organizing map network (SOM) for secondary classification, optimize clustering through dynamic grid mechanism and flow direction constraint, and generate hydrological response units through Voronoi segmentation with hydrological constraints and boundary morphology optimization. Establish a hydrological attribute database for hydrological response units.

[0029] The specific steps of step S2 are as follows:

[0030] Step S2.1: Construct a terrain feature matrix for the natural sub-basin units extracted in Step S1:

[0031]

[0032] In the formula, T i represents the topographic feature matrix of the i-th natural sub-basin unit, and m represents the number of units in the study area where topographic features are calculated or sampled.

[0033] The elevation uses the original elevation, and the slope aspect uses azimuth quantification improvement.

[0034] The slope algorithm uses a third-order inverse distance squared weighting algorithm:

[0035]

[0036] Where x, y, and Z represent the east-west direction, north-south direction, and elevation of the geographic coordinate system, respectively; This represents the rate of change of elevation Z in the east-west direction; This represents the rate of change of elevation Z in the north-south direction;

[0037] The curvature is calculated as follows:

[0038]

[0039] The Terrain Humidity Index (TWI) is calculated using the following formula:

[0040]

[0041] Among them, A s β represents the upstream catchment area, and β is the slope angle.

[0042] Step S2.2: An improved self-organizing map (SOM) clustering algorithm is used to dynamically mesh natural sub-basin units. The mesh size is adaptively adjusted according to the terrain complexity, and the elevation standard deviation σ is calculated. Z If σ Z If σ < 10, the natural sub-basin unit is a flat area, the grid size is reduced, and the local terrain resolution is improved; if 10 ≤ σ Z If the value is ≤30, then the natural sub-basin unit is a medium-complexity region, a transitional region with terrain complexity between flat and steep areas. The algorithm employs a smooth adaptive strategy; the grid size is based on σ. Z The value of σ is linearly interpolated between the minimum and maximum values; if σ Z If the value is greater than 30, then the natural sub-basin unit is a steep region, and the grid size is increased to avoid over-segmentation, as shown in the following formula:

[0043]

[0044] In the formula, GRIDsize represents the dynamic grid size, and N1 is the number of grid cells in the natural sub-basin cell;

[0045] The dynamic mesh generation mechanism and the flow direction constraint function work together to improve distance calculation by introducing the flow direction constraint function:

[0046]

[0047] In the formula, D′ represents the flow direction constraint function, ω k X k C k The meaning, parameter definition, value range, and hydrological significance of these parameters are summarized in Table 1 below:

[0048] Table 1. ω in the flow direction constraint function k X k C k Meaning, parameter definition, value range, and hydrological significance

[0049]

[0050] In the flow direction constraint function, λ and θ ij φ flow The definitions, meanings, calculation methods, and hydrological significance of the parameters are summarized in Table 2 below:

[0051] Table 2. λ, θ in the flow direction constraint function ij φ flow Parameter definition, meaning, calculation method, and hydrological significance

[0052]

[0053] The flow weight λ is calculated as follows:

[0054] λ=0.4×(1-e -0.1×FlowAcc )

[0055] In the formula, FlowAcc represents the cumulative flow amount;

[0056] River channel area (FlowAcc>100): λ≈0.4 (strong constraint). Slope area (FlowAcc<5): λ<0.1 (weak constraint). The topographic feature distance term calculates the difference between the current grid and the cluster center in five topographic features: elevation, slope, aspect, curvature, and topographic moisture index (TWI) using weighted Euclidean distance. The feature weight coefficients are shown in Table 1. This weighting configuration highlights the dominant role of elevation zonation and soil moisture in hydrological processes. The flow direction constraint introduces a hydrophysical mechanism, adjusting the consistency of flow direction through a dynamic weighting formula. When the cumulative flow (FlowAcc) increases (e.g., in river channel areas), λ approaches 0.4 to force the flow direction to align with the mainstream direction determined by the D8 algorithm. In the slope area (FlowAcc<5), the λ value is below 0.1 to weaken the flow direction constraint. Angle difference calculation is used to eliminate circumferential discontinuities, ensuring equivalence between 0° and 360° directions. By coupling topographic features with hydrological flow direction, contour hydrological response units (HRUs) are formed in mountainous areas, linear hydrological response units are generated in river channels, and soil moisture-dominated hydrological response units are constructed in plains, significantly improving the consistency of water flow direction within hydrological response units.

[0057] Step S2.3: Output the clustering result set P from the K winning neurons of the output layer of the improved self-organizing map network (SOM):

[0058]

[0059] In the formula, P k E represents the k-th clustering result. k S is the average elevation. kFor the overall slope, A k As the dominant slope aspect, C k For characteristic curvature, T k The mean of the topographic humidity index; each neuron corresponds to a topographic homogeneous unit, and its position in the feature space determines the generation base point of the Voronoi diagram;

[0060] Step S2.4: Perform hydrologically constrained Voronoi segmentation in the Geographic Information System (GIS), and transform the clustering results into a set of hydrological response units {HRU1, HRU2, ..., HRU} using the Voronoi diagram. K The hydrological response unit is as follows:

[0061]

[0062] In the formula, d hydro Let P represent the hydrological distance function, and let p represent the representative point of the clustering result. j This represents the j-th clustering result;

[0063] The hydrological distance function is as follows:

[0064] d hydro =α·d Euclidean +(1+α)·d flow

[0065] In the formula, α is a dynamic weighting coefficient between 0 and 1, i.e., α∈[0,1], which is used to adjust the relative contribution of the Euclidean distance of two-dimensional geographic coordinates and the cumulative distance along the D8 flow path in the hydrological distance calculation; d Euclidean d represents the Euclidean distance between two-dimensional geographic coordinates. flow This represents the cumulative distance along the D8 water flow path;

[0066] Step S2.5: Extract the river network Γ generated by the D8 algorithm channel As a hard constraint, specifically To constrain the re-boundary delineation of hydrological response units, the channel attractor algorithm is adopted;

[0067] Step S2.6: For areas smaller than threshold A min The hydrological response units are fused together, as shown in the following formula:

[0068]

[0069] In the formula, HRU small This indicates that the area is less than the threshold A. min Hydrological Response Unit, HRU j Indicates candidate neighboring hydrological response units (i.e., HRUs) small (Potential target HRUs to be merged); in one embodiment, A is setmin =0.01km 2 ΔH is the elevation difference, and ΔL is the flow path length. It is preferentially integrated into the downstream hydrological response unit (HRU) with continuous topography.

[0070] Step S2.7: Boundary morphology optimization is the core step in generating hydrological response units. Physical rationalization of the hydrological boundaries is achieved through two-stage processing. Boundary morphology optimization is performed on the merged hydrological response units, using morphological opening operations to eliminate microscopic jagged edges, as shown in the following formula:

[0071]

[0072] In the formula, B raw The raw hydrological boundary refers to the initial boundary polygon output directly from the hydrological model, which exhibits pixel-level jaggedness due to DEM discretization (typical resolution ≤30 meters); B final The term "morphologically optimized boundary" refers to the intermediate boundary after the expansion-corrosion cascade treatment; SE 3×3 The 3×3 structuring element is defined as a 90m×90m (3 times the DEM pixel scale) square binary template with the origin at the center, used to uniformly process the target pixel and its eight neighboring relationships.

[0073] The above process uses a 3×3 structuring element to dilate the original boundary, filling small fractures and gaps caused by DEM errors, followed by an equal amount of erosion to remove isolated protrusions. This cascaded dilation-erosion process effectively eliminates topographic noise at the pixel scale (typically <90 meters), reducing the boundary length by 12%-18%, while maintaining the topological connectivity between the river channel and the slope.

[0074] For the morphologically optimized hydrological boundary, a Bézier curve fitting based on hydrological constraints is used, as shown in the following formula:

[0075] B smooth =BezierCurveFit(B final tolerance = 0.5δ)

[0076] In the formula, B smooth The term "Hydrologically Smoothed Boundary" refers to the final hydrological boundary optimized by the Bézier curve with flow direction constraints, satisfying the consistency between sub-pixel accuracy and runoff direction. "Tolerance" is the tolerance threshold, and "δ" represents the resolution of the DEM data.

[0077] Control points were set every 60 meters (2 times DEM resolution) along the optimization boundary, and the spacing was increased to 15 meters at points of abrupt curvature change. A cubic Bézier curve parametric reconstruction was employed, with its core innovation being the projection of control points along the water flow direction: by correlating with D8 flow direction data, the angle between the curve tangent and the surface runoff direction was ensured to be ≤5°. A tolerance threshold of 0.5 times DEM resolution (15 meters) was specifically set (i.e., 0.5δ in the formula). When the fitting deviation exceeded this limit, new control points were automatically inserted, maintaining boundary accuracy at the sub-pixel scale.

[0078] Step S2.8: For the hydrological response units obtained in step S2.7, establish a hydrological attribute database as shown in Table 3 below:

[0079] Table 3. Hydrological Attribute Database of Hydrological Response Units

[0080]

[0081] Step S3: Integrate multi-source data, including meteorological elements, geological data, land use type data, and vegetation parameters; supplement attribute interpolation for each multi-source data to construct an interpolation auxiliary parameter system;

[0082] The specific steps of step S3 are as follows:

[0083] Step S3.1: Regarding meteorological elements, raster data (temporal resolution ≤ 1 hour) of precipitation, temperature, and wind speed are generated using interpolation based on the ANUSPLIN algorithm. The input data consists of hourly observations from national meteorological stations within a 30km radius of the target watershed, with DEM data introduced as a covariate to eliminate topographic shading effects. Data quality control is performed before interpolation: outliers exceeding 3 standard deviations are removed, and linear interpolation is used to fill in observation gaps of ≤ 3 hours. The spline function order (3rd order recommended) and tension parameters (range 0.01-0.5) are optimized using cross-validation to ensure the mean absolute error (MAE) of the interpolation results is: precipitation ≤ 5%, temperature ≤ 0.5℃, and wind speed ≤ 0.3m / s. The final output is a time-series dataset of meteorological elements (WGS84 coordinate system, UTM projection) with the same resolution as the target DEM data.

[0084] Step S3.2: Regarding geological data and soil permeability coefficient interpolation, the data substitution adopts the permeability coefficient reference values ​​from the soil survey results and the 1:100,000 soil texture map (10-50 m / d for sandy soil, 1-10 m / d for loamy soil, and 0.1-1 m / d for clay soil). The interpolation optimization uses the co-kriging method, with soil texture as an auxiliary variable. The exponential model is used for sandy soil (spatial correlation range of 2-5 km), and the Gaussian model is used for clay soil (1-3 km). The nugget value is set to 10%-15% of the theoretical value, and the error compared with the measured value at the irrigation station is ≤15%. When interpolating lithological types, data integration is based on extracting lithological boundaries from 1:25,000 geological maps. A lithological coding system of 10 primary categories and 30 secondary categories is established in conjunction with geological reports. Geophysical data is introduced to enhance the identification of concealed lithologies. The method is adjusted to use probabilistic kriging, extending the lithological contact zone width to 100-200m, with a buffer weight of 1.2-1.5 times. Verification with borehole data shows a type matching accuracy of ≥80%. For quality control, data from authoritative institutions over the past 15 years is prioritized, with 20% of samples cross-validated (supplementary corrections when the error is >30%). Output data is labeled with a reliability level (A / B / C) to ensure compliance with DEM data interpolation requirements (resolution matching ≥90%, type consistency ≥85%).

[0085] Step S3.3: Regarding land use type data, land use type data is obtained by performing random forest classification based on Sentinel-2 satellite imagery (10m resolution). Image selection must meet the requirement of cloud cover <10% during the growing season (June-August). Preprocessing includes atmospheric correction (Sen2Cor plugin), topographic correction (C-correction method), and radiometric normalization. The feature variable set includes 10 band reflectance (B2-B8A, B11-B12) and 4 vegetation indices (EVI, SAVI, MSAVI, RVI). Stratified sampling is used to divide the data into training samples (70%) and validation samples (30%). The classifier parameters are set as follows: 500 decision trees, minimum sample size for node splitting (5), and out-of-bag (OOB) error controlled within 8%.

[0086] Step S3.4: Regarding vegetation parameters, the NDVI index (Normalized Difference Vegetation Index) is calculated as the core vegetation parameter. Near-infrared band (B8, 842nm) and red band (B4, 665nm) of Sentinel-2 imagery are used, calculated using the formula NDVI = (B8 - B4) / (B8 + B4), with a value range of [-1, 1]. To eliminate the influence of atmospheric aerosols, the original bands are first corrected using BRDF (Bidirectional Reflectance Distribution Function), and then noise is smoothed using median filtering (3×3 window). For seasonal vegetation changes, a monthly maximum NDVI composite (MVC) product needs to be constructed. Cloud-polluted pixels are removed using time-series harmonic analysis (HANTS). The final NDVI dataset needs to be correlated with the concurrent field vegetation survey data (leaf area index LAI) for verification (R). 2 ≥0.75).

[0087] Step S4: Construct a dynamic neural network, divide the target watershed into a neuron grid, build topological connections based on the D8 flow direction, simulate the hydrological response through a dynamic activation function, and train and optimize the parameters of the dynamic neural network using the root mean square error (RMSE) as the loss function.

[0088] The specific steps of step S4 are as follows:

[0089] Step S4.1: Divide the target watershed into regular grids, such as 100×100m. The grid size can be adaptively adjusted according to the original DEM resolution (30m or 10m recommended), ensuring that the topographic slope variation coefficient within a single grid is ≤5%. Treat each grid cell as a neuron. Employ a spatial matching algorithm between hydrological response cells and grids. When a hydrological response cell spans multiple grids, the center point of the grid containing the centroid of the hydrological response cell is used as the reference position for the neuron. Edge grids are buffered (width = 1 / 2 of the grid side length) to ensure complete coverage of the watershed boundary. Define the neuron coordinates as follows:

[0090] X i =(x i ,y i )∈R 2 i = 1, 2, ..., N2

[0091] In the formula, R represents the set of real numbers, N2 is the total number of grid cells after the entire target watershed is divided into regular grids (N2 ≥ 5000 must be satisfied to ensure detailed watershed characterization), X i Let x be the coordinates of the i-th neuron. i ,y i) represents the coordinates of the i-th neuron in the X and Y directions in the unified WGS84-UTM projection coordinate system in step S3. The coordinates of the grid center points of the DEM data are calculated in batches, and missing areas (such as water areas) are filled by interpolation using the inverse distance weighting method (IDW).

[0092] Step S4.2: Each neuron contains dynamic states and static properties, defined as follows:

[0093]

[0094] In the formula, θ i H represents the i-th neuron. i , k i The dynamic state of the i-th neuron is represented by water depth, hydraulic gradient, and conductivity, respectively; A i The static attribute matrix of the i-th neuron contains the multi-source data obtained in step S3, including meteorological elements, geological data, land use type data, and vegetation parameters; static attribute matrix A i It includes 8-dimensional features: the multi-year average values ​​of meteorological elements in step S3.1 (rainfall, temperature, wind speed), the geological data in step S3.2 (specifically, soil texture coding), the land use type index in step S3.3, and the mean NDVI value, slope, and aspect in step S3.4 (calculated from DEM). All attributes need to be min-max normalized to the [0,1] interval to eliminate the influence of dimensions.

[0095] Step S4.3: Determine the connection topology based on the D8 water flow direction (downstream neurons receive upstream input). By calculating the slope of each grid and its 8 neighboring grids, the water flow direction is directed towards the downstream grid with the steepest slope. For flat areas with a slope <0.5°, the D-infinity algorithm is used to assist in determining the flow direction, avoiding topological breaks caused by ambiguity in the flow direction. The neuron at the watershed outlet is set as the global sink, receiving input from all downstream terminal neurons.

[0096] Step S4.4: Construct the dynamic activation function as follows:

[0097]

[0098] In the formula, f act This represents the dynamic activation function, where β is the hydrological conduction coefficient (determined by soil type and slope), and Rt is the real-time rainfall intensity (input via the meteorological interpolation module). This indicates that the activation intensity varies nonlinearly with both runoff velocity and rainfall intensity, simulating a real hydrological response. The physical meaning of the activation function is: when f... act When f approaches 1, the neuron is in a strong confluence state (the grid's current generation capacity is significantly enhanced); when f...act When the value approaches -1, it is in a weak confluence state (mainly water storage). The threshold response characteristics of rainfall-runoff are simulated through nonlinear transformation.

[0099] Step S4.5: Construct the loss function as follows:

[0100]

[0101] In the formula, Loss represents the loss function, and N3 is the total number of time steps in the hydrological simulation; This represents the actual flow rate observed at time point i, in meters per second (m³). 3 / s; This represents the predicted (simulated) flow rate at time point i, in meters per second (m³). 3 / s.

[0102] The following is an application example of the present invention:

[0103] The Tunxi River Basin was chosen as the study area. As an important part of the Xin'an River Basin, it is located in the southeastern part of the mountainous region of southern Anhui. Situated in the main area of ​​the Huangshan Mountains, it borders the Yangtze River system to the west and north, and is adjacent to the Tianmu and Baiji Mountains to the southeast. The basin contains several sizable intermontane basins and valleys, exhibiting a topographical characteristic of being higher in the west and lower in the east, with significant topographic relief and generally steep river slopes. Climatically, it belongs to the typical subtropical monsoon climate zone, with consistently high annual precipitation and an average annual temperature in the transitional zone between temperate and subtropical. The basin has good vegetation cover, mainly consisting of various forest types and agricultural land, forming a complex ecosystem. Elevation varies widely, showing a clear decreasing trend from the western mountainous areas to the eastern river valleys, creating a significant topographic drop.

[0104] The implementation process begins with the collection of basic data, encompassing multi-dimensional information including hydrology, meteorology, topography, and socioeconomic data. The core dataset includes real-time hydrological monitoring records, soil physical property parameters, precipitation observation data at different time scales, high-precision topographic elevation models, population distribution statistics, and historical flood characteristic values. This data is categorized into two sets based on processing requirements: a routine processing set and a set to be processed, corresponding to stable hydrological elements and dynamically changing elements, respectively. Data collection included: DEM data: 30m resolution SRTM DEM, covering the watershed and a 5km buffer zone (matrix D∈R6150×5200); meteorological data: hourly rainfall (P), temperature (T), and wind speed (W) from 12 national meteorological stations over the past 6 years; hydrological data: hourly flow rate (Q) from Tunxi station; geological data: soil texture map (1:100,000): sandy soil 32%, loamy soil 45%, clay soil 23%, and measured permeability coefficient (Ks) from 25 boreholes; remote sensing data: Sentinel-2 imagery (June-August 2020-2022, cloud cover <10%), and 62 measured LAI points (35 in forest, 18 in farmland, and 9 in grassland). A schematic diagram of the watershed DEM for the study area is provided. Figure 2 .

[0105] Subsequently, data standardization preprocessing was performed to convert all indicators to a unified dimensional system. Based on this, a combined weighting method was used to determine the weight coefficients of each influencing factor: first, a judgment matrix was constructed using an expert decision-making model to calculate the subjective weight distribution; simultaneously, based on information entropy theory, the discrete characteristics of the data itself were analyzed to derive an objective weight allocation scheme. The two weighting results were then weighted and fused to form a comprehensive weighting system, providing a quantitative basis for subsequent analysis. This dual-weighting mechanism incorporates domain knowledge and experience while fully respecting the objective laws of the data, ensuring the scientific validity and reliability of the evaluation system.

[0106] In the spatial discretization stage, the grid partitioning scheme is dynamically adjusted according to the terrain complexity, while a flow direction constraint mechanism is introduced to optimize the boundaries of hydrological response units. Morphological operations and curve fitting techniques are used to smooth the generated boundaries, significantly improving the spatial continuity of the hydrological network. For extreme events such as torrential rains, a dynamic activation function is used to adjust parameter sensitivity in real time, enabling the model to adapt to changes in different hydrological conditions. The entire process ultimately constructs a multi-level nested hydrological discretization system, forming a well-defined three-dimensional structural framework from sub-basin partitioning to hydrological response unit definition and refined grid calculation.

[0107] The specific calculations are as follows:

[0108] Step S1: ArcGIS processes the watersheds to obtain the primary classification, and then uses the ArcGIS Hydrological Analysis Toolkit for further processing.

[0109] 1. Depression filling treatment: Repair M depressions in the DEM;

[0110] 2. Flow direction analysis: The D8 algorithm generates the water flow direction matrix;

[0111] 3. Sub-basin division: based on the cumulative runoff threshold F grid;

[0112] U = {U1, U2, ..., U14}

[0113] Output: N natural sub-basins (areas ranging from 82 to 168 km²) 2 ).

[0114] Step S2:

[0115] S2.1: Taking a typical sub-basin U i For example, a spatial matrix containing 5-dimensional terrain features is constructed. First, the slope S is calculated using a third-order inverse distance squared weighting algorithm to accurately characterize the steepness of the land surface; then, the curvature is solved using the second derivative to describe the terrain's unevenness; finally, the terrain humidity index TWI is derived based on the upstream catchment area As and the slope angle. The final terrain feature matrix is ​​then formed.

[0116]

[0117] The slope is calculated using a third-order inverse distance squared weighting algorithm:

[0118]

[0119] Curvature:

[0120]

[0121] Terrain Humidity Index (TWI):

[0122]

[0123] Where A s Let β represent the upstream catchment area and β represent the slope angle. This matrix integrates five types of terrain attributes from F grids, providing a data foundation for cluster analysis.

[0124] S2.2: Integrating dynamic grid mechanism with hydrological flow direction constraints:

[0125] Grid Adaptive: Based on the standard deviation of elevation σ Z =a dynamically adjusts the clustering resolution, refines the mesh in flat regions (a<10), and expands the mesh in steep regions (a>30), and calculates the dynamic mesh size:

[0126] (Original 30m → Adjusted bm)

[0127] Flow direction constraint: Constructing a mixed distance function:

[0128]

[0129] Where the flow direction weight λ = 0.4 × (1 - e -0.1×F Dynamic changes with the cumulative flow:

[0130] Slope area (F<5): λ slope =c = 0.4(1 - e - 0.1 × 3) ≈d

[0131] River channel area (F>100): λ river =e = 0.4(1 - e - 0.1 × 120) ≈ f

[0132] Where c is the expression for the flow direction weight in the slope area; e is the expression for the flow direction weight in the river channel area; d is the numerical value of the flow direction weight in the slope area; and f is the numerical value of the flow direction weight in the river channel area.

[0133] This design ensures that contour-zoned HRUs are generated in mountainous areas and linear HRUs are formed in riverine areas.

[0134] S2.3: Terrain Prototype Generation

[0135] Clustering outputs K terrain prototypes P k =(E k ,S k A k C k ,T k For example, P1 = (e1,s1,a1,c1,t1) represents the average elevation e1, comprehensive slope s1, dominant aspect a1, characteristic curvature c1, and TWI mean t1 of a certain type of homogeneous terrain unit. These prototypes constitute the benchmark points for HRU division.

[0136] Output K terrain prototypes, for example:

[0137] P5 = (g, h, i, j, k) = (1, 024 m, 31.2°, Southeast, 0.045, 9.8)

[0138] P12 = (l, m, n, o, p) = (632m, 8.7°, northeast, -0.018, 12.3)

[0139] S2.4~S2.7: HRU space generation:

[0140] 1. Hydrologically Constrained Voronoi Segmentation: Define the hydrological distance function d hydro =α·d Euclidean +(1-α)·d flow Comprehensive geographical distance d Euclidean Distance d from the water flow pathflow ;

[0141] 2. Topology optimization: Based on the river network Γ channel For hard constraints, the river attractor algorithm is used to reconstruct the boundary;

[0142] 3. Micro-unit fusion: Perform downstream fusion on HRUs with an area of ​​less than 0.01 km². Prioritize incorporating downstream units with continuous topography. Integrate D units with an elevation range of <0.01km. 2 HRU;

[0143] 4. Boundary optimization:

[0144] Morphological treatment: Eliminate jagged edges qkm;

[0145] Bézier curve fitting: B smooth =BezierFit(B final (0.5δ) Maintaining consistency in water flow direction (tangential angle ≤ 5°) The final output is a set of hydrological response units {HRU1, HRU2, ..., HRU}. K}; Schematic diagram of first-level watershed and hydrological response unit (refer to...) Figure 3 .

[0146] S2.8: The hydrological attribute database using HRU_5 as an example is shown in Table 4 below:

[0147] Table 4. Hydrological Attribute Database (HRU_5 Example)

[0148]

[0149] Step S3:

[0150] 1. ANUSPLIN interpolation parameters: 3rd order spline, tension coefficient 0.25;

[0151] Accuracy verification:

[0152] 2. Geological interpolation:

[0153] Spatial distribution of soil permeability coefficient:

[0154] Lithological classification accuracy: z% (>80%);

[0155] 3. Land use classification: Random forest classification: 500 decision trees; OOB error = AN%;

[0156] Main types: Forest land AO%, Farmland AP%, Water area AQ%, Construction land AR%;

[0157] 4. Vegetation parameters:

[0158] LAI correlation: R 2 =AS(>0.75);

[0159] Step S4:

[0160] S4.1: Neuron Division:

[0161] Grid size: 100×100m (adaptive adjustment);

[0162] Total number of neurons: N² = c (>5000);

[0163] Grid size: AU × AV m;

[0164] Coordinates: X i =(x i y i (UTM projection);

[0165] S4.2: Neuron Properties:

[0166] Static attribute matrix A i (Normalized to [0, 1]):

[0167] Ai = [Pˉ, Tˉ, Soil Code, Land Use Index, NDVI￣, S, A]

[0168] S4.3: Connecting topologies, establishing connections based on the D8 flow direction, and using the D-infinity algorithm in flat areas (slope < 0.5°);

[0169] S4.4: Example of a rainstorm scenario using dynamic activation functions (2022-07-23):

[0170]

[0171] When Rt = emm / h, f act →1 (Strong abortion);

[0172] When Rt = fmm / h, f act →-1 (Water storage is the dominant factor);

[0173] S4.5: Accuracy Verification:

[0174] RMSE = 1 / T ∑t = 1 / T (Q obst -Q simt )2=gm 3 / s;

[0175] The conceptual diagram of the dynamic neural network constructed in step S4 is referenced. Figure 4 .

[0176] The embodiments of the present invention have been described in detail above with reference to the accompanying drawings. However, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those skilled in the art, various changes can be made without departing from the spirit of the present invention.

Claims

1. A dynamic DEM spatial interpolation method, characterized in that, The process includes the following steps S1-S4, which complete the spatial interpolation of the DEM data and the hydrological simulation of the target watershed: Step S1: For the target watershed, collect raw DEM data, perform primary segmentation of the raw DEM data based on ArcGIS hydrological analysis toolkit, and extract natural sub-watershed units through depression filling, flow direction analysis, runoff accumulation calculation and minimum catchment area threshold method. Step S2: For natural sub-basin units, construct a topographic feature matrix, use an improved self-organizing mapping network for secondary classification, optimize clustering through dynamic grid mechanism and flow direction constraint, and generate hydrological response units through Voronoi segmentation with hydrological constraints and boundary morphology optimization. Establish a hydrological attribute database for hydrological response units. Step S3: Integrate multi-source data, including meteorological elements, geological data, land use type data, and vegetation parameters; supplement attribute interpolation for each multi-source data to construct an interpolation auxiliary parameter system; Step S4: Construct a dynamic neural network, divide the target watershed into a neuron grid, build topological connections based on the D8 flow direction, simulate hydrological response through a dynamic activation function, and train and optimize the parameters of the dynamic neural network using the root mean square error (RMSE) as the loss function. The specific steps of step S4 are as follows: Step S4.1: Divide the target watershed into a regular grid, treating each grid cell as a neuron; using a spatial matching algorithm between hydrological response units and the grid, define the neuron coordinates as follows: ； In the formula, R represents the set of real numbers. This represents the total number of grid cells after the entire target watershed is divided into regular grids. Let i be the coordinates of the i-th neuron. Let X be the coordinates of the i-th neuron in the X and Y directions in the unified projection coordinate system during step S3. Step S4.2: Each neuron contains dynamic states and static properties, defined as follows: ； In the formula, Represents the i-th neuron. , , The dynamic state of the i-th neuron is represented by water depth, hydraulic gradient, and conductivity, respectively. This represents the static attribute matrix of the i-th neuron, which contains the multi-source data obtained in step S3; Step S4.3: Determine the connection topology based on the D8 water flow direction. By calculating the slope of each grid and its 8 neighboring grids, the water flow direction is directed to the downstream grid with the steepest slope. For flat areas with a slope of <0.5°, the D-infinity algorithm is used to assist in determining the flow direction. The neuron at the watershed outlet is set as the global sink, receiving input from all downstream terminal neurons. Step S4.4: Construct the dynamic activation function as follows: ； In the formula, Represents a dynamic activation function. The hydrological conduction coefficient, For real-time rainfall intensity, This indicates that the activation intensity varies with the flow velocity; Step S4.5: Construct the loss function as follows: ； In the formula, Represents the loss function. This represents the total number of time steps in the hydrological simulation. This represents the actual flow rate observed at the i-th time point, in m³ / s; This represents the flow rate predicted by the model at time point i, in m³ / s.

2. The dynamic DEM spatial interpolation method according to claim 1, characterized in that, The specific steps of step S1 are as follows: Step S1.1: Depression filling process to eliminate depressions in the original DEM data; Step S1.2: Perform flow direction analysis based on the D8 algorithm to generate the D8 flow direction matrix; Step S1.3: Calculate the cumulative runoff to generate the runoff path network; Step S1.4: Automatically extract the boundaries of natural sub-basin units based on the minimum catchment area threshold method; Step S1.5: Output the vectorized set of natural sub-basin units U={u1,u2,...,u n }, where u1, u2, ..., u n This represents each natural sub-basin unit, where n is the total number of natural sub-basin units.

3. The dynamic DEM spatial interpolation method according to claim 1, characterized in that, The specific steps of step S2 are as follows: Step S2.1: Construct a terrain feature matrix for the natural sub-basin units extracted in Step S1: ； In the formula, This represents the topographic feature matrix of the i-th natural sub-basin unit; m represents the number of units in the study area whose topographic features were calculated or sampled; Step S2.2: Using an improved self-organizing map network clustering algorithm, dynamic grid division is performed on natural sub-basin units, and the elevation standard deviation is calculated. ,like Then the natural sub-basin unit is a flat area, and the grid size is reduced; if Then the natural sub-basin unit is a medium-complexity region, and the grid size is based on... The value is linearly interpolated between the minimum and maximum values; if If the natural sub-basin unit is a steep region, the grid size is increased, as shown in the following formula: ； In the formula, Indicates the dynamic mesh size. The number of grid cells within a natural sub-watershed unit; The dynamic mesh generation mechanism and the flow direction constraint function work together to improve distance calculation by introducing the flow direction constraint function: ； In the formula, This represents the flow direction constraint function. The weight of the k-th terrain feature is denoted as , where For elevation weighting, For slope weight, As aspect weight, As curvature weight, Weighting for topographic humidity index; This represents the value of the current raster at the k-th terrain feature. This represents the value of the k-th cluster center; For flow weight, The grid center direction angle. The angle of the water flow direction; Among them, the flow weight The calculation is as follows: ； In the formula, Indicates the cumulative amount of the confluence; Step S2.3: Output the clustering result set from the K winning neurons in the output layer of the improved self-organizing map network. : ； In the formula, This represents the k-th clustering result. The average elevation, To take into account the slope, As the dominant slope aspect, Characteristic curvature, This represents the average topographic humidity index. Step S2.4: Convert the clustering results into a set of hydrological response units using Voronoi diagrams. The hydrological response unit is as follows: ； In the formula, Represents the hydrological distance function. Represents the representative points of the clustering results; This represents the j-th clustering result; The hydrological distance function is as follows: ； In the formula, The dynamic weighting coefficients are between 0 and 1. Represents the Euclidean distance between two-dimensional geographic coordinates. This represents the cumulative distance along the D8 water flow path; Step S2.5: Extract the river network generated by the D8 algorithm As a hard constraint, specifically To constrain the re-delineation of boundaries for hydrological response units; Step S2.6: For areas smaller than threshold A min The hydrological response units are fused together, as shown in the following formula: ； In the formula, This indicates that the area is less than the threshold A. min Hydrological response unit; Indicates the candidate neighboring hydrological response unit; ΔH is the elevation difference, and ΔL is the flow path length; Indicates the distance of the water flow path between HRUs; Step S2.7: Optimize the boundary morphology of the fused hydrological response units by using morphological opening operations to eliminate microscopic jagged edges, as shown in the following formula: ； In the formula, Indicates the original hydrological boundary. This represents the morphologically optimized hydrological boundary. Represents a 3×3 structure element; For the morphologically optimized hydrological boundary, a Bézier curve fitting based on hydrological constraints is used, as shown in the following formula: ； In the formula, Indicates a smooth boundary in hydrology. This is the tolerance threshold. Indicates the resolution of the DEM data; Step S2.8: Establish a hydrological attribute database for the hydrological response units obtained in step S2.

7.

4. The dynamic DEM spatial interpolation method according to claim 3, characterized in that, The parameters in the hydrological attribute database established in step S2.8 include the maximum runoff path length, effective runoff gradient, spatial variability of soil moisture, uniformity of solar radiation distribution, infiltration capacity index, drainage density, and runoff coefficient.

5. The dynamic DEM spatial interpolation method according to claim 1, characterized in that, The specific steps of step S3 are as follows: Step S3.1: Collect hourly observation data from meteorological stations in the target watershed and surrounding preset range, and simultaneously introduce DEM data as a covariate to eliminate the terrain shading effect; generate meteorological element raster data based on the ANUSPLIN algorithm, including rainfall, temperature, and wind speed; Step S3.2: Interpolate the soil permeability coefficient using the co-kriging method. Use the exponential model for sandy soil and the Gaussian model for clay soil, and output the geological data. Step S3.3: Obtain land use type data by performing random forest classification based on Sentinel-2 satellite imagery; Step S3.4: Calculate the NDVI index based on Sentinel-2 imagery as a vegetation parameter, and remove noise through BRDF correction and time-series harmonic analysis.

Images