Distributed model-oriented method for sub-basin division of small and medium-sized river basins
The sub-basin delineation method, which utilizes watershed morphology analysis and dynamic threshold optimization, addresses the scientific rigor and adaptability issues in sub-basin delineation under complex terrain, improves the accuracy and stability of hydrological simulation, and is applicable to watershed management under different terrain and hydrological conditions.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- HEFEI UNIV OF TECH
- Filing Date
- 2025-02-13
- Publication Date
- 2026-06-09
Smart Images

Figure CN120087266B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of water resources management technology, and in particular to a method for dividing small and medium-sized watersheds into sub-basins for distributed models. Background Technology
[0002] In the fields of hydrological simulation and watershed management, sub-basin delineation is a key technical step in water resources management, flood forecasting, and ecological environment assessment. The rationality of sub-basin delineation directly affects the accuracy of hydrological model calculations and determines the reliability of watershed hydrodynamic simulations. Currently, traditional sub-basin delineation methods are mainly based on digital elevation model (DEM) data, determining sub-basin boundaries through flow direction calculation, catchment area analysis, and setting hydrological characteristic parameters. However, these methods have certain limitations and are difficult to adapt to complex watershed topographic features and dynamic hydrological processes, resulting in sub-basin delineation results that lack scientific rigor and adaptability.
[0003] Existing methods for sub-basin delineation typically rely on fixed catchment area thresholds. This involves setting an empirical value as the minimum catchment area threshold for a sub-basin; areas exceeding this threshold are identified as sub-basins, while areas below are merged into adjacent sub-basins. While this method provides some reference for delineation in regular terrain, it struggles to accurately depict the morphology of watersheds with complex geomorphological features, such as mountainous, hilly, and unevenly distributed lake areas. Furthermore, the fixed threshold lacks adaptability and cannot be flexibly adjusted for different watershed types, leading to significant deviations between the delineation results and actual hydrological processes.
[0004] In watershed delineation, traditional methods primarily rely on flow direction calculations based on DEM data, employing the D8 (unidirectional flow) algorithm. While these methods can describe flow direction to some extent, they may produce unreasonable flow direction results in flat areas, depressions, or areas with significant elevation errors. This can lead to discontinuous sub-watershed boundaries or even "reverse flow" phenomena, affecting the accuracy of hydrological simulations. Furthermore, existing flow direction calculation methods often neglect the hydrodynamic connectivity between sub-watersheds, considering only topographic features and failing to dynamically optimize and adjust sub-watersheds in conjunction with watershed hydrological processes, further reducing the rationality of the delineation.
[0005] Another key challenge in subbasin delineation lies in boundary optimization. Traditional subbasin delineation methods typically use threshold segmentation based directly on catchment area, without considering the smoothness and coherence of subbasin boundaries. This results in numerous "fragmented" subbasins—regions that are too small but form independent domains. These small subbasins can cause numerical instability in hydrological simulations, affecting the overall stability of hydrodynamic calculations. Some studies have attempted to optimize the delineation results by manually adjusting subbasin boundaries or merging small subbasins, but these methods rely on manual settings and are difficult to adapt to different types of watershed characteristics.
[0006] In recent years, with the development of Geographic Information System (GIS) technology, artificial intelligence, and machine learning, some scholars have attempted to introduce data-driven methods to optimize the sub-basin delineation process. For example, some studies have combined clustering algorithms for sub-basin optimization, but these methods mostly rely on empirical parameter adjustments and still have certain limitations. Furthermore, traditional methods have failed to effectively utilize watershed morphological characteristics, such as tortuosity, fractal dimension, and confluence path complexity, resulting in insufficient consideration of the overall watershed morphology during sub-basin delineation and affecting the scientific validity of the delineation.
[0007] Therefore, how to provide a method for dividing small and medium-sized watersheds into sub-basins for distributed models is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention
[0008] One objective of this invention is to propose a method for subdividing small and medium-sized watersheds for distributed models. This invention fully utilizes technologies such as watershed morphological feature analysis, multi-scale feature fusion, dynamic threshold optimization, hydrological network topology construction, and neural network optimization for subdividing. It describes in detail the calculation method for intelligent subdividing of subdistricts, which has the advantages of strong adaptability, high calculation accuracy, good boundary optimization continuity, and strong hydrological simulation stability.
[0009] The method for subdividing small and medium-sized watersheds according to a distributed model based on an embodiment of the present invention includes the following steps:
[0010] S1. Obtain digital elevation model data and remote sensing image data of the study area, perform depression filling processing on the digital elevation model data, extract preliminary river network information, and identify the main stream area;
[0011] S2. Extract the actual water system based on remote sensing image data, overlay and analyze it with preliminary river network information, identify the deviation area of the main stream grid, and calculate the elevation difference between the main stream grid and the surrounding depression grid.
[0012] S3. Based on the calculated elevation difference, local elevation adjustment is performed on the digital elevation model data, and depression filling, flow direction calculation and cumulative runoff calculation are re-executed to optimize the river network extraction results.
[0013] S4. Based on the optimized river network extraction results, calculate the watershed morphology parameters, and classify the watershed based on the manifold learning method to identify different watershed types;
[0014] S5. Construct a sub-basin optimal threshold calculation model based on watershed morphological parameters, calculate the catchment area threshold, and dynamically adjust the sub-basin division threshold.
[0015] S6. Based on the catchment area threshold, perform flow direction calculation and confluence analysis on the adjusted digital elevation model data, divide the watershed into sub-basins, optimize the boundary by combining watershed morphology parameters, and perform merging processing on small-area sub-basins.
[0016] S7. Import the sub-basin division results into the distributed hydrological model, analyze the impact of sub-basin division on watershed runoff simulation, and adjust the sub-basin division parameters based on the hydrological simulation results.
[0017] Optionally, S2 specifically includes:
[0018] S21. Based on remote sensing image data, extract the actual water system, and use a water body identification model based on graph attention network and watershed hydrodynamic features to detect water bodies. Combine spectral features, topographic features and hydrodynamic simulation information to calculate the water body classification confidence S. w :
[0019]
[0020] Where σ is the normalized activation function, The adjacency weights are calculated based on the attention mechanism. Let ω represent the water body characteristics of node j at the k-th scale, ω represent the hydrodynamic information weights, and F(i,j) represent the hydrodynamic characteristics.
[0021] S22. Spatial interpolation processing is performed on the extracted preliminary river network information, and the elevation value Z of the main stream grid is calculated using a high-order interpolation method with adaptive weights. r (i,j), stored as interpolated river network raster data R(i,j):
[0022]
[0023] Where Z(m,n) is the known elevation value of the surrounding grid, and μ is the adjustment coefficient. For the second-order Laplace smoothing term, w m,n For adaptive weights, M and N are the number of grid cells, and Z(i,j) is the elevation value of the digital elevation model data at grid coordinate (i,j);
[0024] S23. Binarize the water body classification confidence scores to generate water body distribution raster data W(i,j), and perform raster overlay operation with the obtained interpolated river network raster data R(i,j). Calculate the consistency error using cross-entropy error.
[0025]
[0026] Where E is the consistency error value, if E>T e If so, it is determined to be a main stream grid deviation area, T e Let ln be the error threshold, and ln be the logarithmic function.
[0027] S24. Calculate the elevation difference for the identified main stream grid deviation area and solve for the elevation difference ΔZ. max :
[0028] ΔZ max =max(|Z r (i,j)-Z w (i,j)|)+γ·σ Z ;
[0029] Where, σ Z Z represents the standard deviation of elevation in the adjacent area, γ is the adaptive correction coefficient, max is the maximum value calculation, and Z... r (i,j) represents the grid elevation value of the main stream of the river network, Z w (i,j) represents the elevation values of the surrounding depression grid;
[0030] S25. Store the calculated main stream grid deviation area and elevation difference data.
[0031] Optionally, S3 specifically includes:
[0032] S31. Based on the calculated local elevation difference ΔZ max The elevation of the grid cells in the deviation area is adjusted to construct a local elevation adjustment model:
[0033] Z′ r (i,j)=Z r (i,j)-α·ΔZ max ;
[0034] Among them, Z′ r (i,j) represents the adjusted grid elevation values of the main stream of the river network, α is the local elevation adjustment coefficient, and Z r (i,j) represents the grid elevation values of the main stream of the river network before adjustment;
[0035] S32. A Markov random field model is used to smooth the elevation adjustment area, and a local elevation optimization model is constructed:
[0036]
[0037] in, For the optimized raster elevation values of the main stream of the river network, argmin Z′ To find the candidate elevation adjustment value Z′ that minimizes the objective function value, Z′ r (m,n) represents the adjusted elevation value of the neighboring raster (m,n), λ is the regularization parameter, H(Z′) is the regularization term for elevation changes, N(i,j) is the neighborhood set of the raster, Z′ is the candidate elevation adjustment value, m and n are the row index and column index of the neighboring raster, and β is the neighborhood smoothing factor.
[0038] S33, Based on optimized digital elevation model data Perform a depression-filling operation to fill all internal depressions and calculate the adjusted digital elevation model data:
[0039]
[0040] Among them, Z f (i,j) represents the elevation values in the digital elevation model data after filling the depression, Z w (i,j) represents the elevation values of adjacent depression grid cells;
[0041] S34. Based on the generated digital elevation model data Z after filling the depressions. f (i,j), the flow direction matrix D(i,j) is calculated using the direction-weighted flow direction calculation method:
[0042]
[0043] Among them, Z f (m,n) represents the elevation value of the neighboring raster (m,n), d m,n argmax is the Euclidean distance between grid (i,j) and its neighboring grid (m,n), ε is the flow direction adjustment exponent, and argmax is the Euclidean distance between grid (i,j) and its neighboring grid (m,n). (m,n)∈N(i,j) To determine which adjacent grid cell (m,n) the water flows from the current grid cell (i,j);
[0044] The adaptive runoff weighting model is used to calculate the catchment area of the river network grid.
[0045]
[0046] Where A(i,j) is the catchment area of grid (i,j), exp is the exponential function, ∈ is the catchment weight adjustment parameter, A(m,n) is the catchment area of neighboring grid (m,n), and d p,q Let be the Euclidean distance between grid (i,j) and its neighboring grid (p,q);
[0047] S35. Based on the calculated flow direction matrix and cumulative runoff, optimize the river network extraction results:
[0048]
[0049] in, To extract raster data from the optimized river network, T A This is the minimum catchment area threshold extracted from the river network.
[0050] Optionally, S4 specifically includes:
[0051] S41. Based on the optimized river network, extract raster data and calculate watershed morphological parameters, including curvature C, fractal dimension D, confluence path complexity P, and topological centrality T. c Construct a watershed morphological feature set:
[0052]
[0053] Among them, L r L represents the actual length of the main channel of the basin. s The shortest straight-line distance of the main stream of the basin. In scale Minimum number of grids required to cover the watershed boundary. Where Q is the total number of confluence paths and L is the grid scale. i Let ζ be the length of the i-th confluence path within the sub-basin, and Q be the hydrodynamic adjustment factor. max and Q min Let N(i) be the maximum and minimum flow values within sub-basin i, and let W be the set of sub-basins connected to sub-basin i. ij W represents the convergence weight between sub-basin i and sub-basin j. ik The convergence weight between sub-watershed i and sub-watershed K;
[0054] S42. Based on the calculated watershed morphological parameters, construct the watershed morphological feature matrix M. f :
[0055]
[0056] Among them, C n D n ,P n ,T cn The tortuosity, fractal dimension, confluence path complexity, and topological centrality of the nth sub-basin, where N is the total number of sub-basins;
[0057] S43. Based on the watershed morphological feature matrix, a graph neural network based on an adaptive attention mechanism is used for watershed classification to construct a morphological classification model:
[0058]
[0059] in, Let i be the feature vector of watershed i in the (l+1)th layer of the network. Let σ be the feature vector of watershed j in the l-th layer network, and σ be the nonlinear activation function. Here, W represents the adjacency weights calculated based on the attention mechanism, where exp is the exponential function, τ is the weight vector in the attention mechanism, and W is the adjacency weights. (l( Let be the transformation matrix of the l-th layer network, and LeakyReLU be the leakage-corrected linear unit activation function. Let i be the feature vector of watershed i in the l-th layer network. Let be the feature vector of watershed k in the l-th layer network;
[0060] S44. Based on the calculated watershed morphology classification results, construct a watershed classification label matrix based on soft clustering and neighborhood optimization, and use a joint optimization method to calculate the final classification labels.
[0061]
[0062] Wherein d(T) n M f ) represents the morphological characteristics and category T of sub-basin n. n The Euclidean distance of the centroid, where K is the total number of categories. For temperature regulation parameters, d(T) j M f ) represents the sub-basin j and the morphological feature matrix M f Category T j The Euclidean distance of the centroid, d(T) k M f ) represents the sub-basin k and the morphological feature matrix M f Category T k The Euclidean distance of the centroid, ω n,j For hydrological connectivity weights, T n Here, μ represents the classification label, and μ represents the neighborhood influence factor. To extract from category label T n Choose the category that gives the highest total score;
[0063] S45. Based on the obtained optimized watershed classification labels, generate a visual watershed classification map and store the classification results.
[0064] Optionally, S5 specifically includes:
[0065] S51, Classification results based on storage and watershed morphological feature matrix M fA multi-scale feature fusion model was constructed, incorporating watershed morphology, topographic gradient, and hydrodynamic features, to calculate the catchment area threshold of sub-watersheds.
[0066]
[0067] Among them, F k For the comprehensive morphological characteristics of sub-basin k, β k To optimize the weighting coefficients, N(n) is the set of neighboring sub-basins of sub-basin n, ω n,j For the neighborhood influence weight, F j For the comprehensive morphological characteristics of sub-basin j, α o The neighborhood-weighted influence factor;
[0068] S52. Based on the calculated sub-basin catchment area threshold, and combined with the sub-basin information entropy, calculate the final dynamic threshold used for sub-basin delineation.
[0069]
[0070] Where η and α e E is the dynamic adjustment coefficient. A For the sub-basin information entropy, The average catchment area threshold. The standard deviation of the catchment area threshold;
[0071] S53. Optimize sub-basin partitioning based on computational dynamic thresholds:
[0072]
[0073] Where, ω n,j N represents the neighborhood clustering weights. s T represents the number of sub-basin divisions. A,min T is the minimum area threshold for the sub-basin. A,max This is the threshold for the maximum area of the sub-basin. is the catchment area threshold for the neighboring sub-watershed j;
[0074] S54. Optimize sub-basin data based on computation, store the final sub-basin delineation results, and output the sub-basin delineation matrix M. sub :
[0075]
[0076] Among them, X N and Y N Let n be the geographic coordinates of the sub-basin n. The final catchment area threshold for sub-basin n;
[0077] S55. Store the calculated sub-basin partition data to support the sub-basin partition calculation in step S6.
[0078] Optionally, S6 specifically includes:
[0079] S61. Based on the catchment area threshold, calculate the multi-scale hydrological topological flow direction matrix D(i,j) and construct the watershed topological network:
[0080] D(i,j)=argmax (m,n)∈N(i,j) (Z f (i,j)-Z f (m,n))+γ1·Φ(i,j);
[0081] Among them, Z f (i,j) represents the elevation value of the raster (i,j), Z f (m,n) represents the elevation value of the neighboring grid (m,n) after filling depressions, γ1 is the hydrological topology optimization factor, Φ(i,j) is the flow clustering topology optimization function, and argmax is the topology optimization function. (m,n)∈N(i,j) To select the most likely direction of water flow from all candidate flow directions in the neighborhood, N(i,j) is the neighborhood of grid (i,j);
[0082] S62. Based on the calculated flow direction matrix, a cumulative catchment area calculation method based on information entropy optimization is adopted:
[0083] A(i,j)=∑ (m,n)∈N(i,j) A(m,n)+W(i,j)·(1+γ2·E A );
[0084] Where A(i,j) is the catchment area of grid (i,j), W(i,j) is the watershed contribution rate weight, γ2 is the information entropy adjustment factor, and E A Let A(m,n) be the information entropy of the sub-basin, and let A(m,n) be the cumulative catchment area of the upstream neighboring grid (m,n).
[0085] S63. Based on the calculated catchment area, combined with the calculated dynamic catchment area threshold. Optimizing sub-basin partitioning using dynamic neural networks:
[0086]
[0087] Among them, M sub (i,j) represents the sub-basin partitioning result matrix, σ is the adaptive classification activation function, and ω n,j γ is the neighborhood influence weight, and γ3 is the neighborhood smoothing adjustment parameter;
[0088] S64. Based on the calculated sub-basin division results, optimize the sub-basin boundaries and adjust the merging relationships of small-area sub-basins:
[0089]
[0090] Where B(i,j) is the optimized sub-basin boundary matrix, δ(M sub (m,n),M sub (i,j)) is the sub-basin merging function, T A,min This represents the minimum area threshold for the sub-basin.
[0091] S65. Based on the calculated sub-basin division results, store the final optimized sub-basin data.
[0092] The beneficial effects of this invention are:
[0093] This invention proposes an adaptive sub-basin delineation method by introducing watershed morphology constraints. Compared to traditional fixed-threshold delineation methods, this invention can adaptively adjust the catchment area threshold of sub-basins, making the delineation results more consistent with the actual hydrological characteristics of different types of watersheds. This invention employs a dynamic threshold optimization method based on multi-scale feature fusion, comprehensively considering watershed curvature, fractal dimension, confluence path complexity, topological centrality, slope, and hydrodynamic characteristics. This ensures that sub-basin delineation does not rely solely on a single hydrological parameter but incorporates the overall morphological characteristics of the watershed, thereby improving the scientific rigor and adaptability of the delineation.
[0094] This invention optimizes the flow direction calculation method by employing a hydrological network topology optimization technique based on flow clustering. This enhances the hydrodynamic connectivity between sub-basins and reduces calculation errors caused by abrupt or unreasonable flow direction changes. Traditional flow direction calculation methods are prone to flow direction breaks or incorrect delineation when dealing with depressions, flat areas, or complex terrain. This invention introduces a flow direction optimization matrix and combines it with information entropy optimization to calculate the cumulative catchment area. This allows sub-basin boundaries to more reasonably adapt to terrain changes, improving the stability and accuracy of hydrological simulation calculations.
[0095] This invention employs an optimized sub-basin partitioning method based on dynamic neural networks. It introduces Wasserstein distance to calculate neighborhood influence weights, ensuring the smoothness of sub-basin partitioning boundaries and avoiding the fragmented sub-basin problem common in traditional methods. By introducing a dynamic activation function during sub-basin partitioning, the continuity and adaptability of sub-basin partitioning are guaranteed, reducing numerical instability caused by boundary jumps during computation and improving the interpretability of the calculation results.
[0096] Furthermore, this invention incorporates the Wasserstein distance merging method during the sub-basin optimization stage, enabling small-area sub-basins to be intelligently merged into their optimal neighboring sub-basins, thus preventing the instability of hydrological calculations from being affected by excessively small factor basins. Simultaneously, for sub-basins with excessively large areas, this invention employs the minimum cut maximum flow method for subdividing the sub-basins, ensuring that the sub-basin division reflects the true distribution of hydrological processes within the basin and improving the computational accuracy of hydrological simulations.
[0097] In summary, this invention provides a more intelligent and adaptive method for sub-basin delineation, overcoming the limitations of fixed threshold delineation methods in existing technologies and improving the rationality of sub-basin boundaries and computational stability. This invention is applicable to watershed delineation under different topographical and hydrological conditions, providing more accurate delineation schemes for water resource management, flood prediction, and hydrological model optimization, thereby enhancing the reliability and practicality of watershed hydrological calculations. Attached Figure Description
[0098] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings:
[0099] Figure 1 This is a flowchart of the method for dividing small and medium-sized watersheds into sub-basins for a distributed model proposed in this invention.
[0100] Figure 2 This is a schematic diagram illustrating the watershed morphological feature extraction process of the sub-watershed division method for small and medium-sized watersheds proposed in this invention, which is oriented towards a distributed model. Detailed Implementation
[0101] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.
[0102] refer to Figure 1-2 The method for subdividing small and medium-sized watersheds in a distributed model includes the following steps:
[0103] S1. Obtain digital elevation model data and remote sensing image data of the study area, perform depression filling processing on the digital elevation model data, extract preliminary river network information, and identify the main stream area;
[0104] S2. Extract the actual water system based on remote sensing image data, overlay and analyze it with preliminary river network information, identify the deviation area of the main stream grid, and calculate the elevation difference between the main stream grid and the surrounding depression grid.
[0105] S3. Based on the calculated elevation difference, local elevation adjustment is performed on the digital elevation model data, and depression filling, flow direction calculation and cumulative runoff calculation are re-executed to optimize the river network extraction results.
[0106] S4. Based on the optimized river network extraction results, calculate the watershed morphology parameters, and classify the watershed based on the manifold learning method to identify different watershed types;
[0107] S5. Construct a sub-basin optimal threshold calculation model based on watershed morphological parameters, calculate the catchment area threshold, and dynamically adjust the sub-basin division threshold.
[0108] S6. Based on the catchment area threshold, perform flow direction calculation and confluence analysis on the adjusted digital elevation model data, divide the watershed into sub-basins, optimize the boundary by combining watershed morphology parameters, and perform merging processing on small-area sub-basins.
[0109] S7. Import the sub-basin division results into the distributed hydrological model, analyze the impact of sub-basin division on watershed runoff simulation, and adjust the sub-basin division parameters based on the hydrological simulation results.
[0110] In this embodiment, S2 specifically includes:
[0111] S21. Based on remote sensing image data, extract the actual water system, and use a water body identification model based on graph attention network and watershed hydrodynamic features to detect water bodies. Combine spectral features, topographic features and hydrodynamic simulation information to calculate the water body classification confidence S. w :
[0112]
[0113] Where σ is the normalized activation function, The adjacency weights are calculated based on the attention mechanism. Let ω represent the water body characteristics of node j at the k-th scale, ω represent the hydrodynamic information weights, and F(i,j) represent the hydrodynamic characteristics.
[0114] S22. Spatial interpolation processing is performed on the extracted preliminary river network information, and the elevation value Z of the main stream grid is calculated using a high-order interpolation method with adaptive weights. r (i,j), stored as interpolated river network raster data R(i,j):
[0115]
[0116] Where Z(m,n) is the known elevation value of the surrounding grid, and μ is the adjustment coefficient. For the second-order Laplace smoothing term, w m,nFor adaptive weights, M and N are the number of grid cells, and Z(i,j) is the elevation value of the digital elevation model data at grid coordinate (i,j);
[0117] S23. Binarize the water body classification confidence scores to generate water body distribution raster data W(i,j), and perform raster overlay operation with the obtained interpolated river network raster data R(i,j). Calculate the consistency error using cross-entropy error.
[0118]
[0119] Where E is the consistency error value, if E>T e If so, it is determined to be a main stream grid deviation area, T e Let ln be the error threshold, and ln be the logarithmic function.
[0120] S24. Calculate the elevation difference for the identified main stream grid deviation area and solve for the elevation difference ΔZ. max :
[0121] ΔZ max =max(|Z r (i,j)-Z w (i,j)|)+γ·σ Z ;
[0122] Where, σ Z Z represents the standard deviation of elevation in the adjacent area, γ is the adaptive correction coefficient, max is the maximum value calculation, and Z... r (i,j) represents the grid elevation value of the main stream of the river network, Z w (i,j) represents the elevation values of the surrounding depression grid;
[0123] S25. Store the calculated main stream grid deviation area and elevation difference data.
[0124] In this embodiment, S3 specifically includes:
[0125] S31. Based on the calculated local elevation difference ΔZ max The elevation of the grid cells in the deviation area is adjusted to construct a local elevation adjustment model:
[0126] Z′ r (i,j)=Z r (i,j)-α·ΔZ max ;
[0127] Among them, Z′ r (i,j) represents the adjusted grid elevation values of the main stream of the river network, α is the local elevation adjustment coefficient, and Z r (i,j) represents the grid elevation values of the main stream of the river network before adjustment;
[0128] S32. A Markov random field model is used to smooth the elevation adjustment area, and a local elevation optimization model is constructed:
[0129]
[0130] in, For the optimized raster elevation values of the main stream of the river network, argmin Z′ To find the candidate elevation adjustment value Z′ that minimizes the objective function value, Z′ r (m,n) represents the adjusted elevation value of the neighboring raster (m,n), λ is the regularization parameter, H(Z′) is the regularization term for elevation changes, N(i,j) is the neighborhood set of the raster, Z′ is the candidate elevation adjustment value, m and n are the row index and column index of the neighboring raster, and β is the neighborhood smoothing factor.
[0131] S33, Based on optimized digital elevation model data Perform a depression-filling operation to fill all internal depressions and calculate the adjusted digital elevation model data:
[0132]
[0133] Among them, Z f (i,j) represents the elevation values in the digital elevation model data after filling the depression, Z w (i,j) represents the elevation values of adjacent depression grid cells;
[0134] S34. Based on the generated digital elevation model data Z after filling the depressions. f (i,j), the flow direction matrix D(i,j) is calculated using the direction-weighted flow direction calculation method:
[0135]
[0136] Among them, Z f (m,n) represents the elevation value of the neighboring raster (m,n), d m,n argmax is the Euclidean distance between grid (i,j) and its neighboring grid (m,n), ε is the flow direction adjustment exponent, and argmax is the Euclidean distance between grid (i,j) and its neighboring grid (m,n). (m,n)∈N(i,j) To determine which adjacent grid cell (m,n) the water flows from the current grid cell (i,j);
[0137] The adaptive runoff weighting model is used to calculate the catchment area of the river network grid.
[0138]
[0139] Where A(i,j) is the catchment area of grid (i,j), exp is the exponential function, ∈ is the catchment weight adjustment parameter, A(m,n) is the catchment area of neighboring grid (m,n), and dp,q Let be the Euclidean distance between grid (i,j) and its neighboring grid (p,q);
[0140] S35. Based on the calculated flow direction matrix and cumulative runoff, optimize the river network extraction results:
[0141]
[0142] in, To extract raster data from the optimized river network, T A This is the minimum catchment area threshold extracted from the river network.
[0143] In this embodiment, S4 specifically includes:
[0144] S41. Based on the optimized river network, extract raster data and calculate watershed morphological parameters, including curvature C, fractal dimension D, confluence path complexity P, and topological centrality T. c Construct a watershed morphological feature set:
[0145]
[0146] Among them, L r L represents the actual length of the main channel of the basin. s The shortest straight-line distance of the main stream of the basin. In scale Minimum number of grids required to cover the watershed boundary. Where Q is the total number of confluence paths and L is the grid scale. i Let ζ be the length of the i-th confluence path within the sub-basin, and Q be the hydrodynamic adjustment factor. max and Q min Let N(i) be the maximum and minimum flow values within sub-basin i, and let W be the set of sub-basins connected to sub-basin i. ij W represents the convergence weight between sub-basin i and sub-basin j. ik The convergence weight between sub-watershed i and sub-watershed K;
[0147] S42. Based on the calculated watershed morphological parameters, construct the watershed morphological feature matrix M. f :
[0148]
[0149] Among them, C n D n ,P n ,T cn The tortuosity, fractal dimension, confluence path complexity, and topological centrality of the nth sub-basin, where N is the total number of sub-basins;
[0150] S43. Based on the watershed morphological feature matrix, a graph neural network based on an adaptive attention mechanism is used for watershed classification to construct a morphological classification model:
[0151]
[0152] in, Let i be the feature vector of watershed i in the (l+1)th layer of the network. Let σ be the feature vector of watershed j in the l-th layer network, and σ be the nonlinear activation function. Here, W represents the adjacency weights calculated based on the attention mechanism, where exp is the exponential function, τ is the weight vector in the attention mechanism, and W is the adjacency weights. (l) Let be the transformation matrix of the l-th layer network, and LeakyReLU be the leakage-corrected linear unit activation function. Let i be the feature vector of watershed i in the l-th layer network. Let be the feature vector of watershed k in the l-th layer network;
[0153] S44. Based on the calculated watershed morphology classification results, construct a watershed classification label matrix based on soft clustering and neighborhood optimization, and use a joint optimization method to calculate the final classification labels.
[0154]
[0155] Wherein d(T) n M f ) represents the morphological characteristics and category T of sub-basin n. n The Euclidean distance of the centroid, where K is the total number of categories. For temperature regulation parameters, d(T) j M f ) represents the sub-basin j and the morphological feature matrix M f Category T j The Euclidean distance of the centroid, d(T) k M f ) represents the sub-basin k and the morphological feature matrix M f Category T k The Euclidean distance of the centroid, ω n,j For hydrological connectivity weights, T n Here, μ represents the classification label, and μ represents the neighborhood influence factor. To extract from category label T n Choose the category that gives the highest total score;
[0156] S45. Based on the obtained optimized watershed classification labels, generate a visual watershed classification map and store the classification results.
[0157] In this embodiment, S5 specifically includes:
[0158] S51, Classification results based on storage and watershed morphological feature matrix M f A multi-scale feature fusion model was constructed, incorporating watershed morphology, topographic gradient, and hydrodynamic features, to calculate the catchment area threshold of sub-watersheds.
[0159]
[0160] Among them, F k For the comprehensive morphological characteristics of sub-basin k, β k To optimize the weighting coefficients, N(n) is the set of neighboring sub-basins of sub-basin n, ω n,j For the neighborhood influence weight, F j For the comprehensive morphological characteristics of sub-basin j, α o The neighborhood-weighted influence factor;
[0161] S52. Based on the calculated sub-basin catchment area threshold, and combined with the sub-basin information entropy, calculate the final dynamic threshold used for sub-basin delineation.
[0162]
[0163] Where η and α e E is the dynamic adjustment coefficient. A For the sub-basin information entropy, The average catchment area threshold. The standard deviation of the catchment area threshold;
[0164] S53. Optimize sub-basin partitioning based on computational dynamic thresholds:
[0165]
[0166] Where, ω n,j N represents the neighborhood clustering weights. s T represents the number of sub-basin divisions. A,min T is the minimum area threshold for the sub-basin. A,max This is the threshold for the maximum area of the sub-basin. is the catchment area threshold for the neighboring sub-watershed j;
[0167] S54. Optimize sub-basin data based on computation, store the final sub-basin delineation results, and output the sub-basin delineation matrix M. sub :
[0168]
[0169] Among them, X N and Y N Let n be the geographic coordinates of the sub-basin n. The final catchment area threshold for sub-basin n;
[0170] S55. Store the calculated sub-basin partition data to support the sub-basin partition calculation in step S6.
[0171] In this embodiment, S6 specifically includes:
[0172] S61. Based on the catchment area threshold, calculate the multi-scale hydrological topological flow direction matrix D(i,j) and construct the watershed topological network:
[0173] D(i,j)=argmax (m,n)∈N(i,j) (Z f (i,j)-Z f (m,n))+γ1·Φ(i,j);
[0174] Among them, Z f (i,j) represents the elevation value of the raster (i,j), Z f (m,n) represents the elevation value of the neighboring grid (m,n) after filling depressions, γ1 is the hydrological topology optimization factor, Φ(i,j) is the flow clustering topology optimization function, and argmax is the topology optimization function. (m,n)∈N(i,j) To select the most likely direction of water flow from all candidate flow directions in the neighborhood, N(i,j) is the neighborhood of grid (i,j);
[0175] S62. Based on the calculated flow direction matrix, a cumulative catchment area calculation method based on information entropy optimization is adopted:
[0176] A(i,j)=∑ (m,n)∈N(i,j) A(m,n)+W(i,j)·(1+γ2·E A );
[0177] Where A(i,j) is the catchment area of grid (i,j), W(i,j) is the watershed contribution rate weight, γ2 is the information entropy adjustment factor, and E A Let A(m,n) be the information entropy of the sub-basin, and let A(m,n) be the cumulative catchment area of the upstream neighboring grid (m,n).
[0178] S63. Based on the calculated catchment area, combined with the calculated dynamic catchment area threshold. Optimizing sub-basin partitioning using dynamic neural networks:
[0179]
[0180] Among them, M sub (i,j) represents the sub-basin partitioning result matrix, σ is the adaptive classification activation function, and ω n,j γ is the neighborhood influence weight, and γ3 is the neighborhood smoothing adjustment parameter;
[0181] S64. Based on the calculated sub-basin division results, optimize the sub-basin boundaries and adjust the merging relationships of small-area sub-basins:
[0182]
[0183] Where B(i,j) is the optimized sub-basin boundary matrix, δ(M sub (m,n),M sub (i,j)) is the sub-basin merging function, T A,min This represents the minimum area threshold for the sub-basin.
[0184] S65. Based on the calculated sub-basin division results, store the final optimized sub-basin data.
[0185] Example 1:
[0186] To verify the feasibility of this invention in practice, it was applied to a hilly terrain watershed for experimental analysis. Comparative tests were conducted in four aspects: accuracy of sub-watershed division, boundary smoothness, computational efficiency, and hydrological simulation effect.
[0187] The experiment selected a hilly watershed for sub-basin delineation testing. During the experiment, DEM data was first used to obtain the watershed's elevation information, and then the method of this invention was used for sub-basin delineation. Watershed morphology feature analysis was employed to extract key parameters such as watershed tortuosity, fractal dimension, and confluence path complexity. Based on these features, a dynamic catchment area threshold was calculated for adaptive optimization of the sub-basins. Hydrological network topology optimization was used to calculate flow direction, ensuring hydrodynamic connectivity and improving the rationality of the flow direction calculation.
[0188] During the experiment, to compare the advantages and disadvantages of the present invention and traditional methods, sub-basin division was carried out using the fixed threshold method, the GIS-based division method, and the method of the present invention, respectively, and the computational efficiency, sub-basin boundary continuity, and hydrological simulation accuracy were tested. Multiple sub-basins were selected for detailed analysis within each basin, and experimental data were collected for comparison.
[0189] The experimental data comes from recent measured hydrological data, covering the period from 2023 to 2024. The data includes indicators such as rainfall, flow changes, and hydrological simulation errors after sub-basin division at different times.
[0190] Table 1 Comparison of Sub-basin Delineation Effects
[0191]
[0192] Experimental results show that the present invention outperforms traditional methods in terms of sub-basin delineation accuracy, computational efficiency, and hydrological simulation stability. Regarding the optimization of the number of sub-basins, the method of the present invention delineates an average of 187 sub-basins, a reduction of 18.5% compared to traditional methods. This indicates that adaptive dynamic threshold optimization can effectively reduce redundant sub-basins, resulting in a more compact delineation and avoiding the over-delineation problem caused by unreasonable fixed threshold settings in traditional methods. Meanwhile, although the GIS semi-automatic delineation method has been optimized, it still cannot achieve the same effect as the present invention, indicating that traditional methods have limitations in complex terrain.
[0193] Regarding sub-basin boundary optimization, the number of boundary abrupt change points in this invention is reduced to 735, a 41.7% reduction compared to traditional methods and a 32.1% reduction compared to semi-automatic GIS methods. This indicates that this invention, through Wasserstein distance optimization and watershed morphology analysis, improves the continuity of sub-basin boundaries, reduces boundary abrupt changes caused by flow direction calculation errors, and optimizes the coherence between sub-basins. The reduction in boundary abrupt change points means smoother delineation, ensuring the stability of hydrological simulation calculations and avoiding error accumulation caused by unreasonable boundaries in traditional methods.
[0194] In terms of computational efficiency, the method of this invention takes only 46 minutes, which is 51.6% shorter than the traditional method (95 minutes) and 37% faster than the GIS semi-automatic method (73 minutes). This result shows that the present invention significantly reduces computational redundancy and improves the processing capability of large-scale watershed data by optimizing flow direction calculation, adjusting information entropy, and optimizing dynamic neural networks. At the same time, it improves the degree of computational automation and avoids the drawbacks of the GIS semi-automatic method that requires manual intervention.
[0195] Regarding the accuracy of hydrological simulation, the simulation error of this invention is only 8.5%, a reduction of 47.5% compared to traditional methods (16.2%) and a reduction of 34.1% compared to GIS semi-automatic methods (12.9%). This result demonstrates that this invention effectively reduces watershed delineation errors by dynamically optimizing sub-basin thresholds, optimizing flow direction calculations, and using neural networks to optimize boundary delineation, ensuring the accuracy of hydrological simulation calculations and improving the matching degree of sub-basin hydrological characteristics. Lower hydrological simulation errors mean that the delineation results are more consistent with actual hydrological processes, providing reliable support for water resource management, flood prediction, and ecological assessment.
[0196] In summary, this invention achieves significant improvements in sub-basin delineation accuracy, boundary smoothness, computational efficiency, and hydrological simulation stability. Compared to traditional methods, this invention not only adaptively adjusts sub-basin thresholds based on watershed morphological characteristics but also optimizes flow direction calculations, improves computational efficiency, ensures the rationality of sub-basin delineation, and provides more accurate and reliable technical support for watershed hydrological calculations.
[0197] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.
Claims
1. A method for subdividing small and medium-sized watersheds in a distributed model, characterized in that, Includes the following steps: S1. Obtain digital elevation model data and remote sensing image data of the study area, perform depression filling processing on the digital elevation model data, extract preliminary river network information, and identify the main stream area; S2. Extract the actual water system based on remote sensing image data, overlay and analyze it with preliminary river network information, identify the deviation area of the main stream grid, and calculate the elevation difference between the main stream grid and the surrounding depression grid. S3. Based on the calculated elevation difference, local elevation adjustment is performed on the digital elevation model data, and depression filling, flow direction calculation and cumulative runoff calculation are re-executed to optimize the river network extraction results. S4. Based on the optimized river network extraction results, calculate the watershed morphology parameters, and classify the watershed based on the manifold learning method to identify different watershed types; S5. Construct a sub-basin optimal threshold calculation model based on watershed morphological parameters, calculate the catchment area threshold, and dynamically adjust the sub-basin division threshold. S6. Based on the catchment area threshold, perform flow direction calculation and confluence analysis on the adjusted digital elevation model data, divide the watershed into sub-basins, optimize the boundary by combining watershed morphology parameters, and perform merging processing on small-area sub-basins. S7. Import the sub-basin division results into the distributed hydrological model, analyze the impact of sub-basin division on watershed runoff simulation, and adjust the sub-basin division parameters based on the hydrological simulation results.
2. The method for dividing small and medium-sized watersheds into sub-basins for a distributed model according to claim 1, characterized in that, S2 specifically includes: S21. Based on remote sensing image data, extract the actual water system, and use a water body identification model based on graph attention network and watershed hydrodynamic features to detect water bodies. Combine spectral features, topographic features and hydrodynamic simulation information to calculate the water body classification confidence S. w : Where σ is the normalized activation function, For the adjacency weights calculated based on the attention mechanism, Let ω represent the water body characteristics of node j at the k-th scale, ω represent the hydrodynamic information weights, and F(i,j) represent the hydrodynamic characteristics. S22. Spatial interpolation processing is performed on the extracted preliminary river network information, and the elevation value Z of the main stream grid is calculated using a high-order interpolation method with adaptive weights. r (i,j), stored as interpolated river network raster data R(i,j): Where Z(m,n) is the known elevation value of the surrounding grid, and μ is the adjustment coefficient. For the second-order Laplace smoothing term, w m,n For adaptive weights, M and N are the number of grid cells, and Z(i,j) is the elevation value of the digital elevation model data at grid coordinate (i,j); S23. Binarize the water body classification confidence scores to generate water body distribution raster data W(i,j), and perform raster overlay operation with the obtained interpolated river network raster data R(i,j). Calculate the consistency error using cross-entropy error. Where E is the consistency error value, if E>T e If so, it is determined to be a main stream grid deviation area, T e Let ln be the error threshold, and ln be the logarithmic function. S24. Calculate the elevation difference for the identified main stream grid deviation area and solve for the elevation difference ΔZ. max : ΔZ max =max(|Z r (i,j)-Z w (i,j)|)+γ·σ Z ; Where, σ Z Z represents the standard deviation of elevation in the adjacent area, γ is the adaptive correction coefficient, max is the maximum value calculation, and Z... r (i,j) represents the grid elevation value of the main stream of the river network, Z w (i,j) represents the elevation values of the surrounding depression grid; S25. Store the calculated main stream grid deviation area and elevation difference data.
3. The method for dividing small and medium-sized watersheds into sub-basins for a distributed model according to claim 1, characterized in that, S3 specifically includes: S31. Based on the calculated local elevation difference ΔZ max The elevation of the grid cells in the deviation area is adjusted to construct a local elevation adjustment model: Z′ r (i,j)=Z r (i,j)-α·ΔZ max ; Among them, Z′ r (i,j) represents the adjusted grid elevation values of the main stream of the river network, α is the local elevation adjustment coefficient, and Z r (i,j) represents the grid elevation values of the main stream of the river network before adjustment; S32. A Markov random field model is used to smooth the elevation adjustment area, and a local elevation optimization model is constructed: in, For the optimized raster elevation values of the main stream of the river network, argmin Z′ To find the candidate elevation adjustment value Z′ that minimizes the objective function value, Z′ r (m,n) represents the adjusted elevation value of the neighboring raster (m,n), λ is the regularization parameter, H(Z′) is the regularization term for elevation changes, N(i,j) is the neighborhood set of the raster, Z′ is the candidate elevation adjustment value, m and n are the row index and column index of the neighboring raster, and β is the neighborhood smoothing factor. S33, Based on optimized digital elevation model data Perform a depression-filling operation to fill all internal depressions and calculate the adjusted digital elevation model data: Among them, Z f (i,j) represents the elevation values in the digital elevation model data after filling the depression, Z w (i,j) represents the elevation values of adjacent depression grid cells; S34. Based on the generated digital elevation model data Z after filling the depressions. f (i,j), the flow direction matrix D(i,j) is calculated using the direction-weighted flow direction calculation method: Among them, Z f (m,n) represents the elevation value of the neighboring raster (m,n), d m,n argmax is the Euclidean distance between grid (i,j) and its neighboring grid (m,n), ε is the flow direction adjustment exponent, and argmax is the Euclidean distance between grid (i,j) and its neighboring grid (m,n). (m,n)∈N(i,j) To determine which adjacent grid cell (m,n) the water flows from the current grid cell (i,j); The adaptive runoff weighting model is used to calculate the catchment area of the river network grid. Where A(i,j) is the catchment area of grid (i,j), exp is the exponential function, ∈ is the catchment weight adjustment parameter, A(m,n) is the catchment area of neighboring grid (m,n), and d p,q Let be the Euclidean distance between grid (i,j) and its neighboring grid (p,q); S35. Based on the calculated flow direction matrix and cumulative runoff, optimize the river network extraction results: in, To extract raster data from the optimized river network, T A This is the minimum catchment area threshold extracted from the river network.
4. The method for dividing small and medium-sized watersheds into sub-basins for a distributed model according to claim 1, characterized in that, S4 specifically includes: S41. Based on the optimized river network, extract raster data and calculate watershed morphological parameters, including curvature C, fractal dimension D, confluence path complexity P, and topological centrality T. c Construct a watershed morphological feature set: Among them, L r L represents the actual length of the main channel of the basin. s The shortest straight-line distance of the main stream of the basin, N(θ) is the minimum number of grids required to cover the basin boundary at scale θ, where θ is the grid scale, Q is the total number of confluence paths, and L is the total number of confluence paths. i Let ζ be the length of the i-th confluence path within the sub-basin, and Q be the hydrodynamic adjustment factor. max and Q min Let N(i) be the maximum and minimum flow values within sub-basin i, and let W be the set of sub-basins connected to sub-basin i. ij W represents the convergence weight between sub-basin i and sub-basin j. ik The convergence weight between sub-watershed i and sub-watershed K; S42. Based on the calculated watershed morphological parameters, construct the watershed morphological feature matrix M. f : Among them, C n D n ,P n ,T cn The tortuosity, fractal dimension, confluence path complexity, and topological centrality of the nth sub-basin, where N is the total number of sub-basins; S43. Based on the watershed morphological feature matrix, a graph neural network based on an adaptive attention mechanism is used for watershed classification to construct a morphological classification model: in, Let i be the feature vector of watershed i in the (l+1)th layer of the network. Let σ be the feature vector of watershed j in the l-th layer network, and σ be the nonlinear activation function. Here, W represents the adjacency weights calculated based on the attention mechanism, where exp is the exponential function, τ is the weight vector in the attention mechanism, and W is the adjacency weights. (l) Let be the transformation matrix of the l-th layer network, and LeakyReLU be the leakage-corrected linear unit activation function. Let i be the feature vector of watershed i in the l-th layer network. Let be the feature vector of watershed k in the l-th layer network; S44. Based on the calculated watershed morphology classification results, construct a watershed classification label matrix based on soft clustering and neighborhood optimization, and use a joint optimization method to calculate the final classification labels. Wherein d(T) n M f ) represents the morphological characteristics and category T of sub-basin n. n The Euclidean distance of the centroid, where K is the total number of categories. For temperature regulation parameters, d(T) j M f ) represents the sub-basin j and the morphological feature matrix M f Category T j The Euclidean distance of the centroid, d(T) k M f ) represents the sub-basin k and the morphological feature matrix M f Category T k The Euclidean distance of the centroid, ω n,j For hydrological connectivity weights, T n Here, μ represents the classification label, and μ represents the neighborhood influence factor. To extract from category label T n Choose the category that gives the highest total score; S45. Based on the obtained optimized watershed classification labels, generate a visual watershed classification map and store the classification results.
5. The method for dividing small and medium-sized watersheds into sub-basins for a distributed model according to claim 1, characterized in that, S5 specifically includes: S51, Classification results based on storage and watershed morphological feature matrix M f A multi-scale feature fusion model was constructed, incorporating watershed morphology, topographic gradient, and hydrodynamic features, to calculate the catchment area threshold of sub-watersheds. Among them, F k For the comprehensive morphological characteristics of sub-basin k, β k To optimize the weighting coefficients, N(n) is the set of neighboring sub-basins of sub-basin n, ω n,j For the neighborhood influence weight, F j For the comprehensive morphological characteristics of sub-basin j, α o The neighborhood-weighted influence factor; S52. Based on the calculated sub-basin catchment area threshold, and combined with the sub-basin information entropy, calculate the final dynamic threshold used for sub-basin delineation. Where η and α e E is the dynamic adjustment coefficient. A For the sub-basin information entropy, The average catchment area threshold. The standard deviation of the catchment area threshold; S53. Optimize sub-basin partitioning based on computational dynamic thresholds: Where, ω n,j N represents the neighborhood clustering weights. s T represents the number of sub-basin divisions. A,min T is the minimum area threshold for the sub-basin. A,max This is the threshold for the maximum area of the sub-basin. is the catchment area threshold for the neighboring sub-watershed j; S54. Optimize sub-basin data based on computation, store the final sub-basin delineation results, and output the sub-basin delineation matrix M. sub : Among them, X N and Y N Let n be the geographic coordinates of the sub-basin n. The final catchment area threshold for sub-basin n; S55. Store the calculated sub-basin partition data to support the sub-basin partition calculation in step S6.
6. The method for dividing small and medium-sized watersheds into sub-basins for a distributed model according to claim 1, characterized in that, S6 specifically includes: S61. Based on the catchment area threshold, calculate the multi-scale hydrological topological flow direction matrix D(i,j) and construct the watershed topological network: D(i,j)=argmax (m,n)∈N(i,j) (Z f (i,j)-Z f (m,n))+γ1·Φ(i,j); Among them, Z f (i,j) represents the elevation value of the raster (i,j), Z f (m,n) represents the elevation value of the neighboring grid (m,n) after filling depressions, γ1 is the hydrological topology optimization factor, Φ(i,j) is the flow clustering topology optimization function, and argmax is the topology optimization function. (m,n)∈N(i,j) To select the most likely direction of water flow from all candidate flow directions in the neighborhood, N(i,j) is the neighborhood of grid (i,j); S62. Based on the calculated flow direction matrix, a cumulative catchment area calculation method based on information entropy optimization is adopted: A(i,j)=∑ (m,n)∈N(i,j) A(m,n)+W(i,j)·(1+γ2·E A ); Where A(i,j) is the catchment area of grid (i,j), W(i,j) is the watershed contribution rate weight, γ2 is the information entropy adjustment factor, and E A Let A(m,n) be the information entropy of the sub-basin, and let A(m,n) be the cumulative catchment area of the upstream neighboring grid (m,n). S63. Based on the calculated catchment area, combined with the calculated dynamic catchment area threshold. Optimizing sub-basin partitioning using dynamic neural networks: Among them, M sub (i,j) is the sub-basin partitioning result matrix, σ is the adaptive classification activation function, and ω is the sub-basin partitioning result matrix. n,j γ is the neighborhood influence weight, and γ3 is the neighborhood smoothing adjustment parameter; S64. Based on the calculated sub-basin division results, optimize the sub-basin boundaries and adjust the merging relationships of small-area sub-basins: Where B(i,j) is the optimized sub-basin boundary matrix, δ(M sub (m,n),M sub (i,j)) is the sub-basin merging function, T A,min This represents the minimum area threshold for the sub-basin. S65. Based on the calculated sub-basin division results, store the final optimized sub-basin data.