A method for judging similarity of multi-dimensional basins based on geography-meteorology-hydrology

Abstract

This invention discloses a multi-dimensional watershed similarity discrimination method based on geography, meteorology, and hydrology, relating to the fields of hydrological simulation and watershed management technology. By integrating a method for extracting underlying surface morphology features from multi-source data, combined with a meteorological time-series data feature selection method based on high-dimensional dynamic time distortion distance, and utilizing Sentinel-1 / 2 multi-source remote sensing images to invert flood hydrological process features, the method calculates the geographical morphological similarity, meteorological feature similarity, and hydrological process similarity between data-free and data-rich small watersheds. Weighted comprehensive similarity is then used to complete the small watershed similarity discrimination, achieving high-precision transfer of hydrological model parameters in data-free areas. This invention's method provides comprehensive underlying surface discrimination that closely aligns with hydrological mechanisms, effectively preserving dynamic meteorological and hydrological features and improving time-series classification accuracy. It directly matches flood processes through remote sensing inversion, avoiding misjudgments caused by traditional indirect inference. The process is standardized, adaptable to multiple scenarios, and suitable for hydrological simulation in data-free areas.

Description

Technical Field

[0001] This invention relates to the field of hydrological simulation and watershed management technology, and in particular to a multi-dimensional watershed similarity discrimination method based on geography, meteorology and hydrology. Background Technology

[0002] Watershed similarity assessment is a core technology for solving hydrological simulation in data-scarce areas. Essentially, it achieves parameter transfer or runoff pattern transmission in hydrological models by establishing similarity relationships between data-rich and data-scarce watersheds. The "Data-Scarce Watershed Prediction (PUB)" program, launched by the International Association of Hydrological Sciences (IAHS) in 2003, lists watershed similarity assessment as a core research direction. In areas lacking measured hydrological data, the accuracy of similarity assessment directly determines the reliability of daily runoff predictions for data-scarce watersheds.

[0003] Traditional hydrological simulation relies on measured runoff data to calibrate parameters. However, due to factors such as geographical environment, uneven socio-economic development, and environmental changes, about 60% of watersheds worldwide fall into the category of "no data" or "data-deficient". Watershed similarity discrimination, by exploring the pattern of "similar hydrological responses in similar watersheds", has become a key path to break through the bottleneck of data-deficient area simulation. Its technological evolution has always revolved around how to more accurately quantify watershed similarity. Traditional methods are based on single-dimensional index matching as the core logic and are mainly divided into the following methods: (1) spatial similarity method; (2) physical similarity method. In short, spatial and physical similarity methods are mainly based on the similarity of geographical spatial relationships, climate, and underlying surface conditions between watersheds with and without data. They determine one or more watersheds that are most similar to the watersheds without data and obtain the key parameters of the hydrological model driven by the watersheds without data by directly transplanting known single watershed model parameters, weighted or averaged multi-watershed model parameters. In recent years, in order to obtain model parameters for watersheds with insufficient data, different regional parameterization methods have been proposed at home and abroad. The core idea is to transfer hydrological information from one or more watersheds with data to watersheds with insufficient data.

[0004] Existing watershed similarity discrimination and regional parameterization techniques still have many shortcomings, making it difficult to meet the needs of high-precision hydrological simulation in watersheds without data. On the one hand, traditional similarity discrimination methods (spatial similarity method and physical similarity method) are based on matching indicators in a single dimension or a few dimensions, failing to form a multi-dimensional discrimination system linking "geography-meteorology-hydrology," and thus unable to fully capture the complex influencing factors of watershed hydrological response, resulting in insufficient accuracy of similarity discrimination. On the other hand, existing methods, whether direct parameter transfer, weighted averaging, or various regional parameterization methods, do not fully consider the impact of inherent differences between watersheds, resulting in poor rationality and adaptability of parameter transfer, further restricting the accuracy of hydrological simulation in watersheds without data, and failing to fundamentally break through the technical bottleneck of hydrological simulation in areas without data. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a multi-dimensional watershed similarity discrimination method based on geography, meteorology, and hydrology. This method improves model solution efficiency while ensuring scheduling accuracy, specifically solving the problems of "long time consumption and difficulty in real-time" in traditional model solutions. It provides scientific support for the engineering application of joint optimization scheduling of flood control for watershed water conservancy projects.

[0006] This invention is implemented as follows:

[0007] It includes two aspects, the first aspect:

[0008] A multi-dimensional watershed similarity discrimination method based on geography, meteorology, and hydrology includes the following steps:

[0009] Step S1, Extraction of surface morphology features of small watersheds based on multi-source data:

[0010] Collect national basic geographic information data, land use type maps, soil type maps, soil texture type maps, and small watershed vector layers; preprocess the multi-source raster data and remove invalid data.

[0011] The partitioned statistical method is used to overlay and analyze the preprocessed raster data and the small watershed vector layer, calculate the statistical value of the pixels inside each small watershed polygon, and add the results as new attribute fields to the attribute table of the small watershed layer.

[0012] Extract features of the underlying surface from multiple sources, including small watershed area, average slope, longest runoff path length, longest runoff path gradient, forest area, cultivated land area, and soil type area.

[0013] Geographic feature vectors are constructed based on multi-source underlying surface features. The differences between the feature vectors are calculated using Euclidean distance and converted into geographic morphological similarity. ;

[0014] Step S2, based on high-dimensional dynamic time warp (DTW) M Selection of meteorological time series characteristics based on distance:

[0015] Daily precipitation, daily maximum temperature, daily minimum temperature, solar radiation, and near-surface water vapor pressure data are collected to construct a time-series decision information system (TDS) containing meteorological attributes. The meteorological attributes are standardized using Z-score to eliminate the influence of dimensions and retain the unequal length characteristics of the time-series data.

[0016] A local distance matrix is ​​constructed in a high-dimensional feature space, and dynamic programming is used to recursively calculate the optimal matching path between time series. The cumulative minimum distance on the path is taken as the DTW (Distance-to-Way Distance). M distance;

[0017] Temporal neighborhood relationships are introduced to describe the proximity correlation of meteorological sequences in the time dimension, providing a basis for decision-making.

[0018] A heuristic strategy of forward screening and backward elimination was adopted to select key features by external importance and eliminate redundant features by internal importance, thereby obtaining a subset of core meteorological features.

[0019] Meteorological time-series matrices for small watersheds with and without data are constructed based on core meteorological feature subsets. A high-dimensional dynamic time-warp distance metric is used to measure the similarity between the two time series, yielding the meteorological feature similarity. ;

[0020] Step S3, inversion of flood hydrological process characteristics based on multi-source remote sensing images:

[0021] Multi-source data, including satellite remote sensing data and auxiliary data, are collected. The satellite remote sensing data includes Sentinel-1 SAR data, Sentinel-2 optical data, and Hisa-1 high-resolution SAR data. The auxiliary data includes DEM data and land use data. The collected multi-source remote sensing data and auxiliary data are preprocessed to unify the spatiotemporal reference of the data and standardize the data for subsequent processing.

[0022] For different types of remote sensing data, a differentiated water body extraction index is adopted and combined with threshold segmentation to achieve high-precision water body identification; based on the water body extraction results of multi-temporal remote sensing images, time dimension information is integrated to invert the hydrological process characteristics of floods, including inundation area, duration, start and end time.

[0023] Hydrological response feature vectors for small watersheds are constructed based on the hydrological process characteristics of floods. The differences between these feature vectors are calculated using Euclidean distance and converted into hydrological process similarity. ;

[0024] Step S4, Multi-dimensional Coupling Judgment:

[0025] Overall similarity calculation: ;

[0026] in, As a weight for geographical features, As the weight of hydrological processes, The weights for the three dimensions of meteorological characteristics were determined using the Analytic Hierarchy Process (AHP) combined with expert experience.

[0027] Screening for similar small watersheds: by Sort the data in descending order and select the top 3 watersheds as similar watersheds to the target watershed with no data, thus completing the watershed similarity judgment.

[0028] Furthermore, in step 1, the method for extracting the underlying surface features includes: for continuous raster data, calculating the average value within each watershed to represent the overall topographic features of the watershed; for categorical raster data, calculating the proportion of each type within the watershed, and finally outputting continuous variables of different categories.

[0029] Furthermore, in step 2, the Time Series Decision Information System (TDS) is calculated using the following formula:

[0030]

[0031] In the formula: Small watershed object set;

[0032] : Meteorological attribute set;

[0033] Decision attributes;

[0034] Ordered time sets ;

[0035] : Attribute range, V a This refers to the range of values ​​that attribute 'a' can take.

[0036] Timing information function;

[0037] : Decision function;

[0038] The Z-score standardization of meteorological attributes is used to eliminate the influence of dimensions. The calculation formula is as follows:

[0039]

[0040] In the formula: Let be the mean of attribute 'a'. Let be the standard deviation of attribute a.

[0041] Furthermore, in step 2, constructing the local distance matrix includes:

[0042] For any two small watersheds to be compared for similarity At any given moment and The meteorological attribute vectors are as follows:

[0043]

[0044]

[0045] The Mahalanobis distance between the two is defined as:

[0046]

[0047] In the formula: for The order Mahalanobis covariance matrix, elements It is used to characterize the correlation between meteorological attributes, avoiding the shortcomings of traditional Euclidean distance in ignoring attribute correlation;

[0048] Construct the cumulative distance matrix (CD) between x and y, calculated using the following formula:

[0049]

[0050] In the formula: p is the time series length of x, and q is the time series length of y.

[0051] Furthermore, in step 2, the DTW M Distance, calculated using the following formula:

[0052]

[0053] in, For any two small watersheds to be compared for similarity, Let C be the Mahalanobis distance sub-distance, CD be the two-dimensional cumulative distance matrix, and the boundary conditions be... .

[0054] Furthermore, in step 2, the introduction of temporal neighborhood relationships is based on DTW. M Distance is used to construct temporal neighborhood relationships, forming temporal neighborhood particles that satisfy reflexivity and symmetry and can cover all small watersheds. Based on the temporal neighborhood particles, lower and upper approximations are constructed for decision partitioning to distinguish between small watersheds with definite and uncertain attributions and to describe classification uncertainty. The lower approximation is defined as the positive decision domain, and attribute dependency is calculated to quantify the explanatory power of conditional attributes on decision classification. This is used to realize the approximate representation of small watershed decision classes, classification uncertainty analysis, and key meteorological attribute screening.

[0055] Furthermore, in step 2, the DTW-based M Distance is used to construct temporal neighborhood relationships, and the calculation formula is as follows:

[0056]

[0057] In the formula: For any two small watersheds to be compared for similarity, It is a subset of meteorological attributes. The neighborhood radius;

[0058] The temporal neighborhood particle of the small watershed x under B is calculated using the following formula:

[0059]

[0060] Neighborhood particles satisfy reflexivity. and symmetry, This constitutes the coverage of U;

[0061] The decision division , Let D be the set of small watersheds of class r. Define the lower and upper approximations of D with respect to B, and calculate them as follows:

[0062]

[0063]

[0064] The positive domain is calculated using the following formula:

[0065]

[0066] The dependency ratio is calculated using the following formula:

[0067]

[0068] In the formula: Represents the cardinality of a set. The larger the value, the stronger B's ability to characterize decision D.

[0069] Furthermore, in step 3, the high-precision water body identification is achieved by combining threshold segmentation. The water body difference index SDWI is constructed using the vertical transmit-vertical receive (VV) and vertical transmit-horizontal receive (VH) polarization backscattering coefficients to enhance the difference between the water body and the background. The calculation formula is as follows:

[0070]

[0071] In the formula: The values ​​are calculated for the water body difference index, where VV and VH are the VV and VH polarization backscattering coefficients of the Sentinel-1 image, respectively.

[0072] Using the histogram bimodal method, the valley between the two peaks representing the water body and the background is selected as the water body segmentation threshold. ,when It was determined to be a body of water at that time;

[0073] The difference in reflectance between the green band B3 and the shortwave infrared band B11 is used to highlight water bodies. The calculation formula is as follows:

[0074]

[0075] In the formula: The green band reflectance of Sentinel-2 Reflectivity in the shortwave infrared band;

[0076] Using the histogram bimodal method, when It was determined to be a body of water at that time;

[0077] Since Hisa-1 is single-polarization SAR data, SDWI cannot be calculated. Therefore, the K-means clustering algorithm is used to extract water bodies.

[0078] Initialize cluster centers: Set K=2;

[0079] The Euclidean distance between a pixel and the cluster center is calculated using the following formula:

[0080]

[0081] Where x is the pixel feature vector. Let n be the cluster center of the i-th cluster, and n be the feature dimension;

[0082] The cluster centers are iteratively updated using the following formula:

[0083]

[0084] in, Let i be the set of pixels of the i-th class. Size of the set;

[0085] When the change in cluster centers is less than 0.01 or the number of iterations reaches 100, the low backscattering clusters are classified as water bodies.

[0086] Furthermore, in step 3, the hydrological process characteristics of the inverted flood are as follows:

[0087] Dynamic inversion of flood inundation area: For each temporal phase, water body results are extracted, the number of water body pixels is counted, and the inundation area is calculated by combining the image resolution. The calculation formula is as follows:

[0088]

[0089] In the formula: For the submerged area, The number of pixels representing the water body. Image spatial resolution;

[0090] Calculate the inundated area at each time phase to obtain the dynamic curve of the flood area and identify the flood process;

[0091] Flood duration and inundation timeline inversion: For each pixel, record the time when it was first identified as a body of water. and the last time it was identified as a body of water The flood duration for that pixel is calculated using the following formula:

[0092]

[0093] Based on pixel level A spatial distribution map of flood duration was drawn; combined with ASTER GDEM data, the inundation depth was estimated using the "water level-elevation correlation method," and the mean elevation of the uninundated area was extracted. Assuming the water level in the flooded area The lowest elevation of the unsubmerged area is used to calculate the pixel-level flooding depth, using the following formula:

[0094]

[0095] In the formula: Let x be the flooding depth. is the DEM elevation value of pixel x.

[0096] Furthermore, in step 3, the difference in feature vectors is calculated using Euclidean distance, and the calculation formula is as follows:

[0097]

[0098] in, The maximum inundation area after standardization for watershed T with no data and watershed R with data; For watershed T with no data and watershed R with data, and for the standardized flood duration; Let T be the watershed without data, R be the watershed with data, and the normalized average inundation depth.

[0099] The beneficial effects of this invention are as follows: This invention integrates the three-dimensional features of "geographical morphology—meteorological characteristics—hydrological processes," avoiding the shortcomings of traditional methods that rely solely on a single dimension or static attribute. The underlying surface serves as the basis for hydrological processes, meteorology as the driving factor, and hydrological processes as the response verification. The coupling of these three elements ensures that similar watersheds exhibit "similar attributes, similar driving forces, and similar responses." Compared to existing technologies, this invention has the following significant advantages:

[0100] 1. More comprehensive underlying surface identification, in line with hydrological response mechanisms:

[0101] It covers three dimensions: topography, soil, and vegetation. The weighting is in line with the hydrological mechanism of "topography dominates confluence, soil affects infiltration, and vegetation affects interception". Compared with the traditional physical similarity method, which only involves 2 to 3 indicators, it more comprehensively reflects the overall impact of the underlying surface on the hydrological response.

[0102] 2. Dynamic features are preserved to depict the temporal relationship between meteorology and hydrology:

[0103] It can better handle the problems of "unequal length" and "morphological shift" in meteorological time series, and avoid the problem of losing key dynamic information by using "annual average" in traditional methods. In the validation of public time series datasets such as ArabicDigits and AUSLAN, the meteorological time series classification accuracy is improved by 14.2% to 21.7% compared with the Euclidean distance-based method.

[0104] 3. Compensate for missing hydrological processes and achieve direct similarity matching of flood response processes:

[0105] High-resolution remote sensing inversion directly characterizes the similarity of hydrological processes. Traditional methods completely ignore the similarity of hydrological processes themselves and infer indirectly only through geographical and meteorological static attributes, leading to misjudgments that have similar attributes but large differences in hydrological responses. For example, under the same terrain, the flood inundation processes of short-duration heavy precipitation and long-duration weak precipitation are significantly different.

[0106] 4. Standardized processes, adaptable to multiple application scenarios:

[0107] The entire process of data collection, feature extraction, similarity judgment, and result verification is clearly defined. The parameters of each step can be adapted to different regions through cross-validation. The output format is compatible with hydrological model software and can be directly applied to hydrological simulation in areas without data.

[0108] The present invention will be explained in detail below with reference to the accompanying drawings and specific embodiments. Attached Figure Description

[0109] Figure 1 This is a flowchart of the multidimensional similarity discrimination process of the present invention;

[0110] Figure 2 To create a spatial distribution map of flood duration based on multi-source remote sensing imagery;

[0111] Figure 3 The following is a comparison chart of the simulation results of rainfall-discharge process and parameter migration in small watershed 1 in Example 2:

[0112] Figure 4 Here is a comparison chart of the simulation results of rainfall-discharge processes and parameter migration in small watershed 2 in Example 2:

[0113] Figure 5 Comparison chart of rainfall-discharge process and parameter migration simulation results in small watershed 3 in Example 2:

[0114] Figure 6 Comparison chart of rainfall-discharge process and parameter migration simulation results in small watershed 4 in Example 2:

[0115] Figure 7 Comparison chart of rainfall-discharge process and parameter migration simulation results in small watershed 5 in Example 2:

[0116] Figure 8 Comparison chart of rainfall-discharge process and parameter migration simulation results in small watershed 6 in Example 2:

[0117] Figure 9 Comparison chart of rainfall-discharge process and parameter migration simulation results in small watershed 7 in Example 2:

[0118] Figure 10 Comparison chart of rainfall-discharge process and parameter migration simulation results in small watershed 8 in Example 2:

[0119] Figure 11 Comparison chart of rainfall-discharge process and parameter migration simulation results in small watershed 9 in Example 2:

[0120] Figure 12 This is a comparison chart of the simulation results of rainfall-discharge process and parameter migration in small watershed 10 in Example 2. Detailed Implementation

[0121] Example 1:

[0122] This embodiment presents a multi-dimensional watershed similarity discrimination method based on geography, meteorology, and hydrology, such as... Figure 1 As shown, it includes the following steps:

[0123] Step S1, Extraction of surface morphology features of small watersheds based on multi-source data:

[0124] S11, Data collection and organization of underlying surface:

[0125] Collect national basic geographic information data (including digital elevation model DEM), land use type map, soil type map, soil texture type map, and small watershed vector layer, including attributes such as small watershed area, average slope, longest confluence path length, longest confluence path gradient, forest land ratio, cultivated land ratio, and soil type ratio.

[0126] S12, Data Preparation and Preprocessing:

[0127] Before implementing zonal statistics, rigorous preprocessing of multi-source raster data (excluding small watershed vector layers) is required to ensure spatial consistency, scale uniformity, and attribute purity, eliminating the interference of data heterogeneity on subsequent feature extraction. Specific preprocessing methods are as follows:

[0128] 1) Coordinate System One:

[0129] All raster data must be converted to the same geographic coordinate system or projected coordinate system to ensure consistency in the mathematical basis of spatial overlay analysis. Taking the CGCS2000 national geodetic coordinate system or UTM projection as an example, an affine transformation model is used for projection transformation, and its mathematical expression is as follows:

[0130]

[0131] In the formula: These are the original image pixel coordinates. The coordinates are in the target coordinate system after transformation. For rotation and scaling parameters, These are the translation parameters. The conversion process uses bilinear interpolation or cubic convolution interpolation for pixel resampling to preserve ground feature edge information.

[0132] 2) Uniform spatial resolution:

[0133] Differences in resolution between multi-source data (e.g., DEM 30m, land use 100m) need to be normalized to the target resolution (e.g., 30m) through resampling. Resolution normalization is achieved using the nearest neighbor assignment method (for categorical data) and bilinear interpolation (for continuous data). For categorical data (land use, soil type), it is necessary to ensure that categorical attributes do not produce mixed pixels after resampling; therefore, the nearest neighbor method is used, and the calculation formula is as follows:

[0134]

[0135] In the formula: For the target pixel value, , These are the original resolution and the target resolution, respectively. This is for rounding down.

[0136] For continuous data (DEM, slope), the bilinear interpolation method is used, and the calculation formula is as follows:

[0137]

[0138] In the formula: The values ​​are the values ​​of the four neighboring pixels in the original image. These are the interpolation weights.

[0139] 3) Invalid data removal:

[0140] Raster data may contain null values ​​(NoData) or outliers (such as negative DEM values ​​or abnormal water reflectance), which need to be removed using masking or thresholding methods. The specific steps are as follows:

[0141] Null value handling: Identify pixels marked as NoData in the image and exclude them from subsequent statistics. For continuous data, neighborhood interpolation can be used to fill in the gaps (such as the inverse distance weighting method), but the filling ratio must be controlled to not exceed 5%; otherwise, the small watershed should be directly removed.

[0142] Outlier detection: For continuous data (such as DEM) The principle is to eliminate outliers, i.e., pixel values ​​exceeding the mean. A value more than one standard deviation is considered an anomaly and replaced with the neighborhood mean.

[0143] Categorical data purity check: If the percentage of valid pixels in categorical data in a small watershed is less than 90%, then the watershed is marked as missing for that attribute, and the watershed will be imputed or removed in subsequent analysis.

[0144] The mathematical expression for outlier removal is as follows:

[0145]

[0146] In the formula: The global mean. This represents the global standard deviation. Invalid regions, after being removed, are automatically skipped during partitioned statistics.

[0147] To obtain the overall underlying surface features of each small watershed, a zonal statistical method was adopted. The preprocessed raster data of various types were overlaid and analyzed with the small watershed vector layer. The statistical values ​​of the pixels inside the polygon of each small watershed were calculated, and the results were added as new attribute fields to the attribute table of the small watershed layer.

[0148] The method for extracting features of the underlying surface in a small watershed is as follows:

[0149] 1) For continuous raster data (such as DEM), calculate the average value within each small watershed to represent the overall topographic features of that small watershed.

[0150] 2) For categorical raster data (such as land use and soil type), calculate the proportion of each type within the small watershed and finally output continuous variables of different categories.

[0151] To eliminate the influence of different feature units and numerical ranges, range normalization (Min-MaxNormalization) is used to uniformly map all feature values ​​to the (0,1) interval. For any feature... Its normalized value is calculated as follows:

[0152]

[0153] In the formula: and These are the minimum and maximum values ​​of this feature in all data-rich small watersheds, respectively.

[0154] S13, Geographical morphology similarity extraction:

[0155] Multi-source underlying surface data contain continuous features (such as small watershed area). Average slope Longest Convergence Path Length Longest confluence path ratio ) and continuous variables derived from categorical data (such as the proportion of forest land) arable land ratio (Soil type proportion, etc.). The above p-item core underlying surface characteristics are unified into a geographic feature vector, and the feature vectors for watersheds with and without data are denoted as follows:

[0156]

[0157] In the formula: The core underlying surface features are represented by T, which indicates a small watershed with no data; and R, which indicates a small watershed with data. , These are the i-th underlying surface feature values ​​for the two types of watersheds, respectively.

[0158] The Euclidean distance is used to measure the difference between feature vectors and converted into similarity. First, the Euclidean distance is calculated using the following formula:

[0159]

[0160] In the formula: These are the normalized feature values. The distance is converted to similarity using the following formula:

[0161]

[0162] In the formula: For geographical morphological similarity.

[0163] Step S2, based on high-dimensional dynamic time warp (DTW) M Selection of meteorological time series characteristics based on distance:

[0164] S21, Meteorological Data Acquisition and Preprocessing:

[0165] We collected time-series data on daily precipitation, daily maximum temperature, daily minimum temperature, solar radiation, and near-surface water vapor pressure to construct a time-series decision information system (TDS). We used Z-score standardization to eliminate the influence of dimensions on meteorological attributes and preserved the unequal length characteristics of the time-series data.

[0166] The storage structure for meteorological time-series data is defined as a Time Series Decision Information System (TDS), and the calculation formula is as follows:

[0167]

[0168] In the formula: Small watershed object set;

[0169] : Meteorological attribute set;

[0170] Decision attributes (such as small watershed type labels);

[0171] Ordered time sets ( );

[0172] Attribute range (V) a (The range of values ​​for attribute a).

[0173] Timing information function;

[0174] : Decision function.

[0175] The length of the meteorological sequence may vary due to differences in observation periods, so no truncation or interpolation is needed; it can be achieved through DTW. M Automatic alignment; for each meteorological attribute 'a', Z-score normalization is used to eliminate the influence of dimensions, and the calculation formula is as follows:

[0176]

[0177] In the formula: Let be the mean of attribute 'a'. Let be the standard deviation of attribute a.

[0178] S22, DTW M Distance calculation:

[0179] By constructing a local distance matrix in a high-dimensional feature space, dynamic programming is used to recursively calculate the optimal matching path between time series, and the cumulative minimum distance on the path is taken as the final DTW. M distance.

[0180] For any two small watersheds to be compared for similarity At any given moment and The meteorological attribute vectors are as follows:

[0181]

[0182]

[0183] The Mahalanobis distance between the two is defined as:

[0184]

[0185] In the formula: for The order Mahalanobis covariance matrix, elements It is used to characterize the correlation between meteorological attributes, avoiding the shortcomings of traditional Euclidean distance in ignoring attribute correlation.

[0186] Construct the cumulative distance matrix (CD) between x and y, calculated using the following formula:

[0187]

[0188] In the formula: p is the time series length of x, and q is the time series length of y.

[0189] DTW is calculated recursively using dynamic programming. M Distance, calculated using the following formula:

[0190]

[0191] Boundary conditions: .

[0192] S23, Construction of Temporal Neighborhood Rough Set Model:

[0193] Introducing temporal neighborhood relations to describe the proximity correlation of meteorological sequences in the time dimension, fully exploring the temporal change patterns, improving the accuracy and robustness of similarity matching of meteorological sequences in small watersheds, and avoiding interference from single-point fluctuations on similarity discrimination results;

[0194] Based on DTW M Distance-defined temporal neighborhood relation (DNR) aggregates small watersheds with similar meteorological temporal characteristics into temporal neighborhood particles, characterizing the proximity correlation of small watersheds in the time dimension and attribute space. The calculation formula is as follows:

[0195]

[0196] In the formula: It is a subset of meteorological attributes. is the neighborhood radius.

[0197] The temporal neighborhood particle of watershed x under B is calculated using the following formula:

[0198]

[0199] Neighborhood particles satisfy reflexivity and symmetry This constitutes the coverage of U.

[0200] A decision partitioning mechanism is constructed, introducing lower and upper approximations. These represent the set of watersheds within a neighborhood granularity where all watersheds can be definitively assigned to a certain decision class, and the set of watersheds within a neighborhood granularity where only some watersheds of the same class exist and their assignments cannot be completely determined. This describes the uncertainty in watershed classification based on temporal similarity. The decision partitioning is then discussed. ( Let D be the set of small watersheds of class r). Define the lower and upper approximations of D with respect to B, and calculate them as follows:

[0201]

[0202]

[0203] We further define positive domain and attribute dependency to quantify the explanatory power and importance of meteorological attribute subsets for decision classification.

[0204] The positive domain and dependency of D on B are calculated using the following formula:

[0205]

[0206] In the formula: Represents the cardinality of a set. The larger the value, the stronger B's ability to characterize decision D.

[0207] In real-world scenarios with complex meteorological time-series characteristics and ambiguous small watershed classification boundaries, it can complete approximate representations of small watershed decision classes, analyze classification uncertainties, and evaluate and screen key meteorological attributes.

[0208] S24, Selection of meteorological features:

[0209] A heuristic feature selection strategy combining forward filtering and backward elimination is adopted:

[0210] Meteorological features that contribute significantly to similarity discrimination are added sequentially using a forward filtering method, while features that are highly correlated with the selected features and have redundant information are removed sequentially using a backward elimination method. The feature's ability to distinguish differences between small watersheds is judged based on its external importance; its redundancy within the feature set is judged based on its internal importance; and a core meteorological feature subset is obtained through filtering.

[0211] S25, Meteorological Feature Similarity Extraction:

[0212] The selected core meteorological feature subset includes daily precipitation. Daily maximum temperature Daily minimum temperature Solar radiation Near-surface water vapor pressure Equal time series variables. Meteorological data for each small watershed are multidimensional time series, and their lengths may vary. Small watersheds with no data are recorded. and data for small watersheds The meteorological time series matrix is ​​as follows:

[0213]

[0214] In the formula: The number of meteorological attributes. , These represent the lengths of the time series.

[0215] Employing high-dimensional dynamic time warp distance ( Measuring the similarity between two time series:

[0216]

[0217] In the formula: For meteorological feature similarity, As a scale parameter, it can be taken from all available watersheds. The median distance.

[0218] Step S3, inversion of flood hydrological process characteristics based on multi-source remote sensing images:

[0219] S31, Multi-source remote sensing data acquisition and preprocessing:

[0220] Collect Sentinel-1 SAR data, Sentinel-2 optical data, Hisa-1 high-resolution SAR data, and auxiliary data (DEM, land use).

[0221] 1) Sentinel-1 SAR data preprocessing:

[0222] Orbit correction is performed based on ESA SNAP, automatically updating satellite orbital status information and eliminating orbital drift errors. Thermal noise interference in the imagery is removed; the calculation formula is as follows:

[0223]

[0224] In the formula: The original grayscale value. This represents the grayscale value of thermal noise.

[0225] Convert grayscale values ​​to backscattering coefficients For radiation calibration, the calculation formula is as follows:

[0226]

[0227] In the formula: K is the scaling factor, The angle of incidence is denoted as .

[0228] Lee filtering is used to suppress speckle noise in multi-view filtering. The calculation formula is as follows:

[0229]

[0230] In the formula: For the window mean The variance within the window. This represents the noise variance.

[0231] Based on the WGS84-U TM 49N projection, DEM data is used to eliminate terrain distortion and complete geometric correction.

[0232] 2) Sentinel-2 optical data preprocessing:

[0233] The sensor's DN values ​​are converted to surface reflectance for radiometric calibration. The SMAC (Simplified Method for Atmospheric Correction) model is used to remove the effects of atmospheric scattering and absorption; the formula is as follows:

[0234]

[0235] In the formula: Reflectance of the top layer of the atmosphere, Atmospheric scattering reflectance, Atmospheric transmittance, This is the correction factor for the solar incidence angle.

[0236] Remove cloud-covered areas based on cloud mask files (QA60 band).

[0237] 3) Hisa-1 SAR data preprocessing:

[0238] The original image was resampled to 3m resolution, and Gamma filtering was used to reduce noise. Latitude and longitude coordinates were matched with DEM data and spatially registered with Sentinel-1 / 2 imagery to ensure that the positional error was less than 0.5 pixels.

[0239] S32, Water Extraction:

[0240] Differentiated water body extraction indices are used for different types of remote sensing data, and threshold segmentation is combined to achieve high-precision water body identification.

[0241] The water body difference index SDWI is constructed using the VV (vertical transmit-vertical receive) and VH (vertical transmit-horizontal receive) polarized backscattering coefficients to enhance the difference between the water body and the background. The calculation formula is as follows:

[0242]

[0243] In the formula: The values ​​are calculated for the water body difference index, and VV and VH are the VV and VH polarization backscattering coefficients of the Sentinel-1 image, respectively.

[0244] Using the histogram bimodal method, the valley between the two peaks representing the water body and the background is selected as the water body segmentation threshold. ,when It was determined to be a body of water at that time.

[0245] The difference in reflectance between the green band (B3) and the shortwave infrared band (B11) is used to highlight water bodies. The calculation formula is as follows:

[0246]

[0247] In the formula: The green band reflectance of Sentinel-2 This refers to the reflectivity in the shortwave infrared band.

[0248] Using the histogram bimodal method, when It was determined to be a body of water at that time.

[0249] Since Hisa-1 is single-polarization SAR data, SDWI cannot be calculated. Therefore, the K-means clustering algorithm is used to extract water bodies.

[0250] Initialize cluster centers: Set K=2 (water body, non-water body);

[0251] Calculate the Euclidean distance between the pixel and the cluster center: Where x is the pixel feature vector. Let n be the cluster center of the i-th cluster, and n be the feature dimension.

[0252] Iteratively update cluster centers: ,in, Let i be the set of pixels of the i-th class. The size of the set.

[0253] When the change in cluster centers is less than 0.01 or the number of iterations reaches 100, the low backscattering clusters are ultimately identified as water bodies.

[0254] S33, Inversion of flood hydrological process characteristics:

[0255] Based on the water body extraction results from multi-temporal remote sensing images, and by integrating temporal dimension information, the core hydrological process characteristics of the flood are inverted, including the inundated area, duration, and start and end times.

[0256] 1) Dynamic inversion of flood inundation area:

[0257] For each time phase, the number of water body pixels is counted, and the flooded area is calculated by combining the image resolution. The calculation formula is as follows:

[0258]

[0259] In the formula: For the submerged area, The number of pixels representing the water body. This refers to the spatial resolution of the image.

[0260] The inundation area at each time phase is calculated to obtain the dynamic curve of the flood area and identify the flood process.

[0261] 2) Inversion of flood duration and inundation sequence:

[0262] For each pixel, record the time when it was first identified as a body of water. and the last time it was identified as a body of water The duration of the flood for that pixel is: .

[0263] Based on pixel level Draw a spatial distribution map of the flood duration.

[0264] Using ASTER GDEM data, the inundation depth was estimated using the "water level-elevation correlation method." First, the mean elevation of the uninundated area was extracted. Assuming the water level in the flooded area The lowest elevation of the unsubmerged area is used to calculate the pixel-level flooding depth, using the following formula:

[0265]

[0266] In the formula: Let x be the flooding depth. is the DEM elevation value of pixel x.

[0267] S34, Hydrological Process Similarity Extraction

[0268] Based on the characteristics of each flood hydrological process, a hydrological response feature vector is constructed for each small watershed, including the maximum inundation area. Duration of the flood and average flood depth No data available for small watersheds and data for small watersheds The hydrological feature vector is:

[0269]

[0270] The Euclidean distance is used to measure the difference between feature vectors. First, the features are normalized (same as S12, range normalization method), and the normalized Euclidean distance is calculated using the following formula:

[0271]

[0272] The distance is converted to similarity, and the calculation formula is as follows:

[0273]

[0274] In the formula: This represents the similarity of hydrological processes.

[0275] Step S4, Multi-dimensional Coupling Judgment:

[0276] Using the Analytic Hierarchy Process (AHP) combined with expert experience, the weights of the three dimensions were determined: geographical morphology weight. Hydrological process weights Meteorological feature weights ,satisfy ;

[0277] Overall similarity calculation: ;

[0278] Screening for similar small watersheds: by Sort the data in descending order and select the top N (N<=5) watersheds as similar watersheds to the target watershed with no data, thus completing the watershed similarity judgment.

[0279] Example 2:

[0280] This embodiment is an application example of the multi-dimensional watershed similarity discrimination method based on geography, meteorology, and hydrology, as described in Embodiment 1. The studied watersheds are all located in Shaanxi Province, China, covering areas such as Hanzhong City, Xi'an City, and Shangluo City. Specific details are as follows: The Foping Watershed is located in Foping County, Hanzhong City, with a total area of ​​340 km², and its outlet coordinates are 107°58′44.464800″E, 33°30′11.296800″N. The Laoyukou Watershed is located in Huxian County, Xi'an City, with a total area of ​​321.48 km², and its outlet coordinates are 108.520862E, 33.963853N. The Taipingyu Watershed is located in Huyi District, Xi'an City, with a total area of ​​148.58 km², and its outlet coordinates are 108.678392E, 33.964396N. The Dayu Watershed is located in Chang'an District, Xi'an City, with a total area of ​​55.58 km². The Shipo watershed, located in Chang'an District, Xi'an City, has a total area of ​​259.52 km², with its outlet coordinates at 110.281611 E, 34.001000 N. The Tiesuoguan watershed, located in Ningqiang County, Hanzhong City, has a total area of ​​442.53 km², with its outlet coordinates at 106.410328 E, 32.872310 N. The Yuandun watershed, located in Ningqiang County, Hanzhong City, has a total area of ​​475.94 km², with its outlet coordinates at 106.646550 E, 33.046888 N. The Changtancun watershed, located in Ningqiang County, Hanzhong City, has a total area of ​​219.57 km². The Shipo watershed, located in Chang'an District, Xi'an City, has a total area of ​​259.52 km², with its outlet coordinates at 110.281611 E, 34.228580 N. The Tiesuoguan watershed, located in Ningqiang County, Hanzhong City, has a total area of ​​442.53 km², with its outlet coordinates at 106.410328 E, 32.872310 N. The Majie Basin, located in Majie Town, Shangluo City, covers a total area of ​​315.8 km², with its outlet coordinates at 107.424837 E, 33.301188 N. The Qinghuai Basin, located in Zhen'an County, Shangluo City, covers a total area of ​​166.09 km², with its outlet coordinates at 109.081340 E, 33.442261 N. These basins are situated in the transition zone between the Qinling-Bashan Mountains and the Guanzhong Plain, characterized by complex topography (including mountains, hills, and valleys) and a climate exhibiting characteristics of transition from the northern subtropical to the warm temperate zone. They exhibit significant spatial heterogeneity and seasonal precipitation variations, providing a typical research platform for validating multidimensional watershed similarity discrimination methods.

[0281] Based on the watershed similarity discrimination results, 10 watersheds, including Foping, Laoyukou, and Taipingyu, were classified as pseudo-data-free watersheds from 50 watersheds with available data. Parameters from similar watersheds were transferred to the pseudo-data-free watersheds. Rainfall data from the validation watersheds were input to simulate flood events. The model performance of the transferred parameters was compared to the actual flood simulation, verifying the reliability of the method. The flood simulation results are as follows: Figures 3-12 : .

[0282] As can be seen from the 10 sets of rainfall-discharge process comparison charts, after migrating similar watershed parameters to pseudo-data-free watersheds, the simulated discharge curves and measured discharge curves are highly consistent in the overall trends of the rising water phase, peak time, and receding process. Most of the scenarios can accurately capture the flood peak shape and magnitude, with only minor deviations in peak height and rising and receding rates in some scenarios. Overall, this parameter migration method has high reliability in flood process simulation and can effectively support flood forecasting and risk assessment in data-free watersheds.

[0283] The validation results show that the average Nash efficiency coefficient (NSE) after model parameter migration reaches 0.79, which is 32% higher than that of traditional watershed similarity methods. In summary, the multi-dimensional watershed similarity discrimination method that integrates "geographical morphology, meteorological characteristics and hydrological processes" can effectively perform watershed similarity discrimination. The validation results show that it can achieve high-precision migration of hydrological model parameters in areas without data.

[0284] The above description is only used to illustrate the technical solutions of the present invention and is not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention (such as the application of various formulas, the order of steps, etc.) without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A multi-dimensional watershed similarity discrimination method based on geography, meteorology, and hydrology, characterized in that, Includes the following steps: Step S1, Extraction of surface morphology features of small watersheds based on multi-source data: Collect national basic geographic information data, land use type maps, soil type maps, soil texture type maps and small watershed vector layers, preprocess the above multi-source raster data, remove invalid data, excluding small watershed vector layers; The partitioned statistical method is used to overlay and analyze the preprocessed raster data and the small watershed vector layer, calculate the statistical value of the pixels inside each small watershed polygon, and add the results as new attribute fields to the attribute table of the small watershed layer. Extract features of the underlying surface from multiple sources, including small watershed area, average slope, longest runoff path length, longest runoff path gradient, forest area, cultivated land area, and soil type area. Geographic feature vectors are constructed based on multi-source underlying surface features. The differences between the feature vectors are calculated using Euclidean distance and converted into geographic morphological similarity. ; Step S2, based on high-dimensional dynamic time warp (DTW) M Selection of meteorological time series characteristics based on distance: Daily precipitation, daily maximum temperature, daily minimum temperature, solar radiation, and near-surface water vapor pressure data are collected to construct a time-series decision information system (TDS) containing meteorological attributes. The meteorological attributes are standardized using Z-score to eliminate the influence of dimensions and retain the unequal length characteristics of the time-series data. A local distance matrix is ​​constructed in a high-dimensional feature space, and dynamic programming is used to recursively calculate the optimal matching path between time series. The cumulative minimum distance on the path is taken as the DTW (Distance-to-Way Distance). M distance; Temporal neighborhood relations are introduced to describe the proximity correlation of meteorological sequences in the time dimension, providing a basis for decision-making. A heuristic strategy of forward screening and backward elimination was adopted to select key features by external importance and eliminate redundant features by internal importance, thereby obtaining a subset of core meteorological features. Meteorological time-series matrices for small watersheds with and without data are constructed based on core meteorological feature subsets. A high-dimensional dynamic time-warp distance metric is used to measure the similarity between the two time series, yielding the meteorological feature similarity. ; Step S3, inversion of flood hydrological process characteristics based on multi-source remote sensing images: Multi-source data, including satellite remote sensing data and auxiliary data, are collected. The satellite remote sensing data includes Sentinel-1 SAR data, Sentinel-2 optical data, and Hisa-1 high-resolution SAR data. The auxiliary data includes DEM data and land use data. The collected multi-source remote sensing data and auxiliary data are preprocessed to unify the spatiotemporal reference of the data and standardize the data for subsequent processing. For different types of remote sensing data, a differentiated water body extraction index is adopted and combined with threshold segmentation to achieve high-precision water body identification; based on the water body extraction results of multi-temporal remote sensing images, time dimension information is integrated to invert the hydrological process characteristics of floods, including inundation area, duration, start and end time. Hydrological response feature vectors for small watersheds are constructed based on the hydrological process characteristics of floods. The differences between these feature vectors are calculated using Euclidean distance and converted into hydrological process similarity. ; Step S4, Multi-dimensional Coupling Judgment: Overall similarity calculation: ; in, As a weight for geographical features, As the weight of hydrological processes, The weights for the three dimensions of meteorological characteristics were determined using the Analytic Hierarchy Process (AHP) combined with expert experience. Screening for similar small watersheds: by Sort the data in descending order and select the top 3 watersheds as similar watersheds to the target watershed with no data, thus completing the watershed similarity judgment.

2. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology as described in claim 1, characterized in that, In step 1, the method for extracting the underlying surface features includes: for continuous raster data, calculating the average value within each watershed to represent the overall topographic features of the watershed; for categorical raster data, calculating the proportion of each type within the watershed, and finally outputting continuous variables of different categories.

3. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology as described in claim 1, characterized in that, In step 2, the Time Series Decision Information System (TDS) is calculated using the following formula: In the formula: Small watershed object set; : Meteorological attribute set; Decision attributes; Ordered time sets ; : Attribute range, V a This refers to the range of values ​​that attribute 'a' can take. Timing information function; : Decision function; The Z-score standardization method is used to eliminate the influence of dimensions, and the calculation formula is as follows: In the formula: Let be the mean of attribute 'a'. Let be the standard deviation of attribute a.

4. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology as described in claim 1, characterized in that, Step 2, the construction of the local distance matrix, includes: For any two small watersheds to be compared for similarity At any given moment and The meteorological attribute vectors are as follows: The Mahalanobis distance between the two is defined as: In the formula: for The order Mahalanobis covariance matrix, elements It is used to characterize the correlation between meteorological attributes, avoiding the shortcomings of traditional Euclidean distance in ignoring attribute correlation; Construct the cumulative distance matrix (CD) between x and y, calculated using the following formula: In the formula: p is the time series length of x, and q is the time series length of y.

5. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology as described in claim 1, characterized in that, In step 2, the DTW M Distance, calculated using the following formula: in, For any two small watersheds to be compared for similarity, Let C be the Mahalanobis distance sub-distance, CD be the two-dimensional cumulative distance matrix, and the boundary conditions be... .

6. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology as described in claim 1, characterized in that, In step 2, the introduction of temporal neighborhood relations is based on DTW. M Distance is used to construct temporal neighborhood relationships, forming temporal neighborhood particles that satisfy reflexivity and symmetry and can cover all small watersheds. Based on the temporal neighborhood particles, lower and upper approximations are constructed for decision partitioning to distinguish between small watersheds with definite and uncertain attributions and to describe classification uncertainty. The lower approximation is defined as the positive decision domain, and attribute dependency is calculated to quantify the explanatory power of conditional attributes on decision classification. This is used to realize the approximate representation of small watershed decision classes, classification uncertainty analysis, and key meteorological attribute screening.

7. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology as described in claim 6, characterized in that, In step 2, the DTW-based M Distance is used to construct temporal neighborhood relationships, and the calculation formula is as follows: In the formula: For any two small watersheds to be compared for similarity, It is a subset of meteorological attributes. The neighborhood radius; The temporal neighborhood particle of the small watershed x under B is calculated using the following formula: Neighborhood particles satisfy reflexivity. and symmetry, This constitutes the coverage of U; The decision division , Let D be the set of small watersheds of class r. Define the lower and upper approximations of D with respect to B, and calculate them as follows: The positive domain is calculated using the following formula: The dependency ratio is calculated using the following formula: In the formula: Represents the cardinality of a set. The larger the value, the stronger B's ability to characterize decision D.

8. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology as described in claim 1, characterized in that, In step 3, high-precision water body identification is achieved by combining threshold segmentation. The water body difference index SDWI is constructed by using the vertical transmit-vertical receive (VV) and vertical transmit-horizontal receive (VH) polarization backscattering coefficients to enhance the difference between the water body and the background. The calculation formula is as follows: In the formula: The values ​​are calculated for the water body difference index, where VV and VH are the VV and VH polarization backscattering coefficients of the Sentinel-1 image, respectively. Using the histogram bimodal method, the valley between the two peaks representing the water body and the background is selected as the water body segmentation threshold. ,when It was determined to be a body of water at that time; The difference in reflectance between the green band B3 and the shortwave infrared band B11 is used to highlight water bodies. The calculation formula is as follows: In the formula: The green band reflectance of Sentinel-2, Reflectivity in the shortwave infrared band; Using the histogram bimodal method, when It was determined to be a body of water at that time; Since Hisa-1 is single-polarization SAR data, SDWI cannot be calculated. Therefore, the K-means clustering algorithm is used to extract water bodies. Initialize cluster centers: Set K=2; The Euclidean distance between a pixel and the cluster center is calculated using the following formula: Where x is the pixel feature vector. Let n be the cluster center of the i-th cluster, and n be the feature dimension; The cluster centers are iteratively updated using the following formula: in, Let i be the set of pixels of the i-th class. Size of the set; When the change in cluster centers is less than 0.01 or the number of iterations reaches 100, the low backscattering clusters are classified as water bodies.

9. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology as described in claim 1, characterized in that, In step 3, the hydrological process characteristics of the inverted flood are as follows: Dynamic inversion of flood inundation area: For each temporal phase, water body results are extracted, the number of water body pixels is counted, and the inundation area is calculated by combining the image resolution. The calculation formula is as follows: In the formula: For the submerged area, The number of pixels representing the water body. Image spatial resolution; Calculate the inundated area at each time phase to obtain the flood area dynamic curve and identify the flood process; Flood duration and inundation timeline inversion: For each pixel, record the time when it was first identified as a body of water. and the last time it was identified as a body of water The flood duration for that pixel is calculated using the following formula: Based on pixel level A spatial distribution map of flood duration was drawn; combined with ASTER GDEM data, the inundation depth was estimated using the "water level-elevation correlation method," and the mean elevation of the uninundated area was extracted. Assuming the water level in the flooded area The lowest elevation of the unsubmerged area is used to calculate the pixel-level flooding depth, using the following formula: In the formula: Let x be the flooding depth. is the DEM elevation value of pixel x.

10. The multi-dimensional watershed similarity discrimination method based on geography-meteorology-hydrology according to claim 1, characterized in that, In step 3, the difference in feature vectors is calculated using Euclidean distance, and the calculation formula is as follows: in, The maximum inundation area after standardization for watershed T with no data and watershed R with data; For watershed T with no data and watershed R with data, and for the standardized flood duration; The data represents the basin T (no data), the basin R (data), and the standardized average inundation depth.

Images