CN120974910A

CN120974910A - A dynamic DEM spatial interpolation method

Info

Publication number: CN120974910A
Application number: CN202511101645.9A
Authority: CN
Inventors: 晁丽君; 冉婵娟; 侯婷婷; 李思颖; 邓瑶; 韩笑; 王飞
Original assignee: Hohai University HHU
Current assignee: Hohai University HHU
Priority date: 2025-08-07
Filing date: 2025-08-07
Publication date: 2025-11-18

Abstract

This invention discloses a dynamic DEM spatial interpolation method. It performs depression filling, flow direction analysis, and runoff accumulation calculation on DEM data, and extracts natural sub-basin units using the minimum catchment area threshold method. A topographic feature matrix is constructed, and a refined secondary classification of the first-order watershed is performed using an improved self-organizing mapping network. Hydrological response units are generated after boundary processing, and a hydrological attribute database is established. Multi-source data is fused, and interpolation of underlying surface, rainfall, and other attributes is supplemented to construct an interpolation auxiliary parameter system. The watershed is divided into regular grids as neurons, and a dynamic neural network containing dynamic states and static attributes is constructed. The nonlinear mapping of DEM correction parameters is achieved through optimization using dynamic activation and loss functions. This method overcomes the shortcomings of traditional interpolation algorithms, such as not considering hydrological boundary constraints and having a fixed topology of the neural network model, thus improving the accuracy of DEM interpolation and the hydrological simulation effect in complex topographic areas.

Description

Dynamic DEM spatial interpolation method

Technical Field

The invention relates to the technical field of intersection of a geographic information system and a hydrological model, in particular to a dynamic DEM spatial interpolation method which is suitable for high-precision hydrologic simulation, geological disaster early warning and water resource management.

Background

The Digital Elevation Model (DEM) is used as basic data for describing the relief of the surface topography of the earth, and is a core support in the fields of hydrologic cycle simulation, watershed water resource management, geological disaster early warning, ecological pattern analysis and the like. The accuracy directly determines the reliability of key terrain parameters such as river network extraction, gradient slope calculation, converging path simulation and the like, and further influences the scientificity of flood control engineering design, water resource allocation scheme and ecological protection planning. For example, in flood forecast of small and medium-sized watersheds, the DEM interpolation error is only 0.5m, so that the convergence time deviation is more than 15%, the flood peak early warning misalignment risk is obviously increased, in ecological research of frozen soil areas, the topography relief and characterization distortion can directly influence the simulation of the frozen soil ablation rate, and the evaluation accuracy of hydrologic processes of the cold areas is weakened.

The limitations of conventional spatial interpolation methods are particularly pronounced in complex geographic environments. Taking the kriging method as an example, the kriging method is based on the optimal unbiased estimation assumption of the space autocorrelation, and is difficult to be established in karst landforms or plateau canyon areas with broken terrains, wherein the areas have abrupt terrains due to lithology differences and fault cutting, the space heterogeneity is severe, and the semi-variation function model of the kriging method is easy to be excessively smooth, so that key terrains such as scarps and gullies are blurred. The inverse distance weight rule is influenced by subjective setting of a distance attenuation coefficient, and in mountain areas with sparse data (such as the Qinghai-Tibet plateau edge watershed), the phenomenon of 'isolated point dominance' often occurs due to unbalanced weight distribution, so that the elevation values of adjacent sub-watersheds are erroneously leveled, and the hierarchical structure of a natural water system is damaged.

The existing neural network model has the advantage of nonlinear fitting, but has the problem of data-physical disjoint. Most models take randomly sampled discrete elevation points as training samples, and do not incorporate watershed hydrologic topological relations (such as confluence level of tributary-main flow and hydrologic isolation effect of watershed), so that the learned spatial mode lacks physical rationality. For example, in a soil erosion area caused by rainfall, the terrain weight of the existing model cannot be dynamically adjusted through real-time rainfall data, and the transformation effect of raindrop splash erosion and slope runoff on micro-terrains is difficult to capture, so that the updated DEM still stays at a static terrain characterization level, and the coupling degree with an actual hydrologic process is insufficient.

On the space partitioning level, a conventional regular grid partition (such as 1km multiplied by 1 km) can divide a complete sloping field with the same slope direction into a plurality of grids in a transverse mountain area with severe fluctuation of the terrain to cause the calculation fragmentation of the slope length and the slope, and a partition (such as a county level unit) based on an administrative boundary can split the complete river basin into different management units to cause the simulation of 'boundary crossing fracture' in the converging process and can not reflect the natural continuity of the hydrologic process. This blocking approach is contrary to the principle that the Hydrologic Response Unit (HRU) needs to follow the "terrain-soil-vegetation" synergetics, directly resulting in an increase in the uncertainty of the parameters in the distributed hydrologic model.

In recent years, although research attempts are made to introduce hydrologic constraints (such as correcting interpolation results based on river networks), the method is limited to a single-scale river basin, a multi-stage nesting mechanism of 'full-river basin-sub-river basin-micro-river basin' is not formed, macroscopic convergence trend and microscopic topography details are difficult to be considered, and a part of dynamic interpolation model is excessively dependent on high-resolution observation data, so that the method has limited applicability in underdeveloped areas with data scarcity. Therefore, the construction of the DEM spatial interpolation method with physical constraint, dynamic response and multi-scale adaptability has become a key subject for breaking the hydrologic simulation precision bottleneck of the complex terrain area.

Disclosure of Invention

The invention aims to provide a dynamic DEM spatial interpolation method, which aims to overcome the defects of the existing DEM spatial interpolation method in terms of river basin grading treatment, attribute supplementation, neural network construction and the like, and provides a dynamic DEM spatial interpolation method with higher precision and better dynamic response.

In order to achieve the above functions, the invention designs a dynamic DEM spatial interpolation method, which comprises the following steps S1-S4, and completes the spatial interpolation of DEM data and the hydrologic simulation of a target river basin:

step S1, aiming at a target drainage basin, acquiring DEM original data, performing primary blocking on the DEM original data based on an ArcGIS hydrological analysis kit, and extracting natural sub-drainage basin units through filling processing, flow direction analysis, confluence accumulation calculation and a minimum water collection area threshold method;

step S2, constructing a topography feature matrix aiming at a natural sub-drainage basin unit, adopting an improved self-organizing mapping network to carry out secondary classification, carrying out optimized clustering through a dynamic grid mechanism and a flow direction constraint item, generating a hydrologic response unit through Voronoi segmentation with hydrologic constraint and boundary form optimization, and establishing a hydrologic attribute library aiming at the hydrologic response unit;

Step S3, fusing multi-source data, including meteorological elements, geological data, land utilization type data and vegetation parameters, respectively interpolating the multi-source data by supplementing attribute to construct an interpolation auxiliary parameter system;

and S4, building a dynamic neural network, dividing a target river basin into neuron grids, building topological connection based on the D8 water flow direction, simulating hydrologic response through a dynamic activation function, and training and optimizing parameters of the dynamic neural network by taking Root Mean Square Error (RMSE) as a loss function.

The beneficial effects of the invention are that compared with the prior art, the invention has the following advantages:

The invention designs a dynamic DEM spatial interpolation method, which is not limited to a single-scale drainage basin, forms a multi-level nesting mechanism of 'full drainage basin-sub drainage basin-micro drainage basin', gives consideration to macroscopic convergence trend and microscopic topography detail, and is constructed to have physical constraint, dynamic response and multi-scale adaptability.

Drawings

FIG. 1 is a flow chart of a dynamic DEM spatial interpolation method provided in accordance with an embodiment of the present invention;

fig. 2 is a schematic view of a study area basin DEM provided according to an embodiment of the invention;

FIG. 3 is a schematic diagram of a primary basin and hydrological response unit provided in accordance with an embodiment of the present invention;

fig. 4 is a conceptual diagram of a dynamic neural network provided according to an embodiment of the present invention.

Detailed Description

The invention is further described below with reference to the accompanying drawings. The following examples are only for more clearly illustrating the technical aspects of the present invention, and are not intended to limit the scope of the present invention.

The embodiment of the invention provides a dynamic DEM spatial interpolation method, referring to FIG. 1, comprising the following steps S1-S4, finishing the spatial interpolation of DEM data and the hydrologic simulation of a target river basin:

Step S1, acquiring 30m resolution DEM (digital elevation model) original data aiming at a target river basin, performing primary blocking on the DEM original data based on ArcGIS hydrologic analysis tool kit, and extracting natural sub-river basin units through filling processing, flow direction analysis, confluence accumulation calculation and a minimum water collection area threshold method;

the specific steps of the step S1 are as follows:

Step S1.1, filling processing (Fill) to eliminate DEM original data pits;

s1.2, carrying out Flow Direction analysis (Flow Direction) according to a D8 algorithm to generate a D8 water Flow Direction matrix;

S1.3, calculating the confluence cumulative quantity (Flow Accumulation) to generate a runoff path network;

S1.4, automatically extracting the boundary of a natural sub-drainage basin unit based on a minimum water collecting area threshold method (threshold range: 1-5km ²);

And S1.5, outputting a vectorized natural sub-basin unit set U= { U ₁,u₂,...,u_n }, wherein U ₁,u₂,...,u_n represents each natural sub-basin unit, and n is the total number of the natural sub-basin units.

Step S2, constructing a topography feature matrix aiming at a natural sub-drainage basin unit, adopting an improved self-organizing mapping network (SOM) to carry out secondary classification, carrying out optimized clustering through a dynamic grid mechanism and a flow direction constraint item, generating a hydrologic response unit through Voronoi segmentation with hydrologic constraint and boundary form optimization, and establishing a hydrologic attribute library aiming at the hydrologic response unit;

the specific steps of the step S2 are as follows:

step S2.1, constructing a topography feature matrix aiming at the natural sub-drainage basin unit extracted in the step S1:

Where T _i represents the terrain feature matrix of the ith natural sub-basin cell and m represents the number of cells in the area of investigation where the terrain feature is calculated or sampled.

The elevation adopts an original elevation, and the slope direction adopts azimuth quantization improvement;

slope uses a third-order inverse distance square weighting algorithm:

wherein x, y and Z respectively represent the east-west direction, the north-south direction and the elevation of the geographic coordinate system; Representing the change rate of the elevation Z in the east-west direction; Representing the change rate of the elevation Z in the north-south direction;

The curvature Curvature is calculated as follows:

The topographic humidity index TWI is calculated as follows:

wherein A _s is the upstream catchment area and β is the slope angle.

Step S2.2, adopting an improved self-organizing map network (SOM) clustering algorithm to carry out dynamic grid division on natural sub-basin units, carrying out self-adaptive adjustment on grid size according to terrain complexity, calculating an elevation standard deviation sigma _Z, if sigma _Z is smaller than 10, enabling the natural sub-basin units to be flat areas, reducing the grid size and improving local terrain resolution, if sigma _Z is smaller than or equal to 10 and smaller than or equal to 30, enabling the natural sub-basin units to be medium complexity areas, enabling the algorithm to adopt a smooth self-adaptive strategy, enabling the grid size to carry out linear interpolation between a minimum value and a maximum value according to the numerical value of sigma _Z, and if sigma _Z is larger than 30, enabling the natural sub-basin units to be steep areas, and avoiding over-segmentation, wherein the algorithm comprises the following specific formula:

wherein GRIDsize represents a dynamic grid size, and N ₁ is the grid number in a natural sub-basin unit;

the dynamic meshing mechanism cooperates with the flow direction constraint function, and the flow direction constraint function is introduced to improve the distance calculation:

Wherein D' represents a flow direction constraint function, and the meanings, parameter definitions, value ranges and hydrologic meanings of omega _k、X_k、C_k are summarized in the following table 1:

TABLE 1 meanings of omega _k、X_k、C_k, parameter definitions, value ranges, hydrologic meanings in flow constraint functions

The parameter definitions, meanings, calculation modes and hydrologic meanings of lambda and theta _ij、φ_flow in the flow direction constraint function are summarized in the following table 2:

TABLE 2 parameter definition, meaning, calculation method, hydrologic meaning of lambda, theta _ij、φ_flow in flow constraint function

Wherein, the flow direction weight lambda is calculated as follows:

λ=0.4×(1-e^{-0.1×FlowAcc})

Wherein FlowAcc represents a confluent cumulative amount;

River region (FlowAcc > 100) lambda is about 0.4 (strong constraint). Slope area (FlowAcc < 5) lambda <0.1 (weak constraint). The difference of the current grid and the clustering center in the five topographic features of elevation, gradient, slope direction, curvature and topographic humidity index (TWI) is calculated through the weighted Euclidean distance, the feature weight coefficients are respectively shown in the table 1, and the weight configuration highlights the leading effect of elevation banding and soil humidity on the hydrologic process. The flow direction constraint term introduces a hydrophysics mechanism, the flow direction consistency intensity is regulated and controlled through a dynamic weight formula, when the confluence accumulation amount (FlowAcc) is increased (such as a river channel area), lambda approaches 0.4 to force the water flow direction to be consistent with the main flow direction determined by the D8 algorithm, and the lambda value of the slope surface area (FlowAcc < 5) is lower than 0.1 to weaken the flow direction constraint. And the circumferential discontinuity is eliminated by adopting the calculation of the angle difference, so that the 0-degree direction is equivalent to the 360-degree direction. By coupling the topography features and the hydrologic flow direction, a contour Hydrologic Response Unit (HRU) is formed in a mountain area, a river channel area generates a linear hydrologic response unit, a hydrologic response unit with the soil humidity leading function is constructed in a plain area, and the uniformity of the water flow direction inside the hydrologic response unit is obviously improved.

Step S2.3, outputting a clustering result set P based on K winning neurons of an improved self-organizing map (SOM) output layer:

Wherein P _k represents a kth clustering result, E _k is an average elevation, S _k is a comprehensive gradient, A _k is a main slope direction, C _k is a characteristic curvature, and T _k is a terrain humidity index mean value;

Step S2.4, executing Voronoi segmentation with hydrologic constraint in a Geographic Information System (GIS), and converting a clustering result into a hydrologic response unit set { HRU ₁,HRU₂,…,HRU_K } through a Voronoi diagram, wherein the hydrologic response unit has the following formula:

wherein d _hydro represents a hydrological distance function, P represents a representative point of a clustering result, and P _j represents a j-th clustering result;

the hydrological distance function is as follows:

d_hydro＝α·d_Euclidean+(1+α)·d_flow

Wherein, alpha is a dynamic weight coefficient between 0 and 1, namely alpha E [0,1] which is used for regulating and controlling the relative contribution of the two-dimensional geographic coordinate Euclidean distance and the accumulated distance along the D8 water flow path in hydrological distance calculation, D _Euclidean represents the two-dimensional geographic coordinate Euclidean distance, and D _flow represents the accumulated distance along the D8 water flow path;

S2.5, extracting a river channel network Γ _channel generated by a D8 algorithm as hard constraint, specifically using Adopting a river channel attractor algorithm for restricting the repartition boundary of the hydrological response unit;

step S2.6, fusing the hydrological response units with the area smaller than the threshold A _min, wherein the specific formula is as follows:

Where HRU _small represents a hydrographic response unit with an area less than threshold A _min and HRU _j represents a candidate neighboring hydrographic response unit (i.e., target HRU where HRU _small may be fused), in one embodiment, A _min＝0.01km² is set to a height difference, ΔH is the flow path length, and ΔL is preferentially fused to a topography continuous downstream Hydrographic Response Unit (HRU);

step S2.7, boundary form optimization is a core link generated by the hydrologic response unit, physical rationalization of hydrologic boundaries is realized through two-stage processing, boundary form optimization is carried out on the fused hydrologic response unit, and microscopic saw teeth are eliminated by adopting morphological open operation, wherein the specific formula is as follows:

Wherein B _raw represents an original hydrological boundary (Raw Hydrological Boundary) which refers to an initial boundary polygon directly output from a hydrological model, pixel-level saw teeth (with typical resolution less than or equal to 30 meters) caused by DEM discretization exist, B _final represents a hydrological boundary (Morphologically Optimized Boundary) after morphological optimization and refers to an intermediate boundary after expansion-corrosion cascade treatment, SE _3×3 represents a 3×3 structural element (Structuring Element), a square binary template defined as 90 m×90 m (3 times DEM pixel scale) is defined, and the center is an origin for uniformly treating a target pixel and an eight neighborhood relation thereof;

The above process uses 3 x3 structural elements to expand the original boundary, fill in small breaks and gaps caused by DEM errors, and then perform an equivalent etching operation to remove isolated protrusions. The expansion-corrosion cascade treatment can effectively eliminate topographic noise with pixel dimensions (generally <90 meters), reduce the boundary length by 12% -18%, and simultaneously maintain the topological connectivity of the river channel and the slope.

Aiming at the hydrologic boundary after morphological optimization, a Bezier curve fitting based on hydrologic constraint is adopted, and the following specific formula is as follows:

B_smooth＝BezierCurveFit(B_final,tolerance=0.5δ)

Wherein B _smooth represents a hydrologic smoothing boundary (Hydrologically Smoothed Boundary), which refers to a final hydrologic boundary optimized by a Bezier curve constrained by flow direction, meets consistency of sub-pixel precision and flow direction, tolerance is a tolerance threshold, and delta represents DEM data resolution;

Control points are set every 60 meters (2 times DEM resolution) along the optimization boundary, and are encrypted to 15 meters spacing at the curvature abrupt change. The three Bezier curve parameterized reconstruction is adopted, and the core innovation is that the control point projects along the water flow direction, namely, through correlating D8 flow direction data, the included angle between a curve tangent and the surface runoff direction is ensured to be less than or equal to 5 degrees. A tolerance threshold (i.e. 0.5 delta in the formula) of 0.5 times DEM resolution (15 meters) is specifically set, and when the fitting deviation exceeds this limit, a new control point is automatically inserted, and the boundary accuracy is maintained at the sub-pel scale.

Step S2.8, establishing a hydrologic attribute library according to the hydrologic response unit obtained in the step S2.7, wherein the hydrologic attribute library is shown in the following table 3:

TABLE 3 hydrologic attribute library of hydrologic response units

The specific steps of step S3 are as follows:

And S3.1, regarding the meteorological elements, generating grid data of the meteorological elements such as rainfall, air temperature, wind speed and the like (time resolution is less than or equal to 1 hour) based on interpolation of an ANUSPLIN algorithm. The data input adopts time-by-time observation data of a target river basin and a national weather station within the range of 30km around, and meanwhile, DEM data is introduced as covariates to eliminate the topographic shading effect. And before interpolation, data quality control is required to be carried out, namely abnormal values exceeding 3 times of standard deviation are removed, and a linear interpolation method is adopted to fill in the observation missing less than or equal to 3 hours. The spline function order (recommended 3-order) and the tension parameter (value range is 0.01-0.5) are optimized through a cross validation method, so that the average absolute error (MAE) of interpolation results is ensured, wherein rainfall is less than or equal to 5%, air temperature is less than or equal to 0.5 ℃ and air speed is less than or equal to 0.3m/s. And finally outputting a meteorological element time sequence data set (WGS 84 coordinate system, UTM projection) with the same resolution as the target DEM data.

And S3.2, regarding geological data and soil permeability coefficient interpolation, replacing the data with permeability coefficient reference values (sandy soil 10-50m/d, loam 1-10m/d and clay 0.1-1 m/d) of a soil general investigation result and a 1:10 ten thousand soil map, adopting a collaborative Kriging method for interpolation optimization, taking soil texture as an auxiliary variable, an index model (a space correlation range is 2-5 km) for sandy soil, a Gaussian model (1-3 km) for clay, setting a block gold value as a theoretical value of 10% -15%, and comparing an error with an actual measurement value of an irrigation station to be less than or equal to 15%. During lithology type interpolation, data integration extracts lithology boundaries based on a 1:2.5 geological map, a 10-class primary and 30-class secondary lithology coding system is established by combining geological reports, geophysical prospecting data is introduced to enhance hidden lithology recognition, a probability kriging method is adopted for adjustment, lithology contact zone width is expanded to 100-200m, buffer zone weight is 1.2-1.5 times, and the lithology contact zone width is verified with drilling data, so that type matching accuracy is more than or equal to 80%. In the quality control aspect, the authority mechanism is preferably selected for data of 15 years, 20% of samples are cross-validated (the error is supplemented and corrected when the error is more than 30%), the reliability grade (A/B/C grade) of output data is marked, and the requirement of DEM data interpolation is met (the resolution ratio is more than or equal to 90%, and the type consistency is more than or equal to 85%).

And S3.3, regarding the land utilization type data, carrying out random forest classification based on the Sentinel-2 satellite image (10 m resolution) to obtain the land utilization type data. Image selection requires growth season data (6-8 months) meeting cloud cover <10%, and preprocessing includes atmospheric correction (Sen 2Cor plug-in), topography correction (C correction method) and radiation normalization. The characteristic variable set comprises 10 wave band reflectivities (B2-B8A, B-B12) and 4 vegetation indexes (EVI, SAVI, MSAVI, RVI), a layered sampling method is adopted to divide training samples (70%) and verification samples (30%), classifier parameters are set to be 500 decision trees, the minimum sample quantity of node splitting is 5, and the out-of-bag data (OOB) error is controlled within 8%.

Step S3.4, calculating an NDVI index (normalized difference vegetation index) as a core vegetation parameter with respect to the vegetation parameter. The near infrared band (B8, 842 nm) and the red light band (B4, 665 nm) of the Sentinel-2 image are adopted, and the value range is [ -1,1] is calculated through the formula NDVI= (B8-B4)/(B8+B4). To eliminate the atmospheric aerosol effect, BRDF (bi-directional reflectance distribution function) correction is performed on the original band, and then median filtering (3×3 window) is used to smooth the noise. For seasonal vegetation change, a monthly NDVI maximum synthesis (MVC) product needs to be constructed, cloud pollution pixels are removed through time sequence Harmonic Analysis (HANTS), and finally the generated NDVI data set needs to be subjected to correlation verification with field investigation data (leaf area index LAI) of the contemporaneous vegetation (R ² is more than or equal to 0.75).

The specific steps of step S4 are as follows:

Step S4.1, dividing a target river basin into regular grids, such as 100 multiplied by 100m, wherein the size of the grids can be adaptively adjusted according to the original resolution of the DEM (30 m or 10m is recommended), so that the gradient variation coefficient of the terrain in a single grid is ensured to be less than or equal to 5%, each grid unit is regarded as a neuron, a space matching algorithm of a hydrologic response unit and the grids is adopted, when the hydrologic response unit spans multiple grids, the center point of the grid where the centroid of the hydrologic response unit is located is taken as a reference position of the neuron, and the edge grids are analyzed through a buffer area (the width=1/2 of the side length of the grid) so as to ensure complete coverage of the boundary of the river basin. Neuron coordinates are defined as:

X_i＝(x_i,y_i)∈R²,i＝1,2,…,N₂

Wherein R represents a real number set, N ₂ is the total number of grid units (N ₂ is more than or equal to 5000 to ensure the description of the river basin details to be met) after the whole target river basin is divided into regular grids, X _i is the coordinates of the ith neuron, X _i,y_i is the coordinates of the ith neuron in the X, Y direction in a WGS84-UTM projection coordinate system unified in the step S3, grid center point coordinates of DEM data are calculated in batches, and a missing area (such as a water area) is complemented by interpolation by adopting an inverse distance weighting method (IDW);

Step S4.2, each neuron contains dynamic state and static attribute, which is defined as follows:

wherein θ _i represents the ith neuron, h _i, The dynamic state of the ith neuron is represented by k _i and is sequentially represented by water depth, hydraulic gradient and water conductivity coefficient, the A _i represents a static attribute matrix of the ith neuron and comprises multisource data obtained in the step S3, including meteorological elements, geological data, land utilization type data and vegetation parameters, the static attribute matrix A _i comprises 8-dimensional characteristics, namely a multi-year mean value (rainfall, air temperature and wind speed) of the meteorological elements in the step S3.1, geological data in the step S3.2, in particular soil texture coding, a land utilization type index in the step S3.3, an NDVI mean value and gradient in the step S3.4 and a slope direction (calculated by DEM), and all the attributes need to be normalized to a [0,1] interval through min-max to eliminate dimensional influence;

S4.3, determining connection topology (a downstream neuron receives upstream input) based on the D8 water flow direction, and directing the water flow direction to a downstream grid with the largest gradient by calculating the gradient of each grid and 8 neighborhood grids, wherein aiming at a flat area with the gradient less than 0.5 DEG, a D-infinity algorithm is adopted to assist in judging the flow direction, so that topology fracture caused by flow direction ambiguity is avoided;

step S4.4, constructing a dynamic activation function as follows:

where f _act denotes a dynamic activation function, β is a hydrologic conductivity (determined by soil type and gradient), rt is a real-time rainfall intensity (input by a weather interpolation module), And the activation intensity is expressed as the converging speed, the activation intensity is in nonlinear change along with the converging speed and the rainfall intensity, and the real hydrologic response is simulated. The physical significance of the activation function is that when f _act approaches 1, the neuron is in a strong confluence state (the grid flow production capacity is obviously improved), when f _act approaches-1, the neuron is in a weak confluence state (mainly storing water), and the threshold response characteristic of rainfall-runoff is simulated through nonlinear transformation.

Step S4.5, constructing a loss function as follows:

Where Loss represents a Loss function and N ₃ is the total number of time steps in hydrologic simulation; representing the actual flow value observed at the ith time point in m ³/s; the flow value predicted (simulated) by the model at the ith time point is expressed in m ³/s.

The following is an example of application of the present invention:

The Tunxi river basin is selected as a research area, and the Tunxi river basin is used as an important component of the Xinan river basin and is positioned in the southeast part of the southern Anhui mountain area. The main area of the mountain of the yellow mountain at the river basin is adjacent to the Yangtze river system on the west and north sides, and the mountain of the sky and the mountain of the white mountain are adjacent to the southeast. A plurality of mountain basins and valley lands with considerable scale are distributed in the flow field, the whole body of the river basin presents the characteristics of high west, low east west and low east, the relief of the topography is obvious, and the gradient of the river channel is generally larger. The method belongs to a typical subtropical monsoon climate zone in climate type, annual precipitation is kept at a higher level, and the annual average temperature is in a transition zone from a temperature zone to a subtropical zone. The vegetation coverage condition in the river basin is good, and the river basin mainly comprises various forest types and agricultural lands, so that a complex ecological system combination is formed. The elevation change range is wide, and the elevation change range from the western mountain area to the eastern valley area shows obvious decreasing trend, so that obvious terrain drop is formed.

In the implementation process, basic data acquisition work is completed firstly, and multidimensional information such as hydrological weather, topography, socioeconomic and the like is covered. The core data set comprises real-time hydrologic monitoring records, soil physical characteristic parameters, precipitation observation data of different time scales, high-precision terrain elevation models, demographic data, historical flood characteristic values and the like. The data are divided into a conventional processing set and a to-be-processed set according to processing requirements, and the conventional processing set and the to-be-processed set correspond to the stable hydrologic element and the dynamic change element respectively. The method comprises the steps of data collection, DEM data, 30m resolution SRTM DEM, coverage of a river basin and a 5km buffer area (matrix D epsilon R6150X 5200), meteorological data, weather P, air temperature T and wind speed W of 12 national weather stations at the time of 6 years, hydrological data, time-by-time flow Q of Tunxi stations, geological data, soil map (1:10 ten thousand), sandy soil 32%, loam 45%, clay 23% and 25 measured permeability coefficient Ks of drilling holes, remote sensing data, sentinel-2 images (2020-2022 years 6-8 months, cloud amount < 10%), 62 LAI real measurement points (35 forest lands, 18 farmlands and 9 grasslands). The investigation region basin DEM is schematically shown with reference to fig. 2.

And then carrying out data standardization pretreatment, and uniformly converting all indexes into the same dimension system. On the basis, a combined weighting method is adopted to determine the weight coefficient of each influence factor, namely, firstly, a judgment matrix is constructed through an expert decision model, subjective weight distribution is obtained through calculation, and meanwhile, the discrete characteristics of data are analyzed based on an information entropy theory, and an objective weight distribution scheme is derived. The two weight results are weighted and fused to form a comprehensive weight system, and a quantization basis is provided for subsequent analysis. The double-weight mechanism integrates the knowledge experience of the field, fully honors the objective rule of the data and ensures the scientificity and reliability of an evaluation system.

In the stage of space discretization processing, dynamically adjusting a grid division scheme according to the terrain complexity, and simultaneously introducing a flow direction constraint mechanism to optimize the hydrologic response unit boundary. The generated boundary is smoothed by morphological operation and curve fitting technology, so that the spatial continuity of the hydrologic network is remarkably improved. Aiming at extreme events such as heavy rain, the parameter sensitivity is adjusted in real time by adopting a dynamic activation function, so that the model can adapt to the changes of different hydrologic conditions. The whole process finally builds a multi-stage nested hydrologic discretization system, and the three-dimensional structure framework with distinct layers is formed from sub-watershed division, hydrologic response unit definition and refined grid calculation.

The specific calculation is as follows:

Step S1, arcGIS processing to obtain a first class classification basin, and processing by using ArcGIS hydrologic analysis toolkit:

1. Filling the depressions, namely repairing M concave points in the DEM;

2. The flow direction analysis, namely generating a water flow direction matrix by a D8 algorithm;

3. Sub-basin partitioning based on a confluence cumulative amount threshold value F grid;

U={U1,U2,...,U14}

And outputting N natural sub-watershed (area 82-168km ²).

Step S2:

S2.1, taking a typical sub-drainage basin U _i as an example, constructing a space matrix containing 5-dimensional topographic features. The method comprises the steps of calculating gradient S by adopting a third-order inverse distance square weight algorithm, accurately representing the steep degree of the ground surface, solving curvature Curvature through a second derivative, describing the concave-convex characteristics of the terrain, and deducing a terrain humidity index TWI based on an upstream catchment area As and a gradient angle. Finally forming a topographic feature matrix:

the gradient adopts a third-order inverse distance square weighting algorithm:

Curvature:

Terrain Wetness Index (TWI):

Where A _s is the upstream catchment area and β is the slope angle. The matrix integrates class 5 topographical attributes of the F grids, providing a data base for cluster analysis.

S2.2, fusing a dynamic grid mechanism and hydrologic flow direction constraint:

Grid self-adaptation, namely dynamically adjusting the clustering resolution according to the standard deviation sigma _Z =a of the elevation, encrypting the grid by a flat area (a < 10), expanding the grid by a steep area (a > 30), and calculating to obtain the dynamic grid size:

(original 30m→adjusted bm)

Flow direction constraint, constructing a mixed distance function:

Wherein the flow direction weight λ=0.4× (1-e ^-0.1×F) varies dynamically with the amount of confluence accumulation:

slope area (F < 5): lambda _slope = c = 0.4 (1-e-0.1 x 3) ≡d)

River region (F > 100) λ _river =e=0.4 (1-e-0.1×120) ≡f

Wherein c is a slope area flow direction weight expression, e is a river area flow direction weight expression, d is a slope area flow direction weight value, f is a river area flow direction weight value;

The design ensures that the mountain area generates the equal-height banded HRU and the river area forms the linear HRU.

S2.3 topography prototype Generation

The cluster outputs K terrain prototypes P _k＝(E_k,S_k,A_k,C_k,T_k), for example p1= (e 1, s1, a1, c1, t 1) characterizing the average elevation e1, the integrated slope s1, the dominant slope a1, the characteristic curvature c1 and the TWI mean t1 of a certain class of terrain homogenizing units, which prototypes constitute the reference points of the HRU partition.

Outputting K terrain prototypes, for example:

P5= (g, h, i, j, k) = (1,0248 m,31.2 °, southeast, 0.045,9.8)

P12= (l, m, n, o, P) = (632 m,8.7 °, northeast, -0.018,12.3)

S2.4-S2.7, HRU space generation:

1. The hydrologic constraint Voronoi segmentation, namely defining a hydrologic distance function d _hydro＝α·d_Euclidean+(1-α)·d_flow, and integrating a geographic distance d _Euclidean and a water flow path distance d _flow;

2. Topology optimization, namely reconstructing a boundary by adopting a river channel attractor algorithm by taking a river channel network gamma _channel as a hard constraint;

3. microcell fusion performing downstream fusion on HRU with area less than 0.01km2 Preferably incorporating a topographically continuous downstream unit. Fusing D HRUs with A <0.01km ²;

4. Boundary optimization:

morphological treatment: elimination of saw tooth boundaries qkm;

Bezier curve fitting, B _smooth＝BezierFit(B_final, 0.5 delta) and final outputting the hydrologic response unit set { HRU ₁,HRU₂,…,HRU_K }, wherein the first-order drainage basin and hydrologic response unit schematic diagram is shown in FIG. 3.

S2.8, taking HRU_5 as an example, the hydrologic attribute library is shown in the following table 4:

table 4. Hydrologic attribute library (HRU_5 example)

Step S3:

ANUSPLIN interpolation parameters 3 rd order spline, tension coefficient 0.25;

And (3) precision verification:

2. geological interpolation:

spatial distribution of soil permeability coefficient:

lithology classification accuracy, z% (> 80%);

3. land utilization classification, namely random forest classification, namely 500 decision trees, OOB error=an;

the main types are forest land AO%, farmland AP, water area AQ and construction land AR;

4. vegetation parameters:

LAI correlation: R ² =as (> 0.75);

Step S4:

s4.1, neuronal partitioning:

Grid size 100×100m (adaptive tuning);

total number of neurons, N ₂ = c (> 5000);

Grid size AU×AV m;

Coordinates X _i＝(x_i,y_i) (UTM projection);

S4.2 neuron properties:

static attribute matrix a _i (normalized to [0,1 ]):

Ai= [ P-, T-, soil coding, land use index, NDVI ̄, S, A ]

S4.3, establishing connection based on a D8 flow direction, and adopting a D-definition algorithm in a land leveling area (gradient <0.5 °);

s4.4 dynamic activation function heavy rain scenario example (2022-07-23):

when rt=emm/h, f _act →1 (strong stream);

When rt=fmm/h, f _act → -1 (water storage dominant);

s4.5, verifying the precision:

RMSE=1T∑t=1T(Q_obst-Q_simt)2＝gm³/s;

The conceptual diagram of the dynamic neural network constructed in step S4 refers to fig. 4.

The embodiments of the present invention have been described in detail with reference to the drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge of those skilled in the art without departing from the spirit of the present invention.

Claims

1. The dynamic DEM spatial interpolation method is characterized by comprising the following steps S1-S4, namely finishing the spatial interpolation of DEM data and the hydrologic simulation of a target river basin:

2. A method of dynamic DEM spatial interpolation according to claim 1, wherein step S1 comprises the specific steps of:

step S1.1, filling processing to eliminate DEM original data pits;

s1.2, carrying out flow direction analysis according to a D8 algorithm to generate a D8 water flow direction matrix;

S1.3, calculating the total flow accumulation amount to generate a runoff path network;

S1.4, automatically extracting the boundary of a natural sub-drainage basin unit based on a minimum water collecting area threshold method;

3. A method of dynamic DEM spatial interpolation according to claim 1, wherein step S2 comprises the specific steps of:

Wherein T _i represents the topographic feature matrix of the ith natural sub-basin unit, m represents the number of units in the research area where the topographic feature is calculated or sampled;

step S2.2, adopting an improved self-organizing map network clustering algorithm to perform dynamic grid division on natural sub-drainage basin units, calculating an elevation standard deviation sigma _Z, if sigma _Z is smaller than 10, enabling the natural sub-drainage basin units to be flat areas and reducing grid sizes, if sigma _Z is smaller than or equal to 10 and smaller than 30, enabling the natural sub-drainage basin units to be medium complexity areas, enabling the grid sizes to perform linear interpolation between minimum values and maximum values according to the numerical value of sigma _Z, and if sigma _Z is larger than 30, enabling the natural sub-drainage basin units to be steep areas and enlarging the grid sizes, wherein the specific formula is as follows:

Wherein D' represents a flow direction constraint function, ω _k is a weight of a kth topographic feature, ω ₁ is Gao Chengquan weight, ω ₂ is a slope weight, ω ₃ is a slope weight, ω ₄ is a curvature weight, ω ₅ is a topographic humidity index weight, X _k is a value of a current grid on the kth topographic feature, C _k represents a kth cluster center value, λ is a flow direction weight, θ _ij is a grid-center direction angle, and φ _flow is a water flow direction angle;

Wherein, the flow direction weight lambda is calculated as follows:

λ=0.4×(1-e^{-0.1×FlowAcc})

Wherein FlowAcc represents a confluent cumulative amount;

Step S2.3, outputting a set of clustering results P of K winning neurons based on the improved self-organizing map network output layer:

Wherein P _k represents a kth clustering result, E _k is an average elevation, S _k is a comprehensive gradient, A _k is a dominant slope, C _k is a characteristic curvature, and T _k is a terrain humidity index average value;

step S2.4, converting the clustering result into a hydrologic response unit set { HRU ₁,HRU₂,…,HRU_K } through the Voronoi diagram, wherein the hydrologic response unit has the following formula:

the hydrological distance function is as follows:

d_hydro＝α·d_Euclidean+(1+α)·d_flow

Where α is a dynamic weight coefficient between 0 and 1, D _Euclidean represents the two-dimensional geographic coordinate Euclidean distance, and D _flow represents the cumulative distance along the D8 water flow path;

S2.5, extracting a river channel network Γ _channel generated by a D8 algorithm as hard constraint, specifically using Repartitioning boundaries for the constraint pair hydrologic response units;

Wherein HRU _small represents a hydrographic response unit with an area less than threshold A _min, HRU _j represents a candidate neighboring hydrographic response unit, ΔH is the height difference, ΔL is the flow path length, and D _flow represents the water flow path distance between HRUs;

and S2.7, carrying out boundary morphology optimization on the fused hydrological response unit, and eliminating microscopic saw teeth by adopting morphology open operation, wherein the specific formula is as follows:

wherein, B _raw represents an original hydrological boundary, B _final represents a hydrological boundary after morphological optimization, and SE _3×3 represents a3×3 structural element;

B_smooth＝BezierCurveFit(B_final,tolerance=0.5δ)

Wherein, B _smooth represents a hydrologic smoothing boundary, tolerance is a tolerance threshold, and delta represents DEM data resolution;

And step S2.8, establishing a hydrologic attribute library aiming at the hydrologic response unit obtained in the step S2.7.

4. A method according to claim 3, wherein the parameters in the hydrologic attribute library established in step S2.8 include maximum confluence path length, effective flow slope, spatial soil moisture variability, solar radiation uniformity, downpenetration capacity index, drainage density, runoff coefficient.

5. A method of dynamic DEM spatial interpolation according to claim 1, wherein the specific steps of step S3 are as follows:

step S3.1, collecting time-by-time observation data of a target river basin and weather stations in a peripheral preset range, and simultaneously introducing DEM data as covariates to eliminate a terrain shielding effect;

S3.2, interpolating the soil permeability coefficient by adopting a cooperative Kriging method, using an index model for sandy soil and a Gaussian model for clay soil, and outputting geological data;

S3.3, carrying out random forest classification based on the Sentinel-2 satellite images to obtain land utilization type data;

And S3.4, calculating an NDVI index based on the Sentinel-2 image as a vegetation parameter, and removing noise through BRDF correction and time sequence harmonic analysis.

6. A method of dynamic DEM spatial interpolation according to claim 1, wherein the specific steps of step S4 are as follows:

step S4.1, dividing a target river basin into regular grids, regarding each grid unit as a neuron, and defining the neuron coordinates as the following formula by adopting a space matching algorithm of a hydrological response unit and the grids:

X_i＝(x_i,y_i)∈R²,i＝1,2,…,N₂

Wherein R represents a real number set, N ₂ is the total number of grid cells after the whole target river basin is divided into a regular grid, X _i is the coordinates of the ith neuron, and (X _i,y_i) is the coordinates of the ith neuron in the X, Y direction in the uniform projection coordinate system in step S3;

wherein θ _i represents the ith neuron, h _i, K _i represents the dynamic state of the ith neuron, which is sequentially the water depth, the hydraulic gradient and the water conductivity coefficient, A _i represents the static attribute matrix of the ith neuron and comprises the multi-source data obtained in the step S3;

S4.3, determining a connection topology based on the D8 water flow direction, and directing the water flow direction to a downstream grid with the largest gradient by calculating the gradient of each grid and 8 neighborhood grids, wherein aiming at a flat area with the gradient less than 0.5 DEG, the D-infinity algorithm is adopted to assist in judging the flow direction;

step S4.4, constructing a dynamic activation function as follows:

Wherein f _act represents a dynamic activation function, beta is a hydrologic conduction coefficient, rt is a real-time rainfall intensity, Indicating the activation intensity with the convergence speed;

step S4.5, constructing a loss function as follows:

Where Loss represents a Loss function and N ₃ is the total number of time steps in hydrologic simulation; representing the actual flow value observed at the ith time point in m ³/s; The flow value predicted by the model at the ith time point is expressed as m ³/s.