Dynamic DEM spatial interpolation method

Abstract

The invention discloses a dynamic DEM spatial interpolation method which comprises the following steps: performing depression filling, flow direction analysis and confluence cumulant calculation on DEM data, and extracting a natural sub-basin unit by adopting a minimum catchment area threshold method; constructing a topographic feature matrix, performing refined second-level classification on the first-level drainage basin through an improved self-organizing mapping network, generating a hydrological response unit through boundary processing, and establishing a hydrological attribute library; fusing multi-source data, supplementing attribute interpolation such as underlying surface and rainfall, and constructing an interpolation auxiliary parameter system; a drainage 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 achieved through optimization of a dynamic activation function and a loss function. The method overcomes the defects that a traditional interpolation algorithm does not consider hydrological boundary constraints, a neural network model topological structure is fixed and the like, and the DEM interpolation precision and the hydrological simulation effect of the complex terrain area are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geographic information systems and hydrological models, and particularly relates to a dynamic DEM spatial interpolation method, which is suitable for high-precision hydrological simulation, geological disaster warning and water resource management. BACKGROUND

[0002] As basic data for depicting the terrain undulations of the earth's surface, digital elevation models (DEM) are the core support in the fields of hydrological cycle simulation, watershed water resource management, geological disaster warning, and ecological pattern analysis. The precision of DEM directly determines the reliability of key topographic parameters such as river network extraction, slope and aspect calculation, and confluence path simulation, and further affects the scientificity of flood control engineering design, water resource allocation scheme, and ecological protection planning. For example, in the flood forecasting of small and medium-sized watersheds, an interpolation error of DEM of only 0.5 meters can lead to a confluence time deviation of more than 15%, significantly increasing the risk of flood peak warning inaccuracy; in the ecological research of permafrost regions, the distortion of terrain undulations will directly affect the simulation of permafrost ablation rate, and weaken the accuracy of hydrological process evaluation in cold regions.

[0003] The limitations of traditional spatial interpolation methods are particularly prominent in complex geographical environments. Taking the Kriging method as an example, its optimal unbiased estimation hypothesis based on spatial autocorrelation cannot be established in karst topography or plateau canyon regions with broken terrain. Due to differences in rock properties and fault cutting, the terrain in these regions is abruptly changed, and the spatial heterogeneity is severe. The semi-variogram function model of Kriging is prone to over-smoothing, which blurs key topographic features such as cliffs and gullies. The inverse distance weighting method is affected by the subjective setting of the distance decay coefficient. In mountainous areas with sparse data (such as the edge basins of the Qinghai-Tibet Plateau), the "isolated point domination" phenomenon often occurs due to unbalanced weight distribution, which causes the elevation values of adjacent sub-basins to be incorrectly flattened, destroying the hierarchical structure of the natural drainage system.

[0004] Although existing neural network models have advantages in non-linear fitting, they have the problem of "data-physical" disconnection. Most models use randomly sampled discrete elevation points as training samples, without considering the hydrological topological relationships in the basin (such as the confluence level of tributaries and main streams, and the hydrological isolation effect of drainage divides), resulting in a lack of physical rationality in the learned spatial patterns. For example, in soil erosion areas caused by rainfall, existing models cannot dynamically adjust the terrain weight based on real-time rainfall data, making it difficult to capture the effects of raindrop splash erosion and slope runoff on micro-topography, so that the updated DEM still remains at the level of static terrain representation, and the coupling degree with the actual hydrological process is insufficient.

[0005] At the level of spatial block, the conventional regular grid block (such as 1km*1km) will divide the complete slope land of the same slope into multiple grids in the cross mountainous area with dramatic terrain, resulting in fragmentation of slope length and slope gradient calculation; and the block based on administrative boundary (such as county unit) may split the complete basin into different management units, so that the simulation of confluence process appears 'cross-border fracture', and the natural continuity of hydrological process cannot be reflected. Such block method is contrary to the principle of 'terrain-soil-vegetation' synergy required by hydrological response unit (HRU), and directly leads to the increase of parameter uncertainty in distributed hydrological model.

[0006] In recent years, although some researches have tried to introduce hydrological constraints (such as correction of interpolation results based on river network), they are mostly limited to single scale basin, and have not formed a multi-level nested mechanism of 'whole basin-sub-basin-micro-basin', so it is difficult to balance the macro confluence trend and micro terrain details; and some dynamic interpolation models excessively rely on high-resolution observation data, and the applicability is limited in underdeveloped areas with scarce data. Therefore, to construct a DEM spatial interpolation method with physical constraints, dynamic response and multi-scale adaptability has become a key issue to solve the precision bottleneck of hydrological simulation in complex terrain area. SUMMARY

[0007] The present application aims to provide a dynamic DEM spatial interpolation method, which aims to solve the problems existing in the current DEM spatial interpolation method in the aspects of basin grading processing, attribute supplementing, neural network construction, etc., and to provide a dynamic DEM spatial interpolation method with higher precision and better dynamic response.

[0008] To achieve the above functions, the present application designs a dynamic DEM spatial interpolation method, which comprises the following steps S1-S4 to complete the DEM data spatial interpolation and hydrological simulation of the target basin:

[0009] Step S1: collecting DEM original data for the target basin, performing primary block on the DEM original data based on the ArcGIS hydrological analysis toolkit, and extracting natural sub-basin units through depression filling processing, flow direction analysis, confluence accumulation calculation and minimum catchment area threshold method;

[0010] Step S2: constructing a terrain feature matrix for the natural sub-basin unit, performing secondary classification by using an improved self-organizing mapping network, optimizing clustering by dynamic grid mechanism and flow direction constraint term, and generating hydrological response units through Voronoi segmentation with hydrological constraints and boundary shape optimization, and establishing a hydrological attribute library for the hydrological response units;

[0011] Step S3: fusing multi-source data, including meteorological elements, geological data, land use type data and vegetation parameters; respectively supplementing attribute interpolation for the multi-source data, and constructing an interpolation auxiliary parameter system;

[0012] Step S4: building a dynamic neural network, dividing the target basin into a neuron grid, constructing a topological connection based on the D8 water flow direction, simulating the hydrological response through a dynamic activation function, and training and parameter optimizing the dynamic neural network by taking the root mean square error (RMSE) as a loss function.

[0013] Advantages: Compared with the prior art, the advantages of the present application include:

[0014] The present application designs a dynamic DEM spatial interpolation method, which is not limited to a single scale basin, forms a multi-level nested mechanism of "whole basin-sub-basin-micro-basin", and takes into account the macro convergence trend and micro topographic details; a dynamic DEM spatial interpolation method with physical constraints, dynamic response and multi-scale adaptability is constructed, which overcomes the defects of traditional interpolation algorithms not considering hydrological boundary constraints and neural network model topological structure fixed, improves the DEM interpolation accuracy and hydrological simulation effect in complex terrain area, and is suitable for the technical field of high-precision hydrological simulation, geological disaster warning and water resources management. BRIEF DESCRIPTION OF DRAWINGS

[0015] Figure 1 is a flowchart of a dynamic DEM spatial interpolation method according to an embodiment of the present application;

[0016] Figure 2 is a DEM schematic diagram of a study area basin according to an embodiment of the present application;

[0017] Figure 3 is a schematic diagram of a first-order basin and a hydrological response unit according to an embodiment of the present application;

[0018] Figure 4 is a conceptual diagram of a dynamic neural network according to an embodiment of the present application. DETAILED DESCRIPTION

[0019] The present application will be further described below in conjunction with the drawings. The following embodiments are only used to more clearly illustrate the technical solutions of the present application, and cannot be used to limit the protection scope of the present application.

[0020] A dynamic DEM spatial interpolation method provided by an embodiment of the present application, with reference to Figure 1 , includes the following steps S1-S4, which completes the DEM data spatial interpolation and the hydrological simulation of the target basin:

[0021] Step S1: For the target basin, collect 30m resolution DEM (Digital Elevation Model) original data, based on the ArcGIS hydrological analysis toolkit to divide the DEM original data into a first block, and through the processing of filling the depression, flow direction analysis, convergence accumulation calculation and minimum catchment area threshold method, extract the natural sub-basin unit;

[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; representing the rate of change of elevation Z in the east-west direction; representing the rate of change of elevation Z in the north-south direction;

[0037] The curvature Curvature is calculated as follows:

[0038]

[0039] The terrain wetness index TWI is calculated as follows:

[0040]

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

[0042] Step S2.2: Using the improved self-organizing mapping network (SOM) clustering algorithm, dynamic grid division is performed for the natural sub-basin unit, and the grid size is adaptively adjusted according to the terrain complexity. The elevation standard deviation σ Z is calculated, and if σ Z < 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 ≤ 30, the natural sub-basin unit is a medium complexity area, the terrain complexity is between the flat and steep transition areas, and the algorithm uses a smooth adaptive strategy; the grid size is linearly interpolated between the minimum and maximum values according to the value of σ Z ; if σ Z > 30, the natural sub-basin unit is a steep area, the grid size is expanded to avoid over-segmentation, and the specific formula is as follows:

[0043]

[0044] In the formula, GRIDsize represents the dynamic grid size, and N1 is the number of grids in the natural sub-basin unit.

[0045] The dynamic grid division mechanism and the flow direction constraint function work together to introduce the flow direction constraint function to improve the distance calculation:

[0046]

[0047] In the formula, D' represents the flow direction constraint function, ω k , X k , C k The meanings, parameter definitions, value ranges, and hydrological meanings are summarized in Table 1:

[0048] Table 1. Meanings, parameter definitions, value ranges, and hydrological meanings of ω k , X k , C k in the flow direction constraint function

[0049]

[0050] Parameter definition, meaning, calculation method, and hydrological significance of λ, θ ij , φ in flow direction constraint function flow are summarized in Table 2 as follows:

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

[0052]

[0053] wherein the calculation of the flow direction weight λ is as follows:

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

[0055] In the formula, FlowAcc represents the flow accumulation;

[0056] River area (FlowAcc > 100): λ ≈ 0.4 (strong constraint). Slope area (FlowAcc < 5): λ < 0.1 (weak constraint). Among them, the terrain feature distance term calculates the difference between the current grid and the cluster center in five terrain features of elevation, slope, aspect, curvature and terrain wetness index (TWI) through weighted Euclidean distance, and the feature weight coefficients are shown in Table 1. The weight configuration highlights the dominant role of elevation zoning and soil moisture on hydrological processes. The flow direction constraint term introduces hydrological physical mechanism, and the dynamic weight formula is used to regulate the consistency strength of flow direction. When the flow accumulation (FlowAcc) increases (such as river area), λ tends to 0.4 to force the water flow direction to be consistent with the main flow direction determined by the D8 algorithm, and the λ value in the slope area (FlowAcc < 5) is less than 0.1 to weaken the flow direction constraint. The angle difference calculation eliminates the circumferential discontinuity and ensures that the 0° and 360° directions are equivalent. Through the coupling of terrain features and hydrological flow direction, the isohypse belt hydrological response unit (HRU) is formed in the mountainous area, the linear hydrological response unit is generated in the river area, and the soil moisture dominated hydrological response unit is constructed in the plain area, which significantly improves the consistency of water flow direction within the hydrological response unit.

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

[0058]

[0059] In the formula, P k represents the kth clustering result, E k is the average elevation, and S k ​for comprehensive slope, A k for dominant slope direction, C k for characteristic curvature, T k for terrain humidity index mean value; each neuron corresponds to a terrain homogeneous unit, and its position in the feature space determines the generating base point of the Voronoi diagram;

[0060] Step S2.4: performing Voronoi partitioning with hydrological constraints in a geographic information system (GIS), and converting the clustering result into a hydrological response unit set {HRU1, HRU2, …, HRU K} through the Voronoi diagram, where the hydrological response unit is as follows:

[0061]

[0062] In the formula, d hydro represents a hydrological distance function, p represents a representative point of the clustering result, and P j represents the jth 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 weight coefficient between 0 and 1, that is, α ∈ [0, 1], which is used to regulate the relative contribution of the two-dimensional geographic coordinate Euclidean distance and the cumulative distance along the D8 water flow path in the hydrological distance calculation; d Euclidean represents the two-dimensional geographic coordinate Euclidean distance, and d flow represents the cumulative distance along the D8 water flow path.

[0066] Step S2.5: extracting the river network Γ channel generated by the D8 algorithm as a hard constraint, and specifically taking as a constraint to redivide the boundary of the hydrological response unit, and using a river attractor algorithm;

[0067] Step S2.6: merging the hydrological response units with an area less than a threshold A min , and specifically as follows:

[0068]

[0069] In the formula, HRU small represents the hydrological response unit with an area less than a threshold A min , and HRU j represents a candidate adjacent hydrological response unit (that is, HRU small may be the target HRU to be merged); in an embodiment, Amin = 0.01 km 2 , AH is the height difference, AL is the length of flow path, and the downstream hydrological response unit (HRU) is preferentially fused to the terrain;

[0070] Step S2.7: Boundary morphology optimization is the core step of hydrological response unit generation, and the physical rationalization of the hydrological boundary is achieved through two-stage processing. The boundary morphology of the fused hydrological response unit is optimized, and the morphological opening operation is used to eliminate microscopic sawtooth, which is specifically as follows:

[0071]

[0072] In the formula, B raw represents the original hydrological boundary (Raw Hydrological Boundary), which refers to the initial boundary polygon directly output from the hydrological model, and there are pixel-level sawtooths (typical resolution ≤ 30 meters) caused by DEM discretization; B final represents the morphologically optimized hydrological boundary (Morphologically Optimized Boundary), which refers to the intermediate boundary after the dilation-erosion cascade processing; SE 3×3 represents a 3×3 structuring element (Structuring Element), which is defined as a square binary template of 90 meters × 90 meters (3 times the DEM pixel scale) with the center as the origin, and is used for uniform processing of the target pixel and its eight neighborhood relationships.

[0073] The above process uses a 3×3 structuring element to perform dilation operation on the original boundary to fill small fractures and gaps caused by DEM errors, and then performs equal erosion operation to remove isolated protrusions. This dilation-erosion cascade processing can effectively eliminate the terrain noise at the pixel scale (usually < 90 meters), reduce the boundary length by 12%-18%, and at the same time maintain the topological connectivity of the river channel and the slope surface.

[0074] For the morphologically optimized hydrological boundary, Bézier curve fitting based on hydrological constraints is adopted, which is specifically as follows:

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

[0076] In the formula, B smooth represents the hydrologically smoothed boundary (Hydrologically Smoothed Boundary), which refers to the final hydrological boundary optimized by the Bézier curve with flow direction constraint, and meets the sub-pixel accuracy and runoff direction consistency. Tolerance is the tolerance threshold, and δ represents the resolution of DEM data.

[0077] Control points are set every 60 meters (2 times DEM resolution) along the optimized boundary, and are densified to 15 meters interval at abrupt curvature. The core innovation of the reconstruction with cubic Bézier curve parameterization lies in the projection of control points along the flow direction: by associating D8 flow direction data, the angle between the curve tangent and the surface runoff direction is ensured to be ≤5°. A tolerance threshold of 0.5 times DEM resolution (15 meters) is specially set (i.e. 0.5δ in the formula), and when the fitting deviation exceeds this limit value, a new control point is automatically inserted to maintain the boundary accuracy at sub-pixel scale.

[0078] Step S2.8: For the hydrological response unit obtained in step S2.7, a hydrological attribute library is established as shown in Table 3:

[0079] Table 3. Hydrological attribute library of hydrological response unit

[0080]

[0081] Step S3: Fusion of multi-source data, including meteorological elements, geological data, land use type data, vegetation parameters; respectively for multi-source data supplementary attribute interpolation, construct interpolation auxiliary parameter system;

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

[0083] Step S3.1: Regarding meteorological elements, ANUSPLIN algorithm is used to interpolate to generate meteorological element grid data (time resolution ≤1 hour) such as rainfall, air temperature, wind speed, etc. The data input uses the hourly observation data of the target basin and the surrounding 30 km range of national meteorological stations, and introduces DEM data as a covariate to eliminate the terrain shielding effect. Before interpolation, data quality control is required: remove outliers beyond 3 times the standard deviation, and use linear interpolation method to fill in observation gaps ≤3 hours. Optimize the order of spline function (recommended 3 orders) and tension parameter (value range 0.01-0.5) through cross-validation method, to ensure that the average absolute error (MAE) of the interpolation result is: rainfall ≤5%, air temperature ≤0.5℃, wind speed ≤0.3m / s. Finally, output the meteorological element time series data set with the same resolution as the target DEM data (WGS84 coordinate system, UTM projection).

[0084] Step S3.2: Regarding the geological data, the soil permeability coefficient interpolation aspect, the data replacement uses the soil survey results, the permeability coefficient reference value of the 1:100,000 soil texture map (sandy soil 10-50 m / d, loamy soil 1-10 m / d, clay soil 0.1-1 m / d); the interpolation optimization uses the collaborative Kriging method, with soil texture as an auxiliary variable, the sandy soil uses the exponential model (spatial correlation range 2-5 km), the clay soil uses the Gaussian model (1-3 km), the block gold value is set to the theoretical value 10%-15%, and the comparison error with the measured value of the irrigation station is ≤15%. When interpolating the lithology type, the data integration is based on the 1:25,000 geological map to extract the lithology boundary, combined with the geological report to establish a 10-class first-level and 30-class second-level lithology coding system, and the geophysical prospecting data is introduced to enhance the identification of hidden lithology; the method adjustment uses the probability Kriging method, the lithology contact zone width is expanded to 100-200 m, the buffer zone weight is 1.2-1.5 times, and the type matching accuracy is ≥80% verified with the drilling data. On the quality control, the data of the authoritative agency in the past 15 years is preferred, 20% of the samples are cross-verified (when the error is >30%, supplement and correction), the output data is labeled with reliability level (A / B / C level), and the DEM data interpolation requirements (resolution matching ≥90%, type consistency ≥85%) are ensured.

[0085] Step S3.3: Regarding the land use type data, the land use type data is obtained by random forest classification based on Sentinel-2 satellite image (10 m resolution). The image selection needs to meet the cloud coverage rate <10% of the growing season data (June-August), the pretreatment includes atmospheric correction (Sen2Cor plug-in), terrain correction (C correction method) and radiation normalization. The feature variable set contains 10 band reflectances (B2-B8A, B11-B12) and 4 vegetation indices (EVI, SAVI, MSAVI, RVI), the hierarchical sampling method is used to divide the training samples (70%) and the verification samples (30%), and the classifier parameter setting is: the number of decision trees is 500, the minimum sample size of node splitting is 5, and the out-of-bag data (OOB) error is controlled within 8%.

[0086] Step S3.4: Regarding the vegetation parameter, the NDVI index (Normalized Difference Vegetation Index) is calculated as the core vegetation parameter. The near-infrared band (B8, 842 nm) and red band (B4, 665 nm) of the Sentinel-2 image are used to calculate the NDVI index by the formula NDVI = (B8-B4) / (B8+B4), with a value range of [-1, 1]. To eliminate the influence of atmospheric aerosols, the original band is first corrected by BRDF (Bidirectional Reflectance Distribution Function), and then the median filter (3x3 window) is used to smooth the noise. For seasonal vegetation changes, it is necessary to construct a monthly NDVI maximum value composite (MVC) product, and the cloud-polluted pixels are removed by time series harmonic analysis (HANTS). Finally, the generated NDVI dataset needs to be correlated with the contemporaneous vegetation field survey data (Leaf Area Index LAI) for correlation verification (R 2 ≥ 0.75).

[0087] Step S4: Build a dynamic neural network, divide the target basin into a grid of neurons, build a topological connection based on the D8 water flow direction, simulate the hydrological response through a dynamic activation function, and train and optimize the dynamic neural network with 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 basin into regular grids, such as 100x100m, and the grid size can be adaptively adjusted according to the original resolution of the DEM (recommended 30m or 10m), ensuring that the terrain slope variation coefficient within a single grid is ≤5%; regard each grid cell as a neuron; use the spatial matching algorithm of hydrological response unit and grid, when the hydrological response unit crosses multiple grids, take the center point of the grid where the hydrological response unit centroid is located as the neuron reference position, and the edge grid is ensured to completely cover the basin boundary through buffer analysis (width = 1 / 2 of the grid side length). 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 target basin is divided into regular grids (N2≥5000 is required to ensure the details of the basin), X i is the i-th neuron coordinate, (x i ,y i) is the coordinate of the i-th neuron in the X, Y direction of the WGS84-UTM projection coordinate system unified in step S3, which is batch calculated through the grid center point coordinates of DEM data, and the missing area (such as water area) is interpolated and completed by inverse distance weighting (IDW);

[0092] Step S4.2: Each neuron contains dynamic state and static attribute, which is defined as follows:

[0093]

[0094] wherein θ i represents the i-th neuron, h i , k i represents the dynamic state of the i-th neuron, which is water depth, hydraulic gradient, and hydraulic conductivity in turn; A i represents the static attribute matrix of the i-th neuron, which contains the multi-source data obtained in step S3, including meteorological elements, geological data, land use type data, and vegetation parameters; the static attribute matrix A i contains 8-dimensional features: the multi-year average of meteorological elements in step S3.1 (rainfall, air temperature, wind speed), the geological data in step S3.2, specifically the soil texture code, the land use type index in step S3.3, the average NDVI and slope, and the slope direction (calculated from DEM), all of which are normalized to the [0, 1] interval by min-max normalization to eliminate the dimension effect;

[0095] Step S4.3: Determine the connection topology based on the D8 flow direction (downstream neurons receive upstream input), and point the flow direction to the downstream grid with the largest slope by calculating the slope of each grid and the 8-neighborhood grid; for flat areas with a slope <0.5°, the D-infinity algorithm is used to assist in judging the flow direction to avoid topological fracture caused by ambiguous flow direction; the neurons at the outlet of the watershed are set as global sink points to receive the input of all downstream terminal neurons;

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

[0097]

[0098] wherein f act represents the dynamic activation function, β is the hydrological conduction coefficient (determined by soil type and slope), Rt is the real-time rainfall intensity (input through the meteorological interpolation module), represents that the activation strength changes with the flow velocity, and the activation strength changes nonlinearly with the flow velocity and rainfall intensity, simulating the real hydrological response. The physical meaning of the activation function is: when f act tends to 1, the neuron is in a strong flow concentration state (the grid runoff capacity is significantly improved); when fact When approaching -1, it is in a weak confluence state (mainly water storage), and the threshold response characteristics of rainfall-runoff are simulated by a nonlinear transformation.

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

[0100]

[0101] In the formula, Loss represents the loss function, N3 is the total number of time steps of hydrological simulation; Qobs(i) represents the actual flow value observed at the i th time point, with the unit of m 3 / s; Qpred(i) represents the flow value predicted (simulated) by the model at the i th time point, with the unit of m 3 / s.

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

[0103] The Tunxi River Basin is selected as the research area. The Tunxi River Basin, as an important part of the Xin'an River Basin, is located in the southeast of the Wannan Mountainous Area. The basin is located in the main area of the Huangshan Mountain, adjacent to the Yangtze River system in the west and north, and adjacent to the Tianmu Mountain and the Baiji Mountain in the southeast. There are multiple large-scale intermountain basins and valleys in the basin, and the overall terrain presents a west-high and east-low feature, with significant terrain undulations and generally large river channel slopes. In terms of climate type, it belongs to a typical subtropical monsoon climate zone, with a high annual precipitation and an average annual temperature in the transition zone between temperate and subtropical zones. The vegetation coverage in the basin is good, mainly including various forest types and agricultural land, forming a complex ecological system combination. The elevation changes widely, showing a significant decreasing trend from the western mountainous area to the eastern valley, forming a significant topographic difference.

[0104] In the implementation process, first of all, the basic data collection work is completed, covering hydrological and meteorological, topography and geomorphology, and social and economic information in multiple dimensions. The core data set includes real-time hydrological monitoring records, soil physical property parameters, precipitation observation data at different time scales, high-precision terrain elevation model, population distribution statistical data, and historical flood characteristic values, etc. These data are divided into two categories according to the processing requirements: routine processing set and to-be-processed set, corresponding to stable hydrological elements and dynamic change elements respectively. Data collection, DEM data: 30m resolution SRTM DEM, covering the basin and 5km buffer zone (matrix D∈R6150×5200D∈R6150×5200); meteorological data: 12 national meteorological stations near 6 years of hourly rainfall P, air temperature T, wind speed W; hydrological data: Tunxi station hourly flow Q; geological data: soil texture map (1:100,000): sandy soil 32%, loamy soil 45%, clay soil 23%, 25 borehole measured permeability coefficient Ks; remote sensing data: Sentinel-2 image (2020-2022 June-August, cloud cover <10%), 62 LAI measurement points (35 forest, 18 farmland, 9 grassland). The DEM diagram of the study area basin is referred to Figure 2 .

[0105] Subsequently, data standardization preprocessing is carried out, and all indicators are converted to the same dimension system. On this basis, the combined weighting method is used to determine the weight coefficient of each influencing factor: first, the judgment matrix is constructed through the expert decision model, and the subjective weight distribution is calculated; at the same time, based on the information entropy theory, the discrete characteristics of the data itself are analyzed, and the objective weight distribution scheme is derived. The two kinds of weight results are fused after weighting to form a comprehensive weight system, providing quantitative basis for subsequent analysis. This double weight mechanism not only integrates the field knowledge and experience, but also fully respects the objective law of data, ensuring the scientificity and reliability of the evaluation system.

[0106] In the spatial discretization processing stage, the grid division scheme is dynamically adjusted according to the terrain complexity, and the flow direction constraint mechanism is introduced to optimize the boundary of the hydrological response unit. Through morphological operation and curve fitting technology, the generated boundary is smoothed, which significantly improves the spatial continuity of the hydrological network. For extreme events such as heavy rain, a dynamic activation function is used to adjust the parameter sensitivity in real time, so that the model can adapt to the changes of different hydrological conditions. The whole process finally builds a multi-level nested hydrological discretization system, from sub-basin division to hydrological response unit definition to fine grid calculation, forming a hierarchical three-dimensional structure framework.

[0107] The specific calculation is as follows:

[0108] Step S1: ArcGIS processing to obtain the first-level classified basin, using ArcGIS hydrological analysis tool package processing:

[0109] 1. Filling: Repair M concave points in DEM;

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

[0111] 3. Sub-basin division: based on threshold F of flow accumulation;

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

[0113] Output: N natural sub-basins (area 82-168km 2 ).

[0114] Step S2:

[0115] S2.1: Take a typical sub-basin U i as an example, construct a spatial matrix containing 5-dimensional terrain features. First, use the third order inverse distance square weight algorithm to calculate the slope S, which accurately represents the steepness of the ground surface; solve the curvature Curvature by second derivative, describe the concave-convex characteristics of the terrain; based on the upstream catchment area As and the slope angle to derive the terrain wetness index TWI. Finally, form the terrain feature matrix:

[0116]

[0117] Slope using the third order inverse distance square weight algorithm:

[0118]

[0119] Curvature:

[0120]

[0121] Terrain wetness index (TWI):

[0122]

[0123] Where A s is the upstream catchment area, β is the slope angle. This matrix integrates 5 types of terrain attributes of F grids, providing data basis for clustering analysis.

[0124] S2.2: Fusion of dynamic grid mechanism and hydrological flow direction constraint:

[0125] Grid self-adaptation: dynamically adjust the clustering resolution according to the elevation standard deviation σ Z = a, encrypt the grid in flat areas (a < 10), expand the grid in steep areas (a > 30), and calculate the dynamic grid size:

[0126] (original 30m→ adjusted bm)

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

[0128]

[0129] Where flow direction weight λ = 0.4 × (1 - e -0.1×F ) varies dynamically with the accumulation of confluence:

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

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

[0132] Where c is the slope area flow direction weight expression; e is the river area flow direction weight expression; d is the slope area flow direction weight value; f is the river area flow direction weight value;

[0133] This design ensures that the mountainous area generates contour belt HRU, and the river area forms linear HRU.

[0134] S2.3: Topographic prototype generation

[0135] Cluster output K topographic 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 slope direction a1, characteristic curvature c1 and TWI average t1 of a certain type of homogeneous unit of topography. These prototypes constitute the reference points of HRU division.

[0136] Output K topographic 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) = (632 m, 8.7°, northeast, -0.018, 12.3)

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

[0140] 1. Hydrological constraint Voronoi partition: Define the hydrological distance function d hydro = α · d Euclidean + (1 - α) · d flow , which combines geographic distance d Euclidean and water flow path distance dflow ;

[0141] 2. Topology optimization: Reconstruct boundaries with river network Γ channel as hard constraints using river attractor algorithm;

[0142] 3. Micro-unit fusion: Perform downstream fusion on HRUs with area smaller than 0.01 km2 Preferably incorporate downstream units with terrain continuity. Fuse D HRUs with area A < 0.01 km 2 ;

[0143] 4. Boundary optimization:

[0144] Morphological processing: Eliminate jagged boundaries qkm;

[0145] Bézier curve fitting: B smooth = BezierFit(B final , 0.5δ) Keep water flow direction consistency (tangential angle ≤ 5°) Final output hydrological response unit set {HRU1, HRU2, …, HRU K}; First-order basin and hydrological response unit schematic diagram refer to Figure 3 .

[0146] S2.8: Hydrological attribute library with HRU_5 as an example as shown in Table 4:

[0147] Table 4. Hydrological attribute library (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] Soil permeability coefficient spatial distribution:

[0154] Lithology 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 x 100 m (adaptive adjustment);

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

[0163] Grid size: AU x AV m;

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

[0165] S4.2: Neuron attributes:

[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: Connection topology, based on D8 flow direction to establish connection, flat area (slope <0.5°) using D-infinity algorithm;

[0169] S4.4: Dynamic activation function Rainstorm scenario example (2022-07-23):

[0170]

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

[0172] When Rt = f mm / h, f act → -1 (dominated by water storage);

[0173] S4.5: Accuracy verification:

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

[0175] The concept diagram of the dynamic neural network built in step S4 is referred to Figure 4 .

[0176] The embodiments of the present application are described in detail above with reference to the accompanying drawings, but the present application is not limited to the above-described embodiments, and various changes can be made within the knowledge of those skilled in the art without departing from the spirit of the present application.

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 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.

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, T i 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. Z If σ Z If σ < 10, the natural sub-basin unit is a flat area, and the grid size is reduced; if 10 ≤ σ Z If the value is ≤30, then the natural sub-basin unit is a medium-complexity region, and the mesh 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, as shown in the following formula: In the formula, GRIDsize represents the dynamic grid size, and N1 is the number of grid cells in the natural sub-basin cell; 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, D′ represents the flow direction constraint function, ω k X represents the weight of the k-th terrain feature, where ω1 is the elevation weight, ω2 is the slope weight, ω3 is the aspect weight, ω4 is the curvature weight, and ω5 is the terrain humidity index weight; k C is the value of the current raster at the k-th terrain feature. k This represents the value of the k-th cluster center; λ is the flow direction weight, θ ij φ is the grid-center orientation angle. flow The angle of the water flow direction; The flow weight λ is calculated as follows: λ=0.4×(1-e -0.1×FlowAcc ) In the formula, FlowAcc represents the cumulative flow amount; Step S2.3: Output the clustering result set P from the K winning neurons in the output layer of the improved self-organizing map network: In the formula, P k E represents the k-th clustering result. k S is the average elevation. k For the overall slope, A k As the dominant slope aspect, C k For characteristic curvature, T k This represents the average topographic humidity index. Step S2.4: Convert 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: In the formula, d hydro Let P represent the hydrological distance function, and p represent the representative point of the clustering result; j This represents the j-th clustering result; The hydrological distance function is as follows: d hydro =α·d Euclidean +(1+a)·d flow In the formula, α is a dynamic weighting coefficient between 0 and 1, and d Euclidean d represents the Euclidean distance between two-dimensional geographic coordinates. flow This represents the cumulative distance along the D8 water flow path; Step S2.5: Extract the river network Γ generated by the D8 algorithm channel 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, HRU small This indicates that the area is less than the threshold A. min Hydrological Response Unit; HRU j Indicates the candidate neighboring hydrological response unit; ΔH is the elevation difference, ΔL is the flow path length; D flow 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, B raw B represents the original hydrological boundary. final SE represents the morphologically optimized hydrological boundary. 3×3 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: B smooth =BezierCurveFit(B final ,tolerance=0.5δ) In the formula, B smooth δ represents the hydrological smooth boundary, tolerance is the tolerance threshold, and δ represents the DEM data resolution. 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.

6. The dynamic DEM spatial interpolation method according to claim 1, characterized in that, 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: X i =(x i ,y i )∈R 2 ,i=1,2,…,N2 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, and 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 projection coordinate system during step S3; Step S4.2: Each neuron contains dynamic states and static properties, defined as follows: 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 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, f act This represents the dynamic activation function, where β is the hydrological conduction coefficient and Rt is the 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, 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 flow rate predicted by the model at time point i, in meters. 3 / s.

Images