DEM adaptive multi-algorithm fusion promotion method for high-precision distributed hydrological model

Abstract

The application discloses a DEM adaptive multi-algorithm fusion promotion method for a high-precision distributed hydrological model, and relates to the field of hydrological model basic data construction.The method comprises the following steps: obtaining original DEM data, land use data and soil type data; performing normalization processing on the original DEM data to obtain unified grid DEM; calculating a basic terrain field based on the unified grid DEM; generating at least one enhanced terrain parameter (such as an enhanced slope grid or a terrain index) through a predetermined optimization pipeline based on the basic terrain field; and finally, combining the at least one enhanced terrain parameter, the land use data and the soil type data to divide a hydrological response unit.The method solves the problem of insufficient geomorphological feature reservation of traditional filtering through the optimization pipeline, and improves the physical rationality of hydrological response unit division.

Description

Technical Field

[0001] This invention relates to DEM data processing, and in particular to an adaptive multi-algorithm fusion and enhancement method for high-precision distributed hydrological models. Background Technology

[0002] Digital elevation models (DEMs) are the core foundational data for constructing high-precision distributed hydrological models. The quality and spatial continuity of topographic parameters derived from DEMs, such as slope and topographic indices, directly determine the accuracy of the model's simulation of runoff generation and confluence processes. Therefore, improving the quality of these topographic parameters is of significant research importance for flood early warning, water resource management, and eco-hydrological analysis.

[0003] In existing technologies, to eliminate noise introduced during the acquisition and processing of raw DEM data, researchers typically apply a range of image processing techniques before or after calculating terrain parameters (such as slope and topographic indices). These techniques include common mean filtering and Gaussian filtering, as well as relatively complex morphological filtering (such as opening and closing operations) and multi-scale wavelet transforms (such as wavelet soft thresholding denoising) to smooth raster data and suppress outliers.

[0004] Existing technologies also face several challenges, including: first, the blindness and isotropy of the filtering process regarding hydrogeomorphic features; and second, the static and open-loop nature of parameter settings. Specifically, conventional filters cannot distinguish between random noise and real, detailed topography (such as shallow gullies), and their isotropic smoothing properties can cross and blur the edges of gullies, disrupting the anisotropy of the water flow path. Furthermore, the setting of hyperparameters such as filter kernel size and threshold coefficients relies on human experience and is globally fixed, lacking a closed-loop adaptive mechanism to adjust upstream parameters in reverse based on downstream application effects (such as the fragmentation of hydrological response units, HRUs). Summary of the Invention

[0005] The purpose of this invention is to address the technical problems existing in the prior art by providing a DEM adaptive multi-algorithm fusion enhancement method for high-precision distributed hydrological models.

[0006] The technical solution, a DEM adaptive multi-algorithm fusion and enhancement method for high-precision distributed hydrological models, includes the following steps:

[0007] Step S1: Collect multi-source DEM data of the target watershed, as well as land use and soil attribute data;

[0008] Step S2: Conduct water system matching assessment and data source screening for multi-source DEM data;

[0009] Step S3: Fill the depression DEM, calculate the slope grid, and use morphological denoising and wavelet denoising; combine the slope zoning data before and after denoising, as well as land use and soil attribute data, to divide the hydrological response units.

[0010] Step S4: Calculate the upslope contribution area and other preliminary quantities based on the river network data to obtain the topographic index grid. Combine the slope classification results to adjust the partition weight of the topographic index grid for noise suppression based on the difference in hydrological response between gentle and steep areas.

[0011] An adaptive multi-algorithm fusion and enhancement method for DEM (Distributed Hydrological Model) for high-precision distributed hydrological models includes:

[0012] Obtain raw DEM data, land use data, and soil type data;

[0013] The original DEM data is normalized to obtain a unified grid DEM;

[0014] Calculate the basic topographic field based on a unified grid DEM;

[0015] Based on the base terrain field, at least one enhanced terrain parameter, including a slope raster and a terrain index, is generated through a predefined optimized pipeline.

[0016] Hydrological response units are delineated by combining at least one enhanced topographic parameter, land use data, and soil type data.

[0017] According to one aspect of this application, the method further includes:

[0018] Obtain the original terrain parameters from the basic terrain field;

[0019] Visualize and compare the original and enhanced terrain parameters and perform raster difference analysis.

[0020] Based on the divided hydrological response units, HRU morphological indices are calculated to quantitatively assess the improvement in fragmentation of the hydrological response units.

[0021] According to one aspect of this application, the steps for calculating HRU morphological indices include:

[0022] Extract the area and perimeter of each sub-polygon in the hydrological response unit;

[0023] Based on area and perimeter, the HRU morphological index is calculated from at least one of the following: shape index; compactness; fractal dimension; area-to-perimeter ratio.

[0024] Beneficial effects: This invention solves the problem of insufficient preservation of landform features by traditional filtering through this optimized pipeline, improves the physical rationality of hydrological response unit division, and effectively reduces HRU fragmentation. Attached Figure Description

[0025] Figure 1 This is a flowchart of a DEM adaptive multi-algorithm fusion improvement method for high-precision distributed hydrological models according to the present invention.

[0026] Figure 2 Comparative analysis chart of multi-source DEM data.

[0027] Figure 3 Image showing the results of water system extraction using multi-source DEM.

[0028] Figure 4 The results of D∞ slope tangent for different noise reduction methods.

[0029] Figure 5 Comparison of noise reduction results for mixed terrain indices.

[0030] Figure 6 Comparison of HRU partitioning results before and after noise reduction. Detailed Implementation

[0031] Example 1 provides an adaptive multi-algorithm fusion enhancement method for high-precision distributed hydrological models (DEMs). It details the complete technical process from multi-source data preprocessing and basic field construction to generating enhanced topographic parameters and dividing hydrological response units.

[0032] Step 101: Obtain the original DEM data, land use data, and soil type data.

[0033] Raw DEM data refers to digital elevation model raster data obtained directly from the data source without processing by the method of this invention. To address the interference caused by quality differences between different data sources on subsequent hydrological analysis, this invention requires the acquisition of multi-source DEM data. Exemplarily, these data sources may include, but are not limited to: hydrologically optimized versions (con) and unoptimized raw versions (dem) of HydroSHEDSDEM, ASTERGDEM, NASADEM, FABDEM optimized for vegetation impact, and high-precision dual-benchmark TanDEM. Land use data and soil type data are essential underlying surface auxiliary data for delineating hydrological response units (HRUs). For example, data from the Chinese Academy of Sciences Land Use Data (CNLULC) and the World Soil Database (HWSD) can be used.

[0034] Preferably, after acquiring the original DEM data, this step further includes multi-source DEM water system matching degree assessment and data source screening. This assessment is used to select the data that best reflects the water system morphology of the study area from numerous data sources, which is a key prerequisite for ensuring the accuracy of all subsequent calculations. A specific assessment method is as follows: using a unified tool (such as TauDEM) to extract the water system for each DEM data source to be assessed, this process typically includes depression filling, flow direction calculation, cumulative flow calculation, and threshold-based river network extraction; the river network data extracted from each DEM is visualized and quantitatively compared with recognized standard water system data; considering the water system matching degree, data volume (computational cost), and actual water system characteristics, the optimal data source is selected (for example, in one embodiment, HydroSHEDS con is selected as the optimal data source) for subsequent steps.

[0035] Step 102: Standardize the original DEM data to obtain a unified grid DEM.

[0036] Normalization is a necessary step to resolve spatial inconsistencies in multi-source data (including DEM, CNLULC, and HWSD). Its purpose is to construct a spatially aligned, resolution-consistent analytical grid, i.e., a unified grid DEM. Specifically, this process includes at least: coordinate system one, which transforms all data (DEM, soil, and land use data) to the same projected coordinate system; resolution resampling and raster alignment, which uses the selected optimal DEM (such as the one chosen in step 101) or a certain reference data (such as NASADEM) as a reference, and unifies the remaining data to the same spatial resolution (such as 30m or 90m) through resampling (such as bilinear interpolation, nearest neighbor method, etc.), ensuring that the pixel centers or corners of all rasters are strictly aligned spatially; in addition, it may include data cropping, outlier removal, and invalid value (NoData) filling operations.

[0037] Step 103: Calculate the basic topographic field based on the unified grid DEM.

[0038] The basic topographic field refers to a series of basic parameter rasters derived from a unified grid DEM for hydrological analysis. In this invention, the basic topographic field includes at least: a flow direction field (e.g., D8 or D∞ flow direction), a cumulative runoff raster (or a specific runoff area raster α derived therefrom), a slope raster (e.g., slope angle β or slope tangent tanβ), a topographic curvature raster, and channel masks (i.e., river network locations) and watershed labels for constraints. These basic topographic fields are inputs for subsequent generation of enhanced topographic parameters. In some embodiments, the basic topographic field includes at least a true slope raster and a flow direction raster.

[0039] Step 104: Based on the base terrain field, generate at least one enhanced terrain parameter, including an enhanced slope grid and an enhanced terrain index, through a predetermined optimization pipeline.

[0040] In this step, the enhancement refers to the fact that, compared to the base slope or topographic index calculated in step 103, the parameter raster generated in this step is more robust to noise (such as stripes, artifacts, and local distortions) in the original DEM, and can more realistically reflect hydrogeomorphic features. The predetermined optimization pipeline refers to a series of algorithm combinations provided by this invention to achieve this enhancement objective. This invention provides at least two specific ways to implement this pipeline, which will be detailed in subsequent embodiments: one way is to perform morphological filtering and wavelet transform denoising sequentially on the base slope raster to obtain an enhanced slope raster; the other way is to perform partition-adaptive morphological filtering and conservative wavelet denoising on the base topographic index, and perform weighted fusion of multiple results to obtain an enhanced topographic index. More preferably, there is also a more advanced optimization pipeline, which is implemented through a series of adaptive steps such as topological persistence-gated PDE evolution, pixel-level variable weight fusion, and graph variational solution.

[0041] Step 105: Combine at least one enhanced topographic parameter, land use data, and soil type data to delineate hydrological response units.

[0042] Hydrological Response Units (HRUs) are the smallest spatial units with similar hydrological response characteristics in a distributed hydrological model. The purpose of this step is to use the enhanced terrain parameters generated in step 104 to replace the basic parameters to guide HRU partitioning, thereby avoiding the problem of excessive HRU fragmentation caused by noise.

[0043] In the specific division process: the enhanced slope raster obtained in step 104 is classified according to predetermined rules (such as the quantile method) (for example, divided into two intervals: 0%-10% and >10%); the classified slope raster is spatially overlaid with the normalized land use raster and soil type raster; continuous pixel areas with the same slope level, the same land use type, and the same soil properties are merged and defined as a single HRU. Due to the use of enhanced terrain parameters, the resulting HRUs are more spatially continuous, more regularly shaped, and have a clearer physical meaning.

[0044] Example 2 describes an enhanced slope raster generation method. It is a specific implementation of step 104 in Example 1. Specifically, it describes how an enhanced slope raster is generated using a pipeline combining morphological filtering and wavelet transform denoising.

[0045] Step 201: Obtain the slope grid from the basic topographic field.

[0046] This slope raster is one of the base topographic fields calculated in step 103 of Example 1. Preferably, in order to preserve the subtle undulations of the original terrain and avoid excessive smoothing of the terrain by traditional depression-filling preprocessing, this step uses the D∞ algorithm to process the unified grid DEM (here referring to the unfilled DEM) obtained in step 102 of Example 1 to obtain the terrain slope raster. The D∞ algorithm determines the flow direction and slope by calculating the maximum slope direction from each cell to its 8 neighboring cells. The output slope raster can be stored in the form of tangent (tanβ) or angle (β). This raster (denoted as Slope) raw While preserving true details, it also contains a significant amount of noise (such as stripe noise) introduced by DEM acquisition or interpolation, which requires processing in subsequent steps.

[0047] Step 202: Perform morphological filtering on the slope raster to generate a morphologically denoised slope raster.

[0048] Morphological filtering is used to eliminate the original slope raster. raw High-frequency small-scale noise in the slope raster. Specifically, this step may include opening and closing operations. First, a structuring element (filter kernel) is specified, such as a 3x3 rectangular structuring element. Then, an opening operation (i.e., erosion followed by dilation) is performed to remove small bump noise (such as isolated anomalous high slope points) in the slope raster. Next, a closing operation (i.e., dilation followed by erosion) is performed to fill small depression noise (such as local low slope depressions caused by data acquisition errors) in the slope raster. During processing, invalid value areas (such as water areas) can be left unchanged, and only the effective morphological areas are corrected. The morphologically denoised slope raster generated in this step (denoted as Slope) morph ).

[0049] Step 203: Apply wavelet transform to the morphologically denoised slope raster to obtain an enhanced slope raster.

[0050] This step aims to improve the morphologically denoised slope raster output from step 202. morph Further fine-tuning denoising is performed on the slope to preserve key terrain details that might have been smoothed during morphological filtering. Specifically, Discrete Wavelet Transform (DWT) is used to refine the denoising. morphThe process involves: selecting a suitable wavelet basis (e.g., the Daubechies wavelet basis, such as db4) and a decomposition level (e.g., 3 levels); performing multi-level decomposition on the slope data, separating it into approximation coefficients (low-frequency components representing the overall topographic trend) and detail coefficients (high-frequency components containing residual noise and topographic details); applying a threshold (e.g., soft thresholding) to the detail coefficients to suppress noise while preserving coefficients characterizing key topographic details such as valley edges and terraces; and finally, reconstructing the signal using the processed approximation and detail coefficients through inverse discrete wavelet transform (IDWT) to obtain the final enhanced slope raster. final ).

[0051] In some alternative implementations, the processing objects of steps 202 and 203 can differ. For example, steps 202 (morphological denoising) and 203 (wavelet transform denoising) can both use the original slope raster obtained in step 201. raw The original slope raster, morphological denoising result, and wavelet denoising result can be fused (e.g., weighted average) or selected as the final enhanced slope raster based on the evaluation results of downstream applications (e.g., HRU partitioning) (e.g., the fragmentation index in Example 4).

[0052] Example 3 describes an enhanced terrain index generation method, as another specific implementation of step 104 in Example 1. It specifically explains how to generate an enhanced terrain index through a pipeline of partitioned adaptive filtering, conservative wavelet denoising, and multi-path weighted fusion.

[0053] Step 301: Calculate the initial topographic index based on the basic topographic field.

[0054] The initial topographic index (ATB) is a key parameter characterizing a watershed's runoff generation capacity and catchment trends. In this embodiment, this initial topographic index (denoted as ATB) raw ATB is calculated using the following formula: raw =ln(α / tanβ) or ATB raw =ln(α / β). Where α is the specific catchment area grid; tanβ is the slope tangent grid. A very small positive number (e.g., 1e-6) is added when calculating the quotient to prevent division by zero. Both α and tanβ can be obtained from the base topographic field in step 103 of Example 1. The α grid reflects the total water volume collected upstream of each cell, while tanβ reflects the gradient of the local flow. This ATB... rawAlthough the grid reflects the hydrological characteristics macroscopically, the original DEM often contains unnatural structures, such as diagonal artifacts, due to noise transmission. These artifacts can interfere with subsequent hydrological simulations, so enhancement processing is required.

[0055] Step 302: Based on the slope information in the basic topographic field, establish topographic zones, which include gentle areas and steep areas.

[0056] To identify the steepness of the terrain so that filters of varying intensities can be applied in subsequent steps to achieve adaptive zoning, this embodiment is implemented as follows:

[0057] Specifically, we can first calculate the mean terrain gradient G(i,j) for each pixel (i,j). mean For example, the central difference method can be used to calculate the gradient in the x-direction G_x=(Z(i,j+1)-Z(i,j-1)) / (2d) and the gradient in the y-direction G_y=(Z(i+1,j)-Z(i-1,j)) / (2d), where Z is the elevation and d is the grid spacing. Then, the mean or magnitude of the two values ​​is taken as G(i,j). mean Set a gradient threshold T. split (For example, T) split =0.5), dividing all pixels into two sets: when G(i,j) = 0.5. mean Greater than T split When G(i,j) is in this state, the cell is classified as a steep region; when G(i,j) is in this state... mean Less than or equal to T split At that time, the pixel was classified as a flat area. Steep areas usually correspond to key terrain features such as valleys and ridges, which require careful protection; flat areas are more sensitive to noise and require stronger smoothing processing.

[0058] Step 303: Perform zonal adaptive morphological filtering on the initial terrain index, where a first-size filter kernel is used in flat areas and a second-size filter kernel is used in steep areas. The first size is larger than the second size, thereby generating an adaptively morphologically filtered terrain index.

[0059] This step utilizes the terrain zoning (flat and steep areas) generated in step 302 to refine the initial terrain index ATB generated in step 301. rawDifferentiated processing is performed. Specifically, for pixels in flat areas, a first-size filter kernel (e.g., a 5x5 or 7x7 rectangular structuring element) is used to perform morphological opening and closing operations. Due to the larger kernel size (first size), it can effectively smooth out minute-scale depressions or bumps in the region, preserving the macroscopic water flow trend. For pixels in steep areas, a second-size filter kernel (e.g., a 3x3 rectangular structuring element) is used to perform morphological opening and closing operations. Due to the smaller kernel size (second size), this operation only removes extremely small bumps, minimizing the loss of detail in steep terrain (such as valley edges). The processing results for the two regions are combined to obtain the adaptive morphological filtering topographic index (denoted as ATB). morph ).

[0060] Step 304: Perform conservative wavelet denoising on the initial terrain index to generate the wavelet-denoised terrain index.

[0061] This step is a parallel processing branch to step 303, also using the initial terrain index ATB. raw The input is [value]. Conservative wavelet denoising aims to suppress only high-frequency random noise while preserving, to the maximum extent possible, low-frequency approximation coefficients that characterize the overall topographic trend, avoiding excessive smoothing that blurs valley edges. In practice, discrete wavelet transform (such as the db4 wavelet basis) can be used, but a lower number of decomposition levels (e.g., 2 levels) and smaller soft thresholding coefficients are chosen. For example, the soft threshold T can be calculated using T = fs * σ. d The expression *sqrt(2*log(N)) is used to calculate the value of N, where N is the total number of pixels and σ is the value of N. d is the standard deviation of the detail coefficients, and fs is the threshold coefficient. To achieve conservative noise reduction, fs can be taken as a small value (e.g., 0.5 or 0.7) to set a lower threshold, suppressing only the most significant noise. The reconstructed topographic index (denoted as ATB) is obtained after wavelet denoising. wave ).

[0062] Step 305 is a common intermediate step applicable to various denoising methods such as morphology or wavelet denoising. A buffer mask and normalized distance weight map are generated based on line features for spatially adaptive denoising intensity adjustment. The line feature vector data with the same coordinate system as the input raster data is acquired, and the line feature is buffered and expanded by a preset buffer distance to form a buffer geometry covering the target area. Using raster affine transformation parameters, the buffer geometry is converted into a binary mask of the same size as the input raster, where pixels with a true mask value are located within the buffer. Simultaneously, a very small distance buffer (e.g., 1% of the buffer distance) is applied to the original line feature to generate a thin line mask that approximates the center position of the line, and the Euclidean distance from each raster cell to the nearest line feature is calculated based on this mask. Within the buffer range, the Euclidean distance is divided by the buffer distance for normalization processing to obtain a normalized distance map (normalized distance ∈ [0,1]), where 0 corresponds to the center position of the line feature and 1 corresponds to the buffer boundary. Further, a spatially varying denoising intensity weight is constructed based on the normalized distance map: a lower denoising intensity (high fidelity weight) is assigned near the center of the line feature, and a higher denoising intensity (strong smoothing weight) is assigned at the edge of the buffer, thereby achieving spatial adaptive control of the denoising process.

[0063] Step 306: The initial terrain index, the terrain index after adaptive morphological filtering, and the terrain index after wavelet denoising are weighted and fused to obtain the enhanced terrain index.

[0064] This step combines the advantages of multiple processing branches, balancing noise suppression and detail preservation. Specifically, a weighted average is used to calculate the ATB from step 301. raw Step 303 ATB morph and ATB in step 304 wave Perform the fusion. The fusion formula is: ATB final =W0*ATB raw +W1*ATB morph +W2*ATB wave Where W0, W1, and W2 are the weights of the three, and satisfy W0 + W1 + W2 = 1. W0 retains some original information, W1 (ATB) morph W2 (ATB) contributes to the partition-adaptive morphological smoothing properties. wave This contributes to the conservative wavelet denoising characteristics. By reasonably setting the weights (e.g., W0=0.2, W1=0.4, W2=0.4), the final output ATB... final (i.e., enhanced topographic index) both eliminates the original noise (such as diagonal artifacts) and preserves key hydro-topographic features.

[0065] Example 4 describes a method for evaluating and verifying noise reduction effectiveness. This evaluation is crucial for verifying whether the enhanced terrain parameters are truly superior to the original terrain parameters.

[0066] Step 401: Obtain the original terrain parameters from the basic terrain field; perform a visual comparison and raster difference analysis on the original terrain parameters and the enhanced terrain parameters.

[0067] The original terrain parameters refer to the Slope in step 201 of Example 2. raw Or the ATB in step 301 of Example 3 raw The enhanced terrain parameter refers to the Slope parameter in step 203 of Example 2. final Or the ATB in step 305 of Example 3 final .

[0068] Visual comparison refers to displaying raster images before and after noise reduction side-by-side. For example, when comparing slope raster images, the original raster (Slope) is displayed alongside the original raster image. raw It may show obvious striped noise or scan line artifacts; while in the enhanced grid (Slope) final In this context, these stripes should be effectively suppressed, while topographical details such as the edges of gullies should remain clearly visible.

[0069] Raster difference analysis refers to calculating the difference between rasters before and after noise reduction (e.g., Diff = Slope). raw -Slope final The difference raster (Diff) is then rendered. This difference map visually displays the spatial location and intensity of noise suppression, helping to quantify the spatial characteristics of noise suppression.

[0070] Step 402: Based on the divided hydrological response units, calculate the HRU morphological index to quantitatively assess the improvement in the fragmentation of the hydrological response units.

[0071] To achieve this step, two HRU partitioning operations need to be performed (as in step 105 of Example 1): the first is the baseline partitioning, using the original terrain parameters (such as Slope). raw By combining land use and soil data, a set of original HRU vector datasets (denoted as HRU) is generated. raw The second step is enhanced partitioning, using enhanced terrain parameters (such as slope). final By combining the same land use and soil data, an enhanced HRU vector dataset (denoted as HRU) is generated. final ).

[0072] Quantitatively assessing the improvement in fragmentation of hydrological response units refers to comparing the HRU... raw and HRU finalThe morphological parameters of the HRU were used to verify whether noise reduction alleviated the fragmentation within the HRU. If noise in the original parameters caused the HRU to be fragmented into many small, irregular pieces, then the enhanced HRU... final It should be characterized by: a reduction in the number of sub-polygons, more regular shapes, and fuller forms.

[0073] Step 403: Extract the area and perimeter of each sub-polygon in the hydrological response unit; based on the area and perimeter, calculate an index selected from at least one of the following to form the HRU morphological index.

[0074] This step is detailed in defining the HRU morphological indicators in step 402. For HRU raw and HRU final For each sub-polygon in the dataset, its area S and perimeter P are extracted. Based on S and P, at least one of the following metrics is calculated:

[0075] Shape Index (IS): For example, IS = 0.25 * P / sqrt(S). The larger this value, the more complex the shape or the further it deviates from a circle. After noise reduction, the mean of IS is expected to decrease.

[0076] Compactness (IC): For example, IC = (4*π*S) / (P*P). The closer this value is to 1, the more compact the shape (closer to a circle). After noise reduction, it is expected that the mean of IC will increase (or the variance will decrease, indicating a more consistent shape).

[0077] Fractal dimension (ID): For example, ID = (2 * log(P)) / log(S). The closer this value is to 1, the smoother the boundary. After noise reduction, the mean of ID is expected to decrease.

[0078] Area-to-perimeter ratio (IR): For example, IR = S / P. The larger this value, the fuller the HRU is, rather than a narrow, elongated fragment. After noise reduction, the mean IR is expected to increase.

[0079] Furthermore, the number of sub-polygons n contained in each HRU can be directly counted. After denoising, the mean and variance of the expected n decrease significantly. By comparing the mean and variance of these indicators before and after denoising, the improvement effect of the method of this invention on the quality of HRU can be objectively and quantitatively verified.

[0080] Example 5: This example provides a preferred DEM post-processing workflow for implementation after previous examples, aiming to generate a topographic index product with high structural fidelity, good noise suppression, and high hydrological consistency.

[0081] Step 501: Based on the basic topographic field, generate at least two candidate topographic indices.

[0082] The purpose of this step is to generate different estimates of the terrain index and slope raster through various processing paths, each with different noise characteristics. In this embodiment, three candidates are generated:

[0083] Original Candidate (DVar) raw The slope raster is obtained through the D∞ algorithm. For the terrain index, the specific catchment area raster α and the slope tangent raster tanβ are used, calculated using the formula ATB. raw The result is calculated using ln(α / (tanβ)). This candidate retains the most original terrain information.

[0084] Morphological candidate (DVar) morph ): For DVar raw The grid performs small-scale (e.g., 3x3) morphological opening and closing operations. This candidate excels at removing local, small-scale noise, resulting in smoother local shapes.

[0085] Wavelet Candidate (DVar) wave ): For DVar raw The grid is processed using multi-level discrete wavelet soft thresholding (as in step 203 of Example 2). This candidate is adept at handling multi-scale noise, suppressing noise in specific frequency bands while preserving macroscopic trends.

[0086] Step 502: Define the quality evaluation index system required for the outer objective function. This system is used to quantify the overall quality of the fused final product DVar_u output from the inner processing pipeline, and includes at least the following three categories of indicators:

[0087] Hydrological response unit fragmentation index (M_hru): As in step 403 of Example 4, after the inner pipeline outputs DVar_u, HRU partitioning is performed, and the number of sub-polygons n, shape index IS, compactness IC, etc. of HRU are calculated to evaluate fragmentation.

[0088] River network geometry index (M_river): Based on the ATB_u output from the inner pipeline and the flow direction field and channel mask in the basic topographic field, it evaluates the generated river network geometry, such as calculating the connectivity of the river network, the deviation of the river network density from the standard value, and whether there are unreasonable branching.

[0089] Terrain Artifact Metric (M_topo): Based on the DVar_u output of the inner pipeline, it evaluates the suppression of terrain artifacts, such as the number of diagonal artifacts remaining in the quantized raster.

[0090] Step 503: Construct the outer objective function J(θ). The multiple quality indicators defined in step 502 are mapped to optimizable scalar values ​​through a scoring function, and a unified objective function is formed by weighted combination.

[0091] This step combines the multiple metrics from step 502 into a single, minimized (or maximized) scalar value, namely the outer objective function J(θ). This function can take the form of a weighted combination of multiple metric scores: J(θ) = λ1 * Score_HRU(M_hru) + λ2 * Score_Topo(M_topo); where λ1 and λ2 are the weights of each metric. Score_HRU penalizes high fragmentation, and Score_Topo penalizes artifact persistence. θ represents all the hyperparameters affecting these metrics.

[0092] Step 504: Define one or more hyperparameters (such as morphological kernel size, wavelet decomposition layer number, threshold scaling factor, fusion weight, etc.) involved in the process of generating enhanced terrain parameters as the set of parameters to be optimized θ, and introduce an optimization strategy for automatic optimization.

[0093] The process involves iterative evaluation and optimization. Specifically, the optimizer proposes a set of candidate hyperparameters θ. k The system is based on θ k Configure the inner processing pipeline to generate the fused terrain index DVar_u, use it to execute the steps of Examples 1-3, and then execute step 502 of this example (calculate M_hru, M_topo). The system calculates the outer objective function J(θ) in step 503. k This yields a scalar value. The result (θ_k, J(θ_k)) is fed back to the optimizer. The optimizer updates its surrogate model based on this new information and prepares to propose the next set of parameters θ. k+1 This loop repeats until the value of J(θ) converges (e.g., no significant improvement after multiple iterations) or the maximum number of iterations is reached. The optimal set of parameters θ* output by the optimizer when it stops is the best processing pipeline configuration automatically optimized by this invention for specific DEM data.

[0094] Example 6 describes an application case of DEM slope enhancement and HRU improvement.

[0095] The data preparation in Example 1 was performed. This case study selected a small and medium-sized watershed in a key flood-prone area in China as the study area. Multi-source DEM data were collected, including six types such as HydroSHEDScon (90m resolution), ASTERGDEM (30m), and NASADEM (30m), as well as soil data from HWSDv2.0 and LULC land use data.

[0096] Perform the multi-source DEM river system matching assessment in step 101 of Example 1. Extract the river systems from each DEM using the TauDEM tool and compare them with standard river system data. The assessment revealed that the river network extracted by HydroSHEDScon showed a high degree of agreement with the actual situation and had a smaller data volume. Therefore, HydroSHEDScon (90m) was selected as the original DEM data for this case, and the remaining data (CNLULC, HWSD) were aligned to complete steps 102 and 103 of Example 1.

[0097] Perform the slope enhancement process in Example 2. Step 201: Calculate the slope of the original HydroSHEDS conDEM using the D∞ algorithm to obtain the original slope raster. raw Visualization (Step 401 of Example 4) revealed that the slope raw The grid exhibits noticeable striped noise.

[0098] This case study compares two noise reduction approaches. Approach A (morphological noise reduction): Step 202 is executed to reduce the slope. raw Morphological filtering was performed using a 3x3 rectangular structuring element to obtain the slope. morph Path B (Wavelet Denoising): Execute step 203 (this directly affects the Slope). raw The slope is obtained by performing a three-level decomposition using the Daubechies wavelet basis (db4), followed by soft thresholding and reconstruction. wave .

[0099] Perform the quantitative assessment as described in Example 4. Step 402: Use Slope... raw Slope morph and Slope wave Combining LULC and SOIL data, three sets of HRU datasets were obtained: HRU raw HRU morph HRU wave Step 403: Calculate the morphological indices (number of polygons n, shape index IS, compactness IC, fractal dimension ID, area-to-perimeter ratio IR) for all HRUs in these three datasets, and calculate their mean and standard deviation.

[0100] The evaluation results show that:

[0101] In terms of fragmentation (number of polygons n): HRU raw The mean is 7.32. HRU morph and HRU wave The average value of both noise reduction methods decreased to 4.37 (a decrease of approximately 40.3%), indicating that both noise reduction methods significantly alleviated the fragmentation within the HRU.

[0102] In terms of morphological regularity (shape index IS): HRU raw The mean is 3.27. HRU wave The mean value dropped to 2.86 (a decrease of 12.45%), which is better than HRU. morph (2.82) indicates that the shape of the HRU becomes more regular after noise reduction.

[0103] In terms of morphological fullness (area-to-perimeter ratio IR): HRU raw The mean is 99.86. HRU wave The average score rose to 143.74 (an increase of 43.95%), significantly outperforming HRU. morph (135.23) indicates that wavelet denoising makes the HRU shape more complete and reduces narrow fragments.

[0104] In terms of consistency (standard deviation of compact IC): HRU raw The standard deviation is 0.21. HRU wave The standard deviation decreased to 0.20 (a decrease of 4.23%), while HRU morph The standard deviation actually increased slightly. This indicates that wavelet denoising improves the shape consistency of the HRU.

[0105] Example 7: A DEM adaptive multi-algorithm fusion improvement method for high-precision distributed hydrological models, specifically including the following steps:

[0106] Step S1: First, data collection and preprocessing are carried out for the small and medium-sized watershed.

[0107] Step S1.1 Data Acquisition: For the target watershed, the system collects multi-source digital elevation model (DEM) data and hydrological correlation auxiliary data. The DEM data includes the following six types: HydroSHEDS DEM (90m resolution) in both hydrologically optimized (con) and unoptimized original (dem) versions; ASTER GDEM (30m resolution) for orthographic datum; NASA DEM (30m resolution); FABDEM with optimized vegetation impact; and TanDEM (30m resolution) for high-precision dual datum. Auxiliary data includes soil attribute data (1km resolution) from the Harmonized World Soil Database (HWSDv2.0) and land use data from the Chinese Academy of Sciences (LULC, 2000 period, 30m resolution), providing hydrological-underlying surface correlation basis for subsequent topographic parameter calculations and HRU delineation.

[0108] Step S1.2 Perform quality purification on the collected raw data, crop the TanDEM data, and remove some areas with obvious low quality.

[0109] Step S1.3 Coordinate System 1 and Resampling: To achieve spatial alignment and comparability of multi-source data, all DEM, soil and land use data are uniformly converted to the projected coordinate system to eliminate projection differences between different data sources; using NASADEM as the reference data, other data are resampled, and pixel alignment in each DEM space is ensured, and spatial differences between data are eliminated.

[0110] Figure 2 As shown, each data layer is perfectly matched in coordinate system, spatial range, and resolution, providing a standardized and highly consistent foundation dataset for subsequent steps.

[0111] Step S2: For the six types of DEM data sources to be evaluated (HydroSHEDS con / dem, ASTER GDEM, NASA DEM, FABDEM, TanDEM), assess the water system matching degree and select suitable data sources.

[0112] Step S2.1 DEM water system extraction: A unified workflow based on the TauDEM toolset was used to extract the water system, ensuring the controllability of the comparative experiment. Depression filling was performed on each DEM data point, and the D8 algorithm was used to generate flow direction grids (recording the flow direction of each pixel) and cumulative flow grids (reflecting the water catchment intensity per unit area). Based on the cumulative flow grid usage and the cumulative flow threshold (equivalent to 18 km²), the water system was extracted. 2 (Number of raster cells), extract continuous cells with flow exceeding a threshold to generate initial river network vector data. River system matching assessment and data source selection.

[0113] After extracting the river networks from each DEM in step S2.2, they were compared and analyzed with the standard drainage data of the study area, as shown in Figure 3. The results show that the river networks extracted by HydroSHEDS con and FABDEM highly match the actual situation; while other products exhibited significant disconnections. Furthermore, the data volume of HydroSHEDS con is only 1 / 4 that of other 30m resolution DEMs, significantly reducing storage and computational costs while ensuring extraction accuracy. In conclusion, HydroSHEDS con was ultimately selected as the optimal data source.

[0114] Step S3: Slope data noise reduction.

[0115] Step S3.1: Calculate the true slope (D∞) based on the original DEM. Input the original data of HydroSHEDS con DEM (90m resolution) and use the D∞ algorithm in the TauDEM toolset to calculate the slope, preserving the subtle undulations of the original terrain. The D∞ algorithm generates raster data reflecting the true terrain slope (slope values ​​are stored in tangent form) by calculating the direction of maximum slope from each cell to its neighborhood.

[0116] Step S3.2: Morphological denoising of the slope raster. For the original slope raster output by the D∞ algorithm (slope values ​​are represented by tangent values, containing striped noise and local outliers), morphological filtering is used for denoising: using a 3×3 rectangular structuring element as the core, first perform an opening operation (erosion followed by dilation) to remove small bump noise (such as isolated abnormal high slope points), and then perform a closing operation (dilation followed by erosion) to fill small depression noise (such as local low slope depressions caused by data acquisition errors).

[0117] Step S3.3: Wavelet transform noise reduction process

[0118] For the original slope raster output by D∞, fine noise reduction is performed using Discrete Wavelet Transform (DWT): the Daubechies wavelet basis (db4) is used for three-level decomposition, separating the slope data into approximate coefficients (low frequency, representing the overall terrain trend) and detail coefficients (high frequency, containing noise and terrain details). Soft thresholding is applied to the detail coefficients to suppress noise, preserving coefficients corresponding to significant terrain features; subsequently, wavelet-denoised slope raster is reconstructed through inverse discrete wavelet transform, effectively preserving key terrain details such as valley edges and terraces.

[0119] Step S3.4: HRU segmentation. Using the original slope raster, morphologically denoised slope raster, and wavelet-denoised slope raster, the slope is divided into two intervals (0%-10% and >10%) according to the quantile method. Land use data from the Chinese Academy of Sciences and HWSDv2 soil data are overlaid, and continuous pixels within the same slope interval, with consistent land use types and similar soil permeability, are segmented into Hydrological Response Units (HRUs). Finally, three sets of HRU segmentation results are generated and compared in step S5.

[0120] Step S4: Process cleaning and correction:

[0121] Step S4.1: Topographic index calculation. Based on the HydroSHEDS con DEM selected in S2, a depression filling operation is first performed. Then, according to the river network generated in S2, the elevation of the river network cells is forcibly set to the lowest neighbor value downstream to avoid deviations in the slope calculation of the river network area. The slope gradient β of the cell is calculated by the ratio of the elevation difference to the distance between adjacent cells; and the total catchment area α of each cell is calculated. Finally, a topographic index raster is generated using the formula ATB=ln(α / β), which reflects the comprehensive characteristics of unit catchment area and slope (the larger the value, the stronger the runoff generation capacity).

[0122] Step S4.2: To address the diagonal unnatural noise in the topographic index raster (mostly caused by original DEM error propagation or interpolation artifacts), a multi-strategy hybrid denoising approach is adopted: First, the topographic gradient is calculated, and a gradient threshold (set to 0.5) is used to distinguish between steep and gentle regions: a 3×3 small kernel is used for opening and closing operations in steep regions; a 5×5 large kernel is used for opening and closing operations in gentle regions, completing adaptive morphological filtering. A 3-layer discrete wavelet decomposition (db4 wavelet basis) is used to separate approximation coefficients (overall topographic trend) and detail coefficients (high-frequency noise). A soft threshold (threshold is 0.7 times the standard deviation) is applied to the detail coefficients to suppress high-frequency random noise; after reconstruction, the raster size is adjusted through bilinear interpolation, focusing on preserving the low-frequency trend of the topographic index raster, completing conservative wavelet denoising. Based on the original topographic index raster (weight W0=0.2), the results of mild wavelet denoising (weight W1=0.4) and morphological filtering (weight W2=0.4) are superimposed, and noise suppression and detail preservation are balanced through weighted fusion.

[0123] Step S5: Comparative Analysis of Terrain Features Before and After Noise Reduction

[0124] Step S5.1 Visual comparison and difference analysis of terrain parameters before and after noise reduction: To intuitively evaluate the noise reduction effect, the results of slope and terrain index ATB before and after noise reduction can be visualized and compared visually; calculate the difference between the raster before and after noise reduction, and display the difference distribution through raster rendering.

[0125] Figure 4 shows a comparison between steps S3.2 and S3.3 and the original data. The original slope raster shows obvious striped noise (caused by DEM interpolation errors or scan line artifacts). The middle image shows the morphological denoising result, where the striped noise is largely eliminated, but the slope values ​​of some steep terrains (such as shallow gullies and small ravines) are excessively smoothed. The right image shows the wavelet denoising result, where the striped noise is also suppressed, and the details of steep terrains (such as valley branches) are preserved more completely (circled marked areas). The comparison shows that morphological denoising focuses on global noise suppression, while wavelet denoising is more advantageous in preserving local details. The final output terrain index raster (such as...) Figure 5 As shown in the figure, the original diagonal noise is basically eliminated, but the topographic index value along the valley is slightly lost due to minor fluctuations. Overall, it meets the requirements of the hydrological model for the tolerance of topographic index noise.

[0126] Step S5.2 Quantitative assessment of HRU fragmentation degree and verification of noise reduction effect

[0127] Based on the HRU vector data (multipolygons) generated from the original slope raster, morphologically denoised slope raster, and wavelet-denoised slope raster in S3, the change in the degree of fragmentation within the HRU after denoising was quantified by calculating the average morphological index of the HRU set. The mean results are shown in Table 1, and the variance results are shown in Table 2. The mean number of polygons in the original HRU was 7.32, which decreased to 4.37 (-40.26%) after morphological denoising and 4.37 (-40.28%) after wavelet denoising, indicating that both denoising methods significantly alleviated the fragmentation within the HRU. The mean shape index of the original HRU was 3.27, which decreased to 2.82 (-13.66%) after morphological denoising and 2.86 (-12.45%) after wavelet denoising, indicating that the shape of the HRU became more regular after denoising, reducing abnormal distortion caused by noise. The original mean was 99.86, which increased to 135.23 (+35.42%) after morphological denoising and to 143.74 (+43.95%) after wavelet denoising, indicating that the HRU shape is fuller after denoising and that elongated or irregular fragments are reduced. The original variance of compactness was 0.21, which slightly increased to 0.21 (+0.08%) after morphological denoising and decreased to 0.20 (-4.23%) after wavelet denoising, indicating that the HRU shape is more uniform after wavelet denoising and that abnormally compact / loose units caused by noise are reduced. The original variance of fractal dimension was 0.06, which decreased to 0.05 (-14.22% / -14.17%) after morphological / wavelet denoising, indicating that the smoothness of the HRU boundary is improved and abnormally tortuous boundaries are reduced after denoising.

[0128] Table 1 Comparison of mean HRU index values ​​before and after noise reduction.

[0129] index raw data mean Morphological noise reduction mean Wavelet noise reduction mean Morphological noise reduction mean change rate (%) Wavelet noise reduction mean change rate (%) Number of polygons n 7.32 4.37 4.37 -40.26 -40.28 <![CDATA[Shape Index I S > 3.27 2.82 2.86 -13.66 -12.45 <![CDATA[Compactness I C > 0.20 0.24 0.23 19.88 15.71 <![CDATA[Fractal Dimension I D > 1.35 1.32 1.32 -2.39 -2.46 <![CDATA[Area - perimeter ratio I R > 99.86 135.23 143.74 35.42 43.95

[0130] Table 2 Comparison of standard deviations of HRU indicators before and after noise reduction.

[0131] index Standard deviation of raw data Morphological noise reduction standard deviation Standard deviation of wavelet noise reduction Morphological noise reduction standard deviation change rate (%) Wavelet noise reduction standard deviation change rate (%) Number of polygons n 9.16 4.94 4.92 -46.09 -46.29 <![CDATA[Shape Index I S > 1.73 1.49 1.51 -13.83 -12.8 <![CDATA[Compactness I C > 0.21 0.21 0.2 0.08 -4.23 <![CDATA[Fractal Dimension I D > 0.06 0.05 0.05 -14.22 -14.17 <![CDATA[Area perimeter ratio I R > 95.09 122.25 129.23 28.56 35.91

[0132] According to one aspect of this application, step S2 includes the following sub-steps:

[0133] Step S2.1 Extraction of the water system for each DEM source: For each DEM data source to be evaluated, extract the water system sequentially using the TauDEM tool: eliminate the interference of depressions on the direction of water flow, calculate the flow direction grid, cumulative flow grid, and catchment area grid, and extract the river network according to the set threshold.

[0134] Step S2.2 River system matching degree assessment and data source selection: Compare and analyze the river network data extracted from each DEM in S2.1 with the standard river system data, and select a suitable data source.

[0135] According to one aspect of this application, step S3 includes the following sub-steps.

[0136] Step S3.1 Calculate the slope raster based on the original DEM and obtain the terrain slope raster using the D∞ algorithm. This avoids over-smoothing the terrain by filling depressions and provides a realistic slope base for subsequent noise reduction.

[0137] Step S3.2 Morphological denoising of the slope raster: Morphological filtering is performed on the terrain slope raster. First, specify the filter kernel B. k1 Using a rectangular kernel, k1 is the kernel size. When k1=3, it is a 3×3 structuring element.

[0138] B k1 ={(p,q)∣p,q ∈{- k1,- k1+1,…,0,…,k1-1,k1} (1);

[0139] In the formula, p is the offset in the row direction (vertical relative to the center reference point); q is the offset in the column direction (horizontal relative to the center reference point).

[0140] Using a filter kernel, we perform an erosion-dilation operation on the original two-dimensional matrix S, i.e., an opening operation, denoted as ⊕, to remove small bump noise. Similarly, using a filter kernel, we perform a dilation-erosion operation on the original two-dimensional matrix S, i.e., a closing operation, denoted as ⊕. Fill in small concave noise. Based on the normalized distance weight map generated by the river network buffer, the final output is the linear fusion result S of the two according to geographical weights. OA :

[0141] (2);

[0142] In the formula, i and j are the row and column numbers in the elevation matrix; Z l Z represents a weak denoising path, suitable for use near river networks; s B represents a strong denoising (smoothing) path, suitable for buffer edges; k ω is the filter kernel; ω(i,j) ∈[0,1] is the normalized distance weight map generated based on the river network buffer.

[0143] During the process, the invalid value region remains unchanged, and only the effective topographic region is corrected, so as to eliminate high-frequency small-scale noise while preserving macroscopic topographic features.

[0144] Step S3.3 Wavelet transform denoising of the slope raster: For the slope raster output by D∞, wavelet denoising is performed: a wavelet basis (e.g., db4) and a decomposition level (e.g., 3 levels) are selected to separate the slope data into approximation coefficients (low frequency, representing the overall terrain trend) and detail coefficients (high frequency, containing noise and details); soft thresholding is applied to the detail coefficients to suppress noise, and then the denoised slope raster is reconstructed. Key terrain details such as valley edges and terraces are preserved, and river network spatial information is used to control the regionalized denoising intensity, resulting in the final D∞. WA .

[0145] (3);

[0146] In the formula The slope raster is used to remove invalid values; ψ is the selected wavelet basis; L is the number of decomposition levels; W is the two-dimensional inverse wavelet transform; W -1 For wavelet decomposition operation; A L For W -1 The approximation coefficients of the Lth level of decomposition, Dm k, are obtained after W -1 The k-th layer after decomposition; τ(·;T) represents the soft thresholding function, Tstdk is the original soft threshold, and Tdaptdk is the adjusted soft threshold; m-direction detail coefficients; H, V, and D are the horizontal, vertical, and diagonal directions, respectively; ω(i,j) ∈[0,1] is the normalized distance weight map generated based on the river network buffer; α c ∈[0,1] represents the denoising intensity factor of the central region; α e ∈[0,1] represents the noise reduction intensity factor for the edge region.

[0147] Step S3.4 HRU partitioning: Using the original slope raster and the raster after noise reduction by the two methods, land use and soil attribute data are superimposed according to the slope classification rules to divide the continuous area with consistent land use and soil characteristics within the same slope range into HRUs.

[0148] According to one aspect of this application, step S4 specifically includes:

[0149] Step S4.1 Topographic Index Calculation: Based on the DEM selected in S2, after filling depressions, and according to the river network generated in S2, the river network is forced to the lowest elevation. The slope gradient β of the cell is calculated by the ratio of the elevation difference to the distance between adjacent cells. According to the gradient direction, the catchment area of ​​the current cell is distributed to the downstream cells according to the flow rate ratio, and the total catchment area α of each cell is accumulated. The topographic index raster is generated by the formula ATB=ln(α / tanβ), which reflects the comprehensive characteristics of unit catchment area and slope (the larger the value, the stronger the runoff generation capacity).

[0150] Step S4.2 For the ATB grid, an adaptive region-based, low-interference noise reduction strategy is adopted:

[0151] Adaptive morphological filtering: First, calculate the mean G(i,j) of the terrain gradient in the x and y directions for each grid point. mean ,

[0152] (4);

[0153] Where i and j are the grid point numbers in the x and y directions, respectively, and the gradients in the x and y directions are calculated using the central difference method:

[0154] (5);

[0155] In the formula, d is the grid spacing.

[0156] With the set gradient threshold T split Distinguish between steep regions (gradient greater than a threshold) and gentle regions. Steep regions:

[0157] (6);

[0158] Where M mark This is the classification matrix after thresholding.

[0159] For steep areas, use small kernel B. S (e.g., 3×3) Perform opening and closing operations, removing only small bump noise to avoid destroying details of steep terrain; for flat areas: use a large kernel B. L (e.g., 5×5) Opening and closing operations are performed to smooth small-scale concave noise and preserve macroscopic water flow trends. At the same time, the normalized distance weight map generated based on the river network buffer protects the details around the river network, resulting in the morphologically denoised topographic index matrix ATB. (1) .

[0160] (7);

[0161] In the formula, i and j are the row and column numbers in the elevation matrix, Z l Z represents a weak denoising path, suitable for use near river networks; s This indicates a strong denoising (smoothing) path, suitable for buffer edges. B s B l ω represents the filter kernel corresponding to steep and gentle slope regions. ω(i,j) ∈[0,1] is the normalized distance weight map generated based on the river network buffer.

[0162] Conservative wavelet denoising: The same operation as in equation (3) is used, but a low-level decomposition and a low threshold coefficient are used to suppress only high-frequency random noise and retain the approximate coefficients that characterize the overall trend of the terrain. A soft threshold is applied to the detail coefficients, and a lower threshold is obtained by using the threshold coefficient to avoid excessive smoothing that leads to blurred valley edges. After reconstruction, the grid size is adjusted to ensure consistency with the original terrain index grid. The result after conservative wavelet denoising is denoted as ATB. (2) .

[0163] (8);

[0164] Where T is the adjusted soft threshold; σ d f is the standard deviation of the detail coefficients obtained from wavelet decomposition; s is the threshold coefficient; N is the total number of elements in the terrain index data (two-dimensional matrix).

[0165] Finally, a hybrid weighted denoising method is used: the wavelet and morphological denoising results are fused, and the original topographic index data is used as the basis (weight W0). The results of morphological filtering (weight W1) and conservative wavelet denoising (weight W2) are superimposed, and the noise suppression and detail preservation are balanced by weighted fusion. The final output is a topographic index raster with reduced noise and complete key topographic features.

[0166] (9);

[0167] Where ATB is the result of the original terrain index calculation; ATB (1) The result after morphological denoising; ATB (2) The result is the result after conservative wavelet denoising. W0, W1, and W2 are the weights of the three, and their sum is 1.

[0168] According to one aspect of this application, step S5 specifically includes:

[0169] Step S5.1 Visual comparison and difference analysis of terrain parameters before and after noise reduction: To intuitively evaluate the noise reduction effect, the results of slope and terrain index (ATB) before and after noise reduction can be visualized and compared visually; the difference between the raster before and after noise reduction is calculated and the difference distribution is displayed through raster rendering.

[0170] Step S5.2 Quantitative assessment of HRU fragmentation and verification of noise reduction effect: Based on the HRU vector data (Multipolygon) generated before and after noise reduction, where S is the area of ​​each element in two-dimensional space and P is the perimeter of each element in two-dimensional space, relevant indicators are calculated to verify the effect of noise reduction on HRU quality improvement. Specifically, this includes: ① Number of polygons n: Extracting the number of sub-polygons contained in each HRU (Multipolygon), reflecting the degree of fragmentation within the HRU; ② Shape index I SA larger value indicates a more complex shape; ③ Compactness I C The closer the value is to 1, the more compact the shape; ④ Fractal dimension I D The closer the value is to 1, the smoother the boundary; ⑤ Area to perimeter ratio I R A larger value indicates a fuller HRU.

[0171] Shape Index I S The calculation formula is: (10);

[0172] The formula for calculating compactness IC is: (11);

[0173] The formula for calculating the fractal dimension ID is: (12);

[0174] The formula for calculating the area-to-perimeter ratio (IR) is: (13).

[0175] The average and variance of the above indicators were calculated for the HRU sets before and after noise reduction: a decrease in the average (e.g., a decrease in the average number of polygons and a decrease in the average shape index) indicates that the overall complexity of the HRU is reduced and the fragmentation is alleviated; a decrease in the variance (e.g., a reduction in the compactness variance) indicates that the shape consistency of the HRU is improved and the abnormal units caused by noise are reduced.

[0176] In this embodiment, considering the above indicators, both denoising methods significantly reduce the fragmentation of the HRU: morphological denoising reduces the number of polygons and shape index significantly through global noise suppression, but large kernel filtering in flat areas may cause a few shape anomalies; wavelet denoising, while preserving local details, achieves better area-to-perimeter ratio and compactness through low-level decomposition. Ultimately, wavelet-denoised HRUs were selected as the subsequent hydrological simulation unit, as they better balance shape consistency and topographic detail preservation while reducing fragmentation. A visual comparison of HRUs is shown in Figure 6.

[0177] In summary, the quantization data in this embodiment show that the wavelet denoising path (Slope) wave The method of this invention demonstrates excellent or balanced effects in reducing HRU fragmentation and improving morphological regularity and consistency, confirming that the method of this invention effectively improves the quality of input data for high-precision hydrological models.

[0178] To address the shortcomings of traditional filtering methods, such as the lack of geographic semantic guidance and the tendency to obscure key geomorphic structures during denoising, this invention introduces a spatial adaptive modulation mechanism based on river network line elements. Specifically, a buffer mask is constructed based on the extracted river network vectors, and the distance from each raster cell to the nearest river network line is calculated to generate a normalized distance weight map. On this basis, the denoising intensity is dynamically adjusted, with lower filtering intensity applied to the central region of the line elements to preserve original terrain details, while the denoising effect is gradually enhanced in the edge regions of the buffer. This mechanism achieves a fine balance between noise suppression and terrain fidelity, effectively avoiding excessive smoothing of valley morphology or artifact residue. Furthermore, the design of this strategy fully considers the actual needs of hydrological models: since river networks are usually handled by independent channel modules, and slope runoff calculations rely more on the terrain accuracy of non-channel areas, and given the relatively small area occupied by the river network, preserving its original morphology can minimize the impact on the overall runoff simulation.

[0179] To address the issues of existing terrain processing workflows relying on human experience for hyperparameters and lacking integration with downstream hydrological applications, a closed-loop adaptive optimization framework was constructed. This framework uses the outer objective function as its core, quantifying and fusing key application indicators such as the fragmentation of hydrological response units (HRUs). It iteratively adjusts the hyperparameter combinations of upstream processing modules (such as morphological kernel size, wavelet threshold factor, and multi-candidate fusion weights) through an optimizer, and continuously improves based on feedback from the objective function until the optimal terrain output is obtained.

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

Claims

1. A DEM adaptive multi-algorithm fusion and enhancement method for high-precision distributed hydrological models, characterized in that, Includes the following steps: Step S1: Collect multi-source DEM data of the target watershed, as well as land use and soil attribute data; Step S2: Conduct water system matching assessment and data source screening for multi-source DEM data; Step S3: Fill the depression DEM, calculate the slope grid, and use morphological denoising and wavelet denoising; combine the slope zoning data before and after denoising, as well as land use and soil attribute data, to divide the hydrological response units. Step S4: Calculate the upslope contribution area based on the river network data to obtain the topographic index grid. Combine the slope classification results to suppress noise by adjusting the partition weights of the topographic index grid according to the differences in hydrological response between gentle and steep areas. Step S3 includes: Step S3.1: Based on the original DEM, use the D∞ algorithm to obtain the terrain slope raster; Step S3.2: Specify filter kernel B k1 Using a rectangular kernel, where k1 is the kernel size. B k1 ={(p,q)∣p,q ∈{- k1,- k1+1,…,0,…,k1-1,k1} (1); In the formula, p is the offset in the row direction of the subsequent action; q is the offset in the column direction of the subsequent action. Using a filter kernel, an erosion-dilation operation (also known as an opening operation, denoted as ⊕) is performed on the two-dimensional elevation matrix Z to remove small bump noise. (2)； Using a filter kernel, a dilation-erosion operation is performed on the two-dimensional elevation matrix Z, i.e., a closing operation, denoted as . Fill in small depressions and noise: (3)； In the formula, i and j are the row and column numbers in the elevation matrix, and p, q, u, and v are the offsets in the row and column; Step S3.3: Select wavelet basis and decomposition level to separate slope data into approximation coefficients and detail coefficients; apply soft thresholding to detail coefficients to suppress noise, and then reconstruct to obtain a denoised slope raster, preserving key terrain details including valley edges and terraces; Detail coefficient after noise reduction (4); In the formula, The slope raster is used to remove invalid values; ψ is the selected wavelet basis; L is the number of decomposition levels; W is the two-dimensional inverse wavelet transform; W -1 For wavelet decomposition operation; A L For W -1 The approximation coefficients of the Lth level of decomposition, Dm k, are obtained after W -1 After decomposition, the k-th layer and τ represent the soft thresholding operation, where T is the set soft threshold; m represents the detail coefficients; H, V, and D represent the horizontal, vertical, and diagonal directions, respectively. Step S3.4: Using the original slope raster and the raster after noise reduction by the two methods respectively, land use and soil attribute data are superimposed according to the slope classification rules to divide the continuous area with consistent land use and soil characteristics within the same slope range into HRUs. Step S4 includes: Step S4.1: Based on the filtered DEM, after performing the depression filling operation, according to the generated river network, force the river network to the lowest elevation point, and calculate the cell slope gradient β by the ratio of the elevation difference to the distance between adjacent cells. Based on the gradient direction, the catchment area of ​​the current cell is allocated to the downstream cells according to the flow rate ratio, and the total catchment area α of all cells is obtained by summing them up; the terrain index raster is generated by the formula ATB=ln(α / β); Step S4.2: For the ATB grid, an adaptive region-based, low-interference noise reduction strategy is adopted: First, calculate the mean value G(i,j) of the terrain gradient in the x and y directions for each grid point. mean , (5)； Where i and j are the grid point numbers in the x and y directions, respectively, and the gradients in the x and y directions are calculated using the central difference method: (6)； In the formula, d is the grid spacing; The steep and gentle regions are distinguished by a set gradient threshold Tsplit. Steep region: (7)； Where M mark This is the classification matrix after thresholding. For steep regions, use a small-core BS to perform opening and closing operations to remove small bump noise; For flat regions: use large kernel BL for opening and closing operations to smooth small-scale concave noise; (8)； ATB (1) This is the terrain index matrix after morphological denoising; The same operation as in equation (4) is used, but a low-level decomposition and a low threshold coefficient are used to suppress only high-frequency random noise, retain the approximate coefficients that characterize the overall trend of the terrain, apply a soft threshold to the detail coefficients, and use the threshold coefficient to obtain a lower threshold; after reconstruction, the grid size is adjusted to ensure consistency with the original terrain index grid, and the result after conservative wavelet denoising is denoted as ATB. (2) ; (9)； Where T is the adjusted soft threshold; σ d f is the standard deviation of the detail coefficients obtained from wavelet decomposition; s is the threshold coefficient; N is the total number of elements in the terrain index data; By fusing wavelet and morphological denoising results, based on the original topographic index data W0, the results of morphological filtering W1 and conservative wavelet denoising W2 are superimposed, and the output topographic index raster with reduced noise and complete key topographic features is obtained by weighted fusion to balance noise suppression and detail preservation. (10)； Where ATB is the result of the original terrain index calculation; ATB (1) The result after morphological denoising; ATB (2) The result is the result after conservative wavelet denoising. W0, W1, and W2 are the weights of the three, and their sum is 1.

2. The DEM adaptive multi-algorithm fusion improvement method for high-precision distributed hydrological models according to claim 1, characterized in that, It also includes step S5, comparing the statistical characteristics of slope and topographic index before and after noise reduction to verify the effectiveness of the method and its ability to preserve topographic features, and verifying the quality improvement effect in conjunction with HRU fragmentation, specifically: Step S5.1: Visual comparison and difference analysis of terrain parameters before and after noise reduction. To intuitively evaluate the noise reduction effect, the results of slope and terrain index before and after noise reduction are visualized and compared; the difference between the raster before and after noise reduction is calculated, and the difference distribution is displayed through raster rendering. Step S5.2: Based on the HRU vector data generated before and after noise reduction, S is the area of ​​each element in two-dimensional space, and P is the perimeter of each element in two-dimensional space. Calculate relevant metrics to verify the effect of noise reduction on HRU quality improvement, including: Number of polygons n, Shape Index (11); Compactness (12); fractal dimension (13); Area to perimeter ratio (14).

3. A DEM adaptive multi-algorithm fusion and enhancement method for high-precision distributed hydrological models, characterized in that, include: Obtain raw DEM data, land use data, and soil type data; The original DEM data is normalized to obtain a unified grid DEM; Calculate the basic topographic field based on a unified grid DEM; Based on the base terrain field, at least one enhanced terrain parameter, including a slope raster and a terrain index, is generated through a predefined optimized pipeline. Hydrological response units are delineated by combining at least one enhanced topographic parameter, land use data, and soil type data; Among them, the enhanced terrain parameter is the enhanced slope grid; The steps for generating at least one enhanced terrain parameter include: Obtain the slope grid from the basic topographic field; Morphological filtering is applied to the slope raster to generate a morphologically denoised slope raster. Wavelet transform denoising is applied to the morphologically denoised slope raster to obtain an enhanced slope raster. Among them, the enhanced terrain parameter is the enhanced terrain index; The steps for generating at least one enhanced terrain parameter include: Calculate the initial topographic index based on the basic topographic field; Based on the slope information in the basic topographic field, topographic zones are established, which include gentle zones and steep zones. The initial terrain index is subjected to a zonal adaptive morphological filter, in which a first-size filter kernel is used in flat areas and a second-size filter kernel is used in steep areas. The first size is larger than the second size, thereby generating an adaptive morphologically filtered terrain index. Perform conservative wavelet denoising on the initial terrain index to generate a wavelet-denoised terrain index; The initial terrain index, the terrain index after adaptive morphological filtering, and the terrain index after wavelet denoising are weighted and fused to obtain the enhanced terrain index. The method further includes: Based on the divided hydrological response units, quality indicators are calculated, including at least the hydrological response unit fragmentation index. Define the outer objective function based on quality indicators; One or more hyperparameters used in the step of generating at least one enhanced terrain parameter are defined as the set of parameters to be optimized. An optimizer is employed to optimize the outer objective function. It iteratively executes steps to generate at least one enhanced topographic parameter and to divide the hydrological response units, automatically adjusting the set of parameters to be optimized until the outer objective function converges.

4. The method according to claim 3, characterized in that, Generate at least one enhanced terrain parameter, including: A buffer mask is constructed based on the river network vector, and a normalized distance weight map ω(i,j) is generated. The weight map is embedded into the operation model of morphological filtering and wavelet denoising to perform spatial adaptive denoising on the topographic index grid. The central region of the river network is weakly filtered to preserve the true landform, while the edge region of the buffer is enhanced with denoising intensity, thereby suppressing noise while maintaining the integrity of key hydrological structures.

Images