A hydrological uncertainty analysis method and system based on probability density envelope

Abstract

This invention discloses a method and system for hydrological uncertainty analysis based on probability density envelopes, relating to the field of hydrology and water resources technology. The method includes: collecting time series data of hydrological elements and preprocessing them; generating a uniform grid and performing kernel density estimation at the grid points to obtain density estimates; resampling each series to obtain the density matrix of the samples; calculating the standard error of the grid points; constructing candidate subsets, prioritizing the minimization of the probability density envelope area and minimizing the difference between upper and lower overprobabilities as secondary objectives, and optimizing the selection of the optimal multiplier; using the optimal multiplier to calculate the lower and upper probability density envelopes at the grid points and calculating the overprobabilities; extracting feature parameters; and comparing and analyzing the uncertainties of hydrological elements at different stations based on the feature parameters. This invention does not require a pre-defined distribution, can adaptively construct the minimum area simultaneous confidence band, and provides multi-dimensional uncertainty indicators, providing a scientific basis for watershed hydrological analysis and risk management.

Description

Technical Field

[0001] This invention relates to the field of hydrology and water resources technology, and in particular to a method and system for hydrological uncertainty analysis based on probability density envelope. Background Technology

[0002] Hydrological processes are influenced by multiple factors, including meteorology, underlying surface conditions, and human activities, exhibiting significant randomness, fuzziness, and nonlinearity. The probability distribution of hydrological elements (such as annual runoff, flood season precipitation, and peak flood level) forms the basis for hydrological frequency analysis, engineering design, and risk assessment.

[0003] Existing methods for analyzing hydrological uncertainty are mainly based on uncertainty analysis of known distributions. Specifically, they assume that hydrological variables follow a certain theoretical distribution, such as the Pearson Type III distribution, the log-normal distribution, or the generalized extreme value distribution. Distribution parameters are estimated using methods such as the method of moments and the maximum likelihood method, and then the probability density function is obtained. This method is computationally simple, but it is prone to model misspecification. When actual data deviates from the assumed distribution, the estimation results will produce systematic biases, especially in the tail region (rare floods or low water levels), where uncertainty is difficult to quantify accurately. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a hydrological uncertainty analysis method based on probability density envelope to solve the problems of unstable cross-station and cross-year comparability caused by the difficulty in real-time connection of reference lines and scales under the introduction of new samples, and the inefficiency of similarity comparison and knowledge reuse due to the lack of dimensionless signatures and indexed retrieval paths for event results.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0007] In a first aspect, the present invention provides a method for hydrological uncertainty analysis based on probability density envelopes, comprising,

[0008] S1. Collect time series of hydrological elements, perform preprocessing, and the preprocessed series is a standardized series;

[0009] S2. Generate a uniform grid for each standardized sequence, and perform kernel density estimation at the grid points to obtain the density estimate;

[0010] S3. Perform Bootstrap resampling on each standardized sequence to obtain the density matrix of the Bootstrap samples;

[0011] S4. Calculate the standard error of the grid points, and the standardized upper and lower biases for each Bootstrap sample;

[0012] S5. Construct candidate subsets and select the optimal lower and upper multipliers so that, under the premise that the coverage is close to the set confidence level P, the area of ​​the probability density envelope is minimized and the difference between the upper and lower hyperprobabilities is minimized.

[0013] S6. Calculate the lower and upper envelopes of the probability density at grid points using the optimal multipliers, and calculate the out-of-probability.

[0014] S7. Extract feature parameters;

[0015] S8. Comparative analysis of the uncertainty of hydrological elements at different stations based on characteristic parameters.

[0016] As a preferred embodiment of the hydrological uncertainty analysis method based on probability density envelope described in this invention, the hydrological elements include one or more of precipitation, runoff, water level, sediment content, and evaporation.

[0017] As a preferred embodiment of the hydrological uncertainty analysis method based on probability density envelope described in this invention, in step S1, time series of hydrological elements from at least two hydrological observation stations are collected, denoted as... K is the total number of stations, n k Let be the length of the sequence at the k-th station; preprocessing includes stationarity testing for each station sequence and Z-score standardization for each sequence. The standardized data is denoted as . , z k,i Let i be the standardized value of the i-th data point in the k-th station sequence, where i = 1 to n. k .

[0018] As a preferred embodiment of the hydrological uncertainty analysis method based on probability density envelope described in this invention, step S3 specifically involves: processing the standardized sequence Z... k Perform B resampling cycles with replacement, sampling n samples each time. k 1 sample, to obtain the Bootstrap sample ensemble For each Bootstrap sample in the same g k,j Kernel density estimation is performed on the above to obtain the Bootstrap density estimate. This forms the Bootstrap density vector: b = 1 ~ B, where b is the Bootstrap sample index, and all Bootstrap density vectors are stacked into a density matrix. , It is a real matrix with B rows and M columns.

[0019] As a preferred embodiment of the hydrological uncertainty analysis method based on probability density envelope described in this invention, wherein: if the standardized upper deviation is obtained through calculation... If it is less than 0, then make The standardized lower bias is 0; if the calculated standardization lower bias is 0. If it is less than 0, then make It is 0.

[0020] As a preferred embodiment of the hydrological uncertainty analysis method based on probability density envelope described in this invention, in step S5, multi-objective optimization is used to determine the optimal multiplier, specifically...

[0021] S501. Take the quantile sequence α to generate a candidate set of upper and lower multipliers, and add the symmetric multiplier c. sym ;,

[0022] S502. For each pair of multipliers, construct the current probability density envelope;

[0023] S503. Traverse all multiplier pairs and calculate the coverage Cov of the current probability density envelope. If |Cov-P|≤ε, where ε is the coverage tolerance, then calculate the area A of the current probability density envelope and the difference between the upper and lower hyperprobabilities D.

[0024] S504. Select the multiplier pair that minimizes the area A of the envelope. If there are multiple multiplier pairs that make the area A equal, select the multiplier with the smallest difference between the upper and lower hyperprobabilities D as the optimal solution. If there is no solution that satisfies the coverage tolerance, take the symmetric multiplier as the optimal solution.

[0025] As a preferred embodiment of the hydrological uncertainty analysis method based on probability density envelope described in this invention, the extraction of feature parameters includes extracting the average relative width. Extracting the maximum relative width Extract the tail width ratio Extracting asymmetric coefficients .

[0026] As a preferred embodiment of the hydrological uncertainty analysis method based on probability density envelope described in this invention, wherein: Among them, g j For the j-th grid point, For grid point g j The final probability density envelope value on the surface. The envelope value of the final probability density. For grid point g j The kernel density estimate on, when When it is 0, then the corresponding Items are not included in the summation;

[0027] ;

[0028] Where W(z) represents the width interpolation function at z, obtained through linear interpolation from {g i}and We get z 0.05 and z 0.95 These are the 0.05 and 0.95 quantiles of the standardized data, respectively.

[0029] When the denominator is zero, the corresponding term Take 0.

[0030] As a preferred embodiment of the hydrological uncertainty analysis method based on probability density envelope described in this invention, in step S9, during uncertainty analysis, the average relative width is compared to determine the overall level of uncertainty; the maximum relative width and its position are compared to identify high-risk intervals; the tail width ratio is compared to analyze the asymmetry of uncertainty between low-value and high-value areas; the asymmetry coefficient is compared to determine the direction of systematic deviation in the estimation; and the balance of the estimation is evaluated by combining the upper and lower overprobabilities.

[0031] Secondly, the present invention provides a hydrological uncertainty analysis system based on probability density envelope, and a hydrological uncertainty analysis method based on probability density envelope as described in any one of claims 1 to 9, characterized in that it includes:

[0032] The data acquisition and preprocessing module is used to acquire and preprocess the time series of hydrological elements from at least two hydrological observation stations.

[0033] The kernel density estimation module is used to generate a uniform grid for each standardized sequence and perform kernel density estimation at the grid points to obtain the density estimate.

[0034] The Bootstrap resampling module is used to generate Bootstrap samples and density matrices;

[0035] The sample calculation module is used to calculate the standard error at grid points, and the standardized upper and lower biases for each Bootstrap sample;

[0036] The optimization module is used to determine the optimal multiplier;

[0037] The optimal multiplier calculation module is used to calculate the lower envelope, upper envelope, and excess probability of the probability density at grid points based on the optimal multiplier.

[0038] The feature extraction module calculates feature parameters;

[0039] The uncertainty analysis module is used to compare and analyze the uncertainty of hydrological elements at different stations based on the calculated characteristic parameters.

[0040] The beneficial effects of this invention are as follows: This invention eliminates the need for pre-setting the distribution form, adaptively fitting various hydrological element data through kernel density estimation, thus avoiding systematic biases caused by model misconfiguration in parametric methods. Simultaneously, the use of rotation rules to automatically determine the bandwidth improves the automation and applicability of the method. This invention constructs a simultaneous confidence band rather than a point-by-point confidence interval, obtaining the density estimation distribution through bootstrap resampling and determining the multiplier based on the quantile of the maximum deviation, thereby ensuring that the entire density curve is covered by the probability density envelope with a set probability P, overcoming the defect that point-by-point intervals cannot control the overall coverage probability. A multi-objective optimization method is introduced, with minimizing the area of ​​the probability density envelope as the primary objective and minimizing the difference between the upper and lower hyperprobabilities as the secondary objective, searching while satisfying the coverage tolerance. The optimal multiplier is selected, which makes the obtained probability density envelope as compact as possible while ensuring confidence, improving estimation efficiency and effectively controlling the systematic bias of the estimation. Four characteristic parameters are extracted: average relative width, maximum relative width, tail width ratio, and asymmetry coefficient. These parameters comprehensively characterize the uncertainty of density estimation from multiple dimensions, including overall uncertainty, local high-risk intervals, tail asymmetry, and systematic bias, providing a quantitative tool for multi-site comparison. These parameters all have clear physical or statistical significance, making them easy for hydrological workers to understand and apply. It is applicable to various hydrological elements such as precipitation, runoff, and water level, and has strong versatility and scalability. The multi-site comparison analysis function can provide a scientific basis for watershed hydrological monitoring network optimization, climate change impact assessment, and integrated water resources management. Attached Figure Description

[0041] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0042] Figure 1 This is the overall flowchart of the present invention.

[0043] Figure 2 This is a graph showing the annual precipitation process at Kaifeng Station, Suxian Station, and Huoshan Station in the Huai River Basin from 1960 to 2020.

[0044] Figure 3 This is the probability density envelope diagram of the annual precipitation at Kaifeng station after standardization.

[0045] Figure 4 This is the probability density envelope diagram of the annual precipitation at Suxian Station after standardization.

[0046] Figure 5 This is the probability density envelope diagram of the annual precipitation at Huoshan Station after standardization.

[0047] Figure 6 A radar comparison chart of precipitation characteristic parameters at three stations in the Huai River Basin. Detailed Implementation

[0048] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0049] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0050] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0051] Reference Figure 1 As one embodiment of the present invention, this embodiment provides a hydrological uncertainty analysis method based on probability density envelope, comprising the following steps:

[0052] S1. Collect time series of hydrological elements from at least two hydrological observation stations, denoted as... K is the total number of stations, n k Let n be the length of the sequence (sample) at the k-th station. k ≥20; Preprocess the collected data, including,

[0053] For each station sequence, if a trend or periodicity is found, the stationarity test (such as the ADF test) is performed to make it satisfy the independent and identically distributed or weakly stationary assumptions through differencing or detrending.

[0054] Z-score normalization is performed on each sequence. The standardized data of the sample are denoted as , z k,i Let i be the standardized value of the i-th data point in the k-th station sequence, where i = 1 to n. k μ k Let σ be the mean of the sequence of the k-th station. k Let be the sample standard deviation of the sequence at the k-th station;

[0055] Hydrological elements include one or more of the following: precipitation, runoff, water level, sediment content, and evaporation.

[0056] S2. Generate a uniform grid for each standardized sequence. j=1~M; the number of grid points M is preferably 200;

[0057] Kernel density estimation is performed at grid points. ;

[0058] in, Let J be the coordinates of the j-th grid point at the k-th station (dimensionless). For grid points The probability density function is given by , where min represents the minimum value and max represents the maximum value. , h k For bandwidth, the rotation rule is used to automatically determine it: , For Z k The sample standard deviation; let the density estimates at the grid points be vectors: .

[0059] S3, for the standardized sequence Z k Perform B resampling cycles with replacement, sampling n samples each time. k 1 sample, to obtain the Bootstrap sample ensemble For each Bootstrap sample in the same g k,j Kernel density estimation is performed on the above to obtain the Bootstrap density estimate. This forms the Bootstrap density vector: b=1~B, stack all Bootstrap density vectors into a density matrix. , It is a 1-row, M-column real number matrix, where b is the Bootstrap sample index. It is a real matrix with B rows and M columns.

[0060] S4. Calculate the standard error of the grid points: , Bootstrap sample at grid points The average value is calculated; the standard deviation of each Bootstrap sample is calculated. ; Calculate the standardized bias for each Bootstrap sample: If the standard deviation is calculated using the above formula... If it is less than 0, then make The standard deviation is 0; if the standard deviation calculated using the above formula is... If it is less than 0, then make It is 0.

[0061] S5. Construct candidate subsets and select the optimal lower and upper multipliers such that, while ensuring the coverage is close to the set confidence level P, the area of ​​the probability density envelope is minimized and the difference between the upper and lower hyperprobabilities is minimized. Specifically:

[0062] S501, Obtain the quantile sequence ,calculate , c down c is the submultiplier. up For the superior component, among which, Standardized upper bias d for Bootstrap samples up sequential Quantiles Standardized lower bias d for Bootstrap samples down sequential Quantiles, calculating symmetric multipliers: P is the set confidence level (e.g., 0.90 or 0.95). The standardized maximum deviation (d) of the Bootstrap sample up With d down The P-quantile of the sequence formed by the maximum value in (c) is used to calculate all the above c. up (α), c sym Merge into candidate superior subset C up , and all of the above c down (α), c sym Merge into candidate lower product subset C down ;

[0063] S502, for each pair (c down , c up )∈C down ×C up Construct the current probability density envelope.

[0064] ,

[0065] ;

[0066] The envelope value is the probability density function. The envelope value on the probability density;

[0067] S503, Calculate the coverage of the Bootstrap samples. , Where I{∙} is an indicator function, which takes the value 1 when the condition within the parentheses is true, and 0 otherwise; if |Cov-P|≤ε, where ε is the allowable coverage deviation tolerance, ε is preferably 0.02, calculate the envelope area A and the difference between the upper and lower hyperprobabilities D. ;

[0068] ;

[0069] S504. Select the multiplier pair that minimizes the area A of the envelope. If multiple multiplier pairs exist that result in equal areas A, then select the multiplier pair with the smallest upper and lower hyperprobability difference D as the optimal solution. If no solution satisfies the coverage tolerance, take the symmetric multiplier (c). sym , c sym As the optimal solution, that is, let =(c sym , c sym ).

[0070] S6. Calculate the lower envelope of the probability density at the grid points using the optimal multipliers, specifically:

[0071] ;

[0072] The envelope of the probability density at grid points is calculated using the optimal multipliers, specifically as follows: ;

[0073] Recalculate the out-of-range case for the Bootstrap sample, and you will get...

[0074] Lower probability for, ;

[0075] Overprobability for, ;

[0076] High probability for, , The value after this must not exceed the preset coverage tolerance ε. In actual calculations, when... If the value exceeds ε, increase the number of Bootstrap resampling times B in step S3 (for example, gradually increase B by multiples of the original value), and repeat steps S3 to S6 until the coverage tolerance requirement is met.

[0077] S7. Extract feature parameters, including:

[0078] Extracting the average relative width , Among them, g j For the j-th grid point, For grid point g j The final probability density envelope value on the surface. The envelope value of the final probability density. For grid point g j The kernel density estimate on, when When it is 0, then the corresponding Items are not included in the summation;

[0079] Extract the maximum relative width , ;

[0080] Extract tail width ratio , Where W(z) represents the width interpolation function at z, obtained through linear interpolation from {g i}and We get z 0.05 and z 0.95 These are the 0.05 and 0.95 quantiles of the standardized data, respectively.

[0081] Extracting asymmetric coefficients , When the denominator is zero, the corresponding term Take 0.

[0082] S8. For the z-axis corresponding to the original data points k,i Using the same kernel function and bandwidth as in step S2, the density estimate at that point is calculated. , ;

[0083] The standard error at that point is obtained through linear interpolation. ,

[0084] ;

[0085] Where Lerp(∙) represents based on grid points {g j} and corresponding value For z k,j Interpolation operators for linear interpolation (extrapolation);

[0086] Calculate the probability density envelope value at this point. , , .

[0087] S9. Based on characteristic parameters, a comparative analysis of the uncertainty of hydrological elements at different stations is conducted, specifically comparing the average relative width. To determine the overall level of uncertainty; to compare the maximum relative width. and its position g max Identify high-risk areas; compare tail width ratios Analyze the asymmetry of uncertainty between low-value and high-value regions; compare the asymmetry coefficients. Determine the direction of the systematic bias in the estimate; combine the upper and lower overprobabilities. , To assess the balance of the estimates.

[0088] Between steps S8 and S9, the following steps are also included:

[0089] (1) Output calculation results: Output the result matrix for each station. , for A real matrix with 5 rows and 5 columns, whose column vectors are as follows: Original data X k Standardized data Z k Density estimates Lower envelope value Upper envelope value Optimal multiplier vector: , A real number vector containing 2 elements, exceeding the probability vector: , For a real vector containing 3 elements, the eigenvector is: , It is a real number vector containing 4 elements.

[0090] (2) Plot the standardized data Z k The histogram (normalized to probability density form) is plotted, and the following curves / regions are overlaid on the graph: kernel density estimation curves. Envelope under probability density probability density envelope And the confidence band filling area.

[0091] This embodiment also provides a hydrological uncertainty analysis system based on probability density envelope, including:

[0092] The data acquisition and preprocessing module is used to acquire and preprocess the time series of hydrological elements from at least two hydrological observation stations.

[0093] The kernel density estimation module is used to generate a uniform grid for each standardized sequence and perform kernel density estimation at the grid points to obtain the density estimate.

[0094] The Bootstrap resampling module is used to generate Bootstrap samples and density matrices;

[0095] The sample calculation module is used to calculate the standard error at grid points, and the standardized upper and lower biases for each Bootstrap sample;

[0096] The optimization module is used to determine the optimal multiplier;

[0097] The optimal multiplier calculation module is used to calculate the lower envelope, upper envelope, and excess probability of the probability density at grid points based on the optimal multiplier.

[0098] The feature extraction module calculates feature parameters;

[0099] The uncertainty analysis module is used to compare and analyze the uncertainty of hydrological elements at different stations based on the calculated characteristic parameters.

[0100] This invention eliminates the need for pre-defined distribution patterns, adaptively fitting various hydrological element data through kernel density estimation, thus avoiding systematic biases caused by model misconfiguration in parametric methods. Simultaneously, it employs rotation rules to automatically determine bandwidth, improving the method's automation and applicability. This invention constructs a simultaneous confidence band rather than point-by-point confidence intervals, obtaining the density estimation distribution through bootstrap resampling and determining the multiplier based on the quantile of the maximum deviation. This ensures the entire density curve is covered by the probability density envelope with a set probability P, overcoming the limitation of point-by-point intervals in controlling the overall coverage probability. A multi-objective optimization method is introduced, prioritizing minimizing the area of ​​the probability density envelope and secondarily minimizing the difference between upper and lower hyperprobabilities, searching for the optimal multiplier while satisfying coverage tolerance. This approach ensures that the obtained probability density envelope is as compact as possible while maintaining confidence, improving estimation efficiency and effectively controlling systematic bias. Four characteristic parameters—mean relative width, maximum relative width, tail width ratio, and asymmetry coefficient—are extracted to comprehensively characterize the uncertainty of density estimation from multiple dimensions, including overall uncertainty, local high-risk intervals, tail asymmetry, and systematic bias. This provides a quantitative tool for multi-site comparisons. These parameters all have clear physical or statistical significance, facilitating understanding and application by hydrological workers. Applicable to various hydrological elements such as precipitation, runoff, and water level, it possesses strong versatility and scalability. Its multi-site comparison analysis function can provide a scientific basis for optimizing watershed hydrological monitoring networks, assessing the impact of climate change, and comprehensively managing water resources.

[0101] It should be noted that performing stationarity tests on each sequence is existing technology (the specific methods for stationarity testing are not elaborated in this application, as they are not improvements of this application), such as the Augmented Dickey-Fuller test (ADF), a commonly used method in time series analysis to determine whether a sequence is stationary. If the test finds that the sequence has a trend or periodicity, conventional differencing methods or detrending processing (such as linear regression to remove the trend) are used to convert it into a stationary sequence to meet the basic requirements of independent and identically distributed or weakly stationary data for subsequent bootstrap resampling. The above-mentioned stationarity tests and preprocessing methods are common knowledge to those skilled in the art, and for details, please refer to classic time series analysis literature, such as "Time Series Analysis: Forecasting and Control" (ISBN-978-0-470-27284-8) edited by George E. Box et al.

[0102] Example 2: Refer to Figures 2-6 This is the second embodiment of the present invention. The difference between this embodiment and embodiment 1 is that this embodiment uses a specific case to conduct hydrological uncertainty analysis. Specifically, uncertainty analysis is conducted using the annual precipitation of three rain gauge stations in the Huai River Basin.

[0103] Annual precipitation data (unit: mm) for 61 years (1960–2020) were selected from three rain gauge stations at different latitudes in the Huai River basin (Kaifeng station, 34.78°N; Suxian station, 33.63°N; Huoshan station, 31.40°N). The annual precipitation process at the three stations is as follows: Figure 2 As shown, the multi-year average rainfall at Kaifeng, Suxian, and Huoshan stations is 616 mm, 860 mm, and 1371 mm, respectively. The three annual rainfall series showed no significant trend after ADF testing, but all exhibited strong autocorrelation. Block bootstrap resampling was used (block length 6 for Kaifeng station, and block length 7 for both Suxian and Huoshan stations) to ensure independence.

[0104] Following step S1, Z-score standardization was performed on the annual precipitation series of Kaifeng Station, Suxian Station, and Huoshan Station to obtain their respective standardized data. The confidence level Р=0.95, Bootstrap times B=1000, grid points M=200, and coverage tolerance ε=0.02 were set. The method of this invention was applied to obtain the output and visualization results for each station.

[0105] (1) Kaifeng Station

[0106] The resulting matrix M of the annual precipitation sequence at this station pdAs shown in Table 1, the optimal multiplier vector is C = (2.89, 2.75), where the lower multiplier... =2.89, multiplier =2.75. The probability vector Pe = (0.040, 0.060, 0.032), where the lower probability is... =0.040, overshoot probability =0.060, and also exceeding the probability. =0.032. Eigenvector V c =(1.045,2.029, 1.201, -0.024) T The average relative width =1.045, maximum relative width =2.029 (appears near z=2.41), tail width ratio =1.201, asymmetry coefficient =-0.024. |(p blow +p above -p both The formula )-(1-P)|= |(0.040+0.060-0.032)-(1-0.95)|=0.018<ε indicates that the calculation results of the method of the present invention at Kaifeng station are reasonable. The standardized probability density envelope diagram of annual precipitation at Kaifeng station from 1960 to 2020 is shown below. Figure 3 As shown, Figure 3 It includes histograms, kernel density estimation curves, lower envelope of probability density, upper envelope of probability density, and confidence bands at pre-set confidence levels.

[0107] Table 1. M of the precipitation sequence at Kaifeng Station from 1960 to 2020 pd data

[0108]

[0109] (2) Suxian Station

[0110] The resulting matrix M of the annual precipitation sequence at this station pd As shown in Table 2, the optimal multiplier vector C = (2.69, 2.98), the out-of-probability vector Pe = (0.040, 0.040, 0.016), and the eigenvector V c =(1.338, 2.825, 1.098, 0.080) T The maximum relative width occurs around z=2.95. |(p blow +p above -p bothThe formula )-(1-P)|=|(0.040+0.040-0.016)-(1-0.95)|=0.014<ε indicates that the calculation results of the method of the present invention at the Suxian station are reasonable. The standardized probability density envelope diagram of annual precipitation at the Suxian station from 1960 to 2020 is shown below. Figure 4 As shown, Figure 4 It includes histograms, kernel density estimation curves, lower envelope of probability density, upper envelope of probability density, and confidence bands at pre-set confidence levels.

[0111] Table 2. M of the precipitation sequence at Suxian Station from 1960 to 2020 pd data

[0112]

[0113] (3) Huoshan Station

[0114] The resulting matrix M of the annual precipitation sequence at this station pd As shown in Table 3, the optimal multiplier vector C = (2.73, 2.97), the out-of-probability vector Pe = (0.020, 0.042, 0.006), and the eigenvector V c =(1.097, 2.129, 1.486, 0.042) T The maximum relative width occurs around z=2.60. |(p blow +p above -p both The formula )-(1-P)|=|(0.020+0.042-0.006)-(1-0.95)|=0.006<ε indicates that the calculation results of the method of the present invention at Huoshan Station are reasonable. The standardized probability density envelope diagram of annual precipitation at Huoshan Station from 1960 to 2020 is shown below. Figure 5 As shown, Figure 5 It includes histograms, kernel density estimation curves, lower envelope of probability density, upper envelope of probability density, and confidence bands at pre-set confidence levels.

[0115] Table 3. M of the precipitation sequence at Huoshan Station from 1960 to 2020 pd data

[0116]

[0117] (4) Comparative analysis

[0118] Based on the four feature parameters extracted by the method of this invention and the excess probability vector, a systematic comparative analysis of the uncertainty of annual precipitation at the three stations was conducted, and the results are summarized in Table 4.

[0119] Table 4 Summary of characteristic parameters and exceedance probability for each station

[0120]

[0121] ① Overall Uncertainty Analysis

[0122] Average relative width This reflects the average relative uncertainty of the density estimate across the entire range of values. (Suxian Station) The highest value is (1.338), while Kaifeng Station (1.045) and Huoshan Station (1.097) are relatively close and slightly lower. This indicates that the overall estimation uncertainty of annual precipitation at Suxian Station is the highest, and its probability density curve fluctuates more significantly. Combined with geographical location analysis, Suxian Station is located in the southern part of the Huang-Huai-Hai Plain, in the transitional zone between northern and southern climates, resulting in greater precipitation variability and therefore higher uncertainty. Kaifeng Station is located in the eastern Henan Plain, where precipitation is relatively stable. Although Huoshan Station is located in the Dabie Mountains, it receives abundant precipitation, but its uncertainty falls between the two due to topographical influences.

[0123] ② Identification of high-risk areas

[0124] Maximum relative width The location of these values ​​was used to identify the most unstable local regions, i.e., high-risk areas. The maximum relative widths of the three stations all occurred in regions with standardized values ​​z > 2.4, corresponding to wet years or extreme precipitation events. Among them, the maximum relative width of the Suxian station was... The maximum relative width at (2.825) occurs around z=2.95, indicating that the density estimation at this station in the extreme high-value area is extremely unstable, and there is considerable uncertainty in the probability estimation of extreme precipitation. The maximum relative widths at Huoshan Station (2.129) and Kaifeng Station (2.029) are relatively small, but the maximum relative width at Huoshan Station is... The occurrence of this value near z=2.60 indicates that the uncertainty in its extreme high-value area is also quite prominent. This result reveals that the high-value tail should be the focus of attention in the extreme precipitation frequency analysis of the three stations, with the risk being particularly significant at the Suxian station.

[0125] ③ Distribution symmetry and tail characteristics

[0126] Tail width ratio This measures the degree of asymmetry in uncertainty between low-water-rate areas (dry years) and high-water-rate areas (wet years). (Huoshan Station) (1.486) is greater than 1, indicating that the uncertainty in the low-value area is significantly greater than that in the high-value area, that is, the density estimate of precipitation in a dry year is more uncertain than that in a wet year; Kaifeng station (1.201) is greater than 1, also showing a slightly larger uncertainty in the low-value area, but to a lesser extent; Suxian Station (1.098) is close to 1, indicating that its tail uncertainty is relatively symmetrical, with small differences in uncertainty between low-value and high-value areas. This characteristic is related to the climate background of each station. Specifically, Huoshan station has abundant precipitation and relatively few samples in dry years, leading to unstable estimation in the low-value area; while Suxian station has frequent extreme precipitation events, with more samples in both high-value and low-value areas, resulting in a more balanced tail uncertainty.

[0127] asymmetric coefficient Reflects the direction of systematic bias in density estimation. (Suxian Station) (0.080) is greater than 0, indicating that the upper envelope of the probability density deviates slightly more than the lower envelope of the probability density, that is, the density estimation has a slight upward bias; Huoshan Station (0.042) is greater than 0, also showing a slight upward bias; Kaifeng Station (-0.024) is less than 0, indicating a slight downward bias. Although the biases at all three stations are small (absolute value <0.1), the directional differences are still noteworthy. The upward bias at Suxian and Huoshan stations may be related to a slightly higher estimate of the frequency of extreme high-value events, while the downward bias at Kaifeng station shows the opposite trend. These subtle biases can reveal systematic differences in precipitation distribution across different regions in multi-station comparisons.

[0128] ④ Exceeding probability and estimation equilibrium

[0129] The out-of-bounds probability vector Pe describes the probability distribution of the density curve exceeding the upper and lower bounds of the probability density envelope. At Kaifeng station, the probability of exceeding the upper bound (0.060) is greater than the probability of exceeding the lower bound (0.040), seemingly contradicting the negative value of the asymmetry coefficient (-0.024). However, it should be noted that the asymmetry coefficient is a point-by-point average, while the out-of-bounds probability is the overall event probability; they reflect characteristics at different levels. At Suxian station, the probabilities of exceeding the upper and lower bounds are equal (both 0.040), which does not perfectly correspond to the positive asymmetry coefficient (0.080). This indicates that although the overall probability density envelope deviates significantly, the proportion of samples experiencing exceeding the upper bound is equal to that of experiencing exceeding the lower bound, suggesting a relatively balanced estimation. At Huoshan station, the probability of exceeding the upper bound (0.042) is greater than the probability of exceeding the lower bound (0.020), consistent with the positive direction of the asymmetry coefficient (0.042), indicating that its upward bias is consistent with the increasing probability of exceeding the upper bound. Meanwhile, the out-of-bounds probability P... both The values ​​are relatively small (0.006~0.032) across the three sites, indicating that the probability density envelope is likely to exceed both the upper and lower bounds simultaneously, suggesting that the design of the probability density envelope is reasonable.

[0130] ⑤ Comprehensive comparison

[0131] Figure 6The radar chart visually illustrates the relative differences of the three stations in four characteristic parameters: Suxian Station has a clear advantage in average relative width and maximum relative width, Huoshan Station is most prominent in tail width ratio, and Kaifeng Station has relatively balanced parameters.

[0132] In summary, the main conclusions are as follows:

[0133] a) The overall uncertainty of Kaifeng Station is moderate. The high-risk area is located in a wet year, the uncertainty of the low-value area is slightly higher than that of the high-value area, the density estimate is slightly downward biased, the probability of overshoot is slightly higher than the probability of undershoot, and the overall estimate is relatively stable.

[0134] b) Suxian Station has the highest overall uncertainty, the most prominent risk in the extreme high value area, relatively symmetrical uncertainty at the tail, slightly upward bias in density estimation, and balanced upper and lower probabilities, making it suitable as a key station for extreme precipitation frequency analysis.

[0135] c) The overall uncertainty of Huoshan Station is moderate to low, but the uncertainty in the high-value area is significantly greater than that in the high-value area. There is also a certain risk in the extremely high-value area. The density estimate is slightly upward biased, and the probability of exceeding the upper limit is significantly greater than that of exceeding the lower limit. We need to be wary of the uncertainty in the estimate for dry years.

[0136] The method of this invention successfully quantifies the uncertainty characteristics of annual precipitation at three stations by extracting multidimensional feature parameters, revealing the inherent differences in precipitation distribution under different climatic and geographical backgrounds, and providing a scientific basis for watershed water resources planning, extreme event risk assessment and monitoring station optimization.

[0137] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are 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 without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for analyzing hydrological uncertainties based on probability density envelopes, characterized in that: include, S1. Collect time series of hydrological elements, perform preprocessing, and the preprocessed series is a standardized series; S2. Generate a uniform grid for each standardized sequence, and perform kernel density estimation at the grid points to obtain the density estimate; S3. Perform Bootstrap resampling on each standardized sequence to obtain the density matrix of the Bootstrap samples; S4. Calculate the standard error of the grid points, and the standardized upper and lower biases for each Bootstrap sample; S5. Construct candidate subsets and select the optimal lower and upper multipliers so that, under the premise that the coverage is close to the set confidence level P, the area of ​​the probability density envelope is minimized and the difference between the upper and lower hyperprobabilities is minimized. S6. Calculate the lower and upper envelopes of the probability density at grid points using the optimal multipliers, and calculate the out-of-probability. S7. Extract feature parameters; S8. Comparative analysis of the uncertainty of hydrological elements at different stations based on characteristic parameters.

2. The hydrological uncertainty analysis method based on probability density envelope as described in claim 1, characterized in that: The hydrological elements include one or more of the following: precipitation, runoff, water level, sediment content, and evaporation.

3. The hydrological uncertainty analysis method based on probability density envelope as described in claim 1, characterized in that: In step S1, time series of hydrological elements are collected from at least two hydrological observation stations, denoted as... K is the total number of stations, n k Let be the length of the sequence at the k-th station; preprocessing includes stationarity testing for each station sequence and Z-score standardization for each sequence. The standardized data is denoted as . , z k,i Let i be the standardized value of the i-th data point in the k-th station sequence, where i = 1 to n. k .

4. The hydrological uncertainty analysis method based on probability density envelope as described in claim 3, characterized in that: Step S3 specifically involves processing the standardized sequence Z. k Perform B resampling cycles with replacement, sampling n samples each time. k 1 sample, to obtain the Bootstrap sample ensemble For each Bootstrap sample in the same g k,j Kernel density estimation is performed on the above to obtain the Bootstrap density estimate. This forms the Bootstrap density vector: b = 1 ~ B, where b is the Bootstrap sample index, and all Bootstrap density vectors are stacked into a density matrix. , It is a real matrix with B rows and M columns.

5. The hydrological uncertainty analysis method based on probability density envelope as described in claim 1, characterized in that: If the standardized upper deviation is obtained through calculation If it is less than 0, then make The standardized deviation is 0; if the calculated deviation is... If it is less than 0, then make It is 0.

6. The hydrological uncertainty analysis method based on probability density envelope as described in claim 5, characterized in that: In step S5, multi-objective optimization is used to determine the optimal multiplier, specifically as follows: S501. Take the quantile sequence α to generate a candidate set of upper and lower multipliers, and add the symmetric multiplier c. sym ;, S502. For each pair of multipliers, construct the current probability density envelope; S503. Traverse all multiplier pairs and calculate the coverage Cov of the current probability density envelope. If |Cov-P|≤ε, where ε is the coverage tolerance, then calculate the area A of the current probability density envelope and the difference between the upper and lower hyperprobabilities D. S504. Select the multiplier pair that minimizes the area A of the envelope. If there are multiple multiplier pairs that make the area A equal, select the multiplier with the smallest difference between the upper and lower hyperprobabilities D as the optimal solution. If there is no solution that satisfies the coverage tolerance, take the symmetric multiplier as the optimal solution.

7. The hydrological uncertainty analysis method based on probability density envelope as described in claim 1, characterized in that, Feature parameter extraction includes extracting the average relative width. Extracting the maximum relative width Extract the tail width ratio Extracting asymmetric coefficients .

8. The hydrological uncertainty analysis method based on probability density envelope as described in claim 7, characterized in that: Among them, g j For the j-th grid point, For grid point g j The final probability density envelope value on the surface. The envelope value of the final probability density. For grid point g j The kernel density estimate on the above, when When it is 0, then the corresponding Items are not included in the summation; ； Where W(z) represents the width interpolation function at z, obtained through linear interpolation from {g i }and We get z 0.05 and z 0.95 These are the 0.05 and 0.95 quantiles of the standardized data, respectively. When the denominator is zero, the corresponding term Take 0.

9. The hydrological uncertainty analysis method based on probability density envelope as described in claim 8, characterized in that: In step S9, during uncertainty analysis, the average relative width is compared to determine the overall level of uncertainty; the maximum relative width and its location are compared to identify high-risk areas. By comparing the tail width ratio, we can analyze the asymmetry of uncertainty between the low-value and high-value regions; by comparing the asymmetry coefficients, we can determine the direction of the systematic bias in the estimation; and by combining the upper and lower overprobabilities, we can assess the balance of the estimation.

10. A hydrological uncertainty analysis system based on probability density envelope, based on the hydrological uncertainty analysis method based on probability density envelope as described in any one of claims 1 to 9, characterized in that: include, The data acquisition and preprocessing module is used to acquire and preprocess the time series of hydrological elements from at least two hydrological observation stations. The kernel density estimation module is used to generate a uniform grid for each standardized sequence and perform kernel density estimation at the grid points to obtain the density estimate. The Bootstrap resampling module is used to generate Bootstrap samples and density matrices; The sample calculation module is used to calculate the standard error at grid points, and the standardized upper and lower biases for each Bootstrap sample; The optimization module is used to determine the optimal multiplier; The optimal multiplier calculation module is used to calculate the lower envelope, upper envelope, and excess probability of the probability density at grid points based on the optimal multiplier. The feature extraction module calculates feature parameters; The uncertainty analysis module is used to compare and analyze the uncertainty of hydrological elements at different stations based on the calculated characteristic parameters.

Images