Karst basin runoff segment number and attenuation coefficient optimization method based on dynamic planning

Abstract

According to the karst drainage basin runoff segmentation attenuation coefficient optimization method based on dynamic planning, the optimal segmentation number and segmentation division positions of the water recession curve are determined based on the dynamic planning algorithm, the influence of subjective judgment is avoided, the consistency of results can be kept, and the karst drainage basin runoff segmentation attenuation coefficient optimization method based on dynamic planning can be optimized according to karst drainage basins of different scales and characteristics. And rigid constraints such as the number of segments are flexibly applied, so that scientificity and accuracy of identification of different runoff components in the runoff attenuation process of the karst basin are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of watershed recession or river flow recession coefficient calculation, and particularly relates to a karst watershed runoff segmentation number and attenuation coefficient optimization method based on dynamic programming. BACKGROUND

[0002] Compared with non-karst basins, the unique geological structure (such as the developed fissure and pipe network) of karst basins shapes its complex water flow path and groundwater storage mechanism. This complexity makes it difficult for traditional models based on linear reservoir assumptions to accurately depict the recession process of karst basins. As a nonlinear runoff attenuation analysis tool, Maillet equation can more effectively represent this complexity, and is widely used in the attenuation analysis of river flow recession process and aquifer spring flow in karst areas. The core of Maillet equation is to transform the flow process by logarithm, and decompose the recession curve into multiple straight lines with different slopes. Each straight line corresponds to a specific stage in the flow attenuation process, and its slope is the attenuation coefficient of the stage. By fitting Maillet equation, the development degree of karst can be indirectly evaluated, and the types of karst aquifer (such as the pipe flow dominated fast drainage system and the fissure / matrix flow dominated slow drainage system) can be effectively distinguished. Further, by comparing the attenuation coefficients of different stages of the recession curve, the differences and spatial heterogeneity of the hydrological response of the "pipe-fissure" dual-porosity system in karst basins can be revealed, providing key basis for building and parameterizing the hydrological model of karst basins.

[0003] According to the difference in attenuation coefficients of different types of water-containing media such as surface, slope, and pores, fissures, solution gaps, and pipes in karst aquifer system, the recession curve of karst basin can be divided into multiple attenuation stages. However, how to accurately identify the number of different runoff attenuation stages in the recession curve and their segmentation positions, and then quantitatively calculate the recession coefficients, is still a key technical problem to be solved. In the recession curve analysis of karst basins, the traditional visual estimation method often relies on the subjective judgment of researchers, which not only affects the objectivity and accuracy of the results, but also increases the uncertainty of the analysis. For large-scale data sets, manual division not only consumes time, but also is difficult to ensure consistency. SUMMARY

[0004] To solve the technical problems proposed in the background art, it is necessary to provide a karst watershed runoff segmentation number and attenuation coefficient optimization method based on dynamic programming, to improve the scientificity and accuracy of identifying different runoff components in the runoff attenuation process of karst basins.

[0005] To achieve the above purpose, the technical scheme adopted by the present application is:

[0006] The karst basin runoff segmentation number and attenuation coefficient optimization method based on dynamic programming comprises the following steps.

[0007] S1, define the measured runoff time series data, the measured runoff time series data comprising a horizontal coordinate array t and a vertical coordinate array Q with a length of n, wherein the time index is used as the horizontal coordinate array t, and the measured recession curve flow observation value is used as the vertical coordinate array Q; define the maximum iteration number k_max of the outer loop controlling the segmentation number, k_max corresponding to the maximum medium type number, and set the iteration k as the recession curve segmentation number numSegments;

[0008] S2, establish a two-dimensional DP matrix DP[n+1, numSegments+1], and the element DP[i, j] of the two-dimensional DP matrix DP[n+1, numSegments+1] represents the minimum residual sum of squares of dividing the first i vertical coordinate data points of the measured runoff time series data into j segments;

[0009] S3, boundary condition initialization is performed on the two-dimensional DP matrix DP[n+1, numSegments+1];

[0010] S4, all segmentation numbers k are traversed through the outer loop, k∈[3, k_max], k=3 corresponding to the surface, soil and underground runoff respectively, and in each iteration process of the outer loop, the DP[i, j] is updated by using a multi-stage recursive optimization method;

[0011] S5, according to the obtained minimum residual sum of squares DP[i, j], the average fitting goodness R 2 of the segmentation results corresponding to all k∈[3, k_max] is calculated, the average fitting goodness R 2 with the maximum value is selected as the optimal segmentation point k value, and the segmentation number k and the segmentation position are obtained; after the segmentation number k and the segmentation position are obtained, the slope of each segment curve is obtained by linear fitting the scatter points between the segmentation positions, and the slope is the attenuation coefficient α; the segmentation number, the segmentation position and the attenuation coefficient α are taken as the segmentation result.

[0012] Further, in the boundary condition initialization process of step S3, the minimum residual sum of squares of DP[i, 0] is set to infinity, indicating an invalid state, and only DP[0, 0]=0 is a reasonable initial state; all DP[i, j] with j≥1 are initialized to positive infinity for subsequent recursive optimization.

[0013] Further, the DP[i, j] is updated by using the multi-stage recursive optimization method in step S4, comprising the following steps:

[0014] The DP[i, j] corresponding to the first segment curve of j = 1 is processed separately: directly calculate the first segment curve from the starting point to the current point i to perform linear fitting and calculate RSS to update DP[i, 1];

[0015] For j-segment curves of j≥2, the DP[i, j] is updated by using a multi-stage recursive optimization method, including the following steps:

[0016] Set constraints: ensure that each segment contains at least 2 data points, and the recursive range satisfies i≥2;

[0017] State transition: traverse all possible previous segmentation points k (2*(j-1)-1≤k≤i-3), calculate the RSS from k+1 to i-1, and then update DP[i, j] according to the following recursive formula as the state transition equation:

[0018] DP[i, j] = min k in(DP[i, j-1]+RSS[k+1, i-1]);

[0019] In the formula, RSS[a, b] represents the linear fitting residual sum of squares of data points a to b; m k The in(*) function represents traversing all k that satisfy the condition, calculating the sum of "the minimum residual of the first k segments (DP[k, j-1])" and "the last segment residual (RSS[k+1, i-1])" for each k, and finally taking the minimum value of these sums.

[0020] Further, in step S1, it is also necessary to verify the consistency of t and the array length of Q. If they are not equal, an exception is thrown to ensure data rationality.

[0021] Further, in step S2, a record matrix bestSplit[n+1, numSegments+1] is also established to store the optimal segmentation point k value and the position of the optimal segmentation point k obtained in step S5.

[0022] Further, in step S5, the output of the segmentation result includes the following steps:

[0023] The reverse order backtracking method is used, that is, starting from the endpoint currentIndex = n and currentSegments = numSegments, the segmentation points in the record matrix bestSplit matrix are backtracked segment by segment, and a segmentation list List1 is generated <segment>, List1 <segment>Each segment in the segment set comprises: a runoff attenuation coefficient a of each segment curve reflecting the runoff change trend in unit time, and a fitting goodness R corresponding to the segment curve for evaluating the explanation ability of the linear model to the actual data 2 value;

[0024] List1 <segment> <segment>, List2 <segment>Each segment in the list contains the following key parameters: the runoff decay coefficient a for each segment curve and the goodness of fit R corresponding to this segment curve 2 value;

[0025] Finally output the list of segments List2 <segment>as a result of the segmentation.

[0026] Further, in the process of separately processing the DP[i, j] corresponding to the first segment curve of j = 1, marking bestSplit[i, 1] = -1 indicates that the first segment has no pre-segmentation; in the process of backtracking the split points in the matrix bestSplit matrix and generating the segmentation list List1 by backtracking segment by segment in reverse order <segment>In the step of the above-mentioned method, the first section marker splitIndex=-1 needs to be processed separately, and the first section marker is reset to the starting point 0.

[0027] Further, in generating the section list List1 <segment>after the step of List1 <segment>Segment validity verification: ensure that the length of each segment currentIndex-splitIndex≥2, meet the minimum data requirements of linear fitting.

[0028] By adopting the technical scheme, the application has the following advantages compared with the traditional visual method:

[0029] (1) In terms of objectivity of segment division, the application maintains the consistency of results based on the global optimal solution calculated by dynamic programming; while the traditional visual method relies on the experience of analysts and is highly subjective.

[0030] (2) In terms of calculation efficiency, the application adopts computer program automatic processing, and the time complexity is controllable; while the traditional visual method needs manual analysis of each segment, which is time-consuming.

[0031] (3) In terms of division accuracy, the application uses the least square error to quantitatively evaluate the goodness of fit when fitting each segment; while the traditional visual method is easily affected by subjective judgment when identifying segment points, and lacks quantitative indicators.

[0032] (4) In terms of constraint processing ability, the application can flexibly impose hard constraints such as the number of segments according to different scales and characteristics of karst basins; while the traditional visual method is difficult to systematically handle complex constraints.

[0033] In summary, the karst basin runoff segment attenuation coefficient optimization method based on dynamic programming of the application determines the optimal number of segments and segment division positions of the recession curve based on the dynamic programming algorithm, is not affected by subjective judgment, can maintain the consistency of results, and can flexibly impose hard constraints such as the number of segments according to different scales and characteristics of karst basins, thereby improving the scientificity and accuracy of identification of different runoff components in the karst basin runoff attenuation process. BRIEF DESCRIPTION OF DRAWINGS

[0034] Figure 1 The flowchart of the karst basin runoff segment attenuation coefficient optimization method based on dynamic programming of the preferred embodiment of the application.

[0035] Figure 2 The fitting result of the multi-stage attenuation coefficient for the recession analysis of the 20190707 flood of the DS river basin.

[0036] Figure 3 The fitting result of the multi-stage attenuation coefficient for the recession analysis of the 20190714 flood of the DS river basin.

[0037] Figure 4 The fitting result of the multi-stage attenuation coefficient for the recession analysis of the 20200811 flood of the DS river basin.

[0038] Figure 5 The fitting result of the multi-stage attenuation coefficient for the DS river basin 20210618 flood recession analysis.

[0039] Figure 6 The fitting result of the multi-stage attenuation coefficient for the DS river basin 20220527 flood recession analysis. DETAILED DESCRIPTION

[0040] Please also refer to Figures 1 to 5 The method for optimizing the runoff segmentation number and attenuation coefficient of a karst basin based on dynamic programming provided by the preferred embodiment of the present application comprises the following steps:

[0041] S1, define the measured runoff time series data, which contains an array t of horizontal coordinates with a length of n and an array Q of vertical coordinates, wherein the time index is used as the array t of horizontal coordinates, and the measured recession curve flow observation value is used as the array Q of vertical coordinates; define the maximum iteration number k_max of the outer loop of the control segmentation number, k_max corresponds to the maximum medium type number, and the iteration k is used as the recession curve segmentation number numSegments.

[0042] In step S1, the consistency of the array lengths of t and Q also needs to be verified, and if they are not equal, an exception is thrown to ensure data rationality.

[0043] S2, establish a two-dimensional DP matrix DP[n+1,numSegments+1] and a record matrix bestSplit[n+1,numSegments+1], the element DP[i,j] of the two-dimensional DP matrix DP[n+1,numSegments+1] represents the minimum residual sum of squares of dividing the first i vertical coordinate data points of the measured runoff time series data into j segments.

[0044] In the embodiment, the record matrix bestSplit[n+1,numSegments+1] is also established in step S2, which is used to store the optimal split point k value and the position of the optimal split point k obtained subsequently.

[0045] S3, perform boundary condition initialization on the two-dimensional DP matrix DP[n+1,numSegments+1].

[0046] Specifically, in the boundary condition initialization process, for DP[i,0] (i.e., 0 segment division), the minimum residual sum of squares is set to infinity, indicating an invalid state, and only DP[0,0]=0 is a reasonable initial state; all DP[i,j] for j≥1 are initialized to positive infinity for subsequent optimization.

[0047] S4, traverse all segment numbers k, k ∈ [3, k_max], k = 3 respectively corresponding to 3 types of media, surface, soil and underground runoff, in each iteration process of the outer loop, update DP[i, j] by using the multi-stage recursive optimization method.

[0048] In the embodiment, updating DP[i, j] by using the multi-stage recursive optimization method includes the following steps:

[0049] S41, DP[i, j] corresponding to the first segment curve of j = 1 is processed separately: directly calculate the linear fitting of the first segment curve from the starting point to the current point i and calculate the RSS (residual sum of squares of linear fitting) to update DP[i, 1]; at the same time, in the process of processing DP[i, j] corresponding to the first segment curve of j = 1 separately, mark bestSplit[i, 1] = -1 to represent that the first segment has no previous segment;

[0050] S42, for j segment curve of j ≥ 2, update DP[i, j] by using the multi-stage recursive optimization method, including the following steps:

[0051] Set the constraint condition: ensure that each segment contains at least 2 data points, and the recursive range satisfies i ≥ 2 (that is, the previous j-1 segments each contain at least 2 points);

[0052] State transition: traverse all possible previous segment points k (2*(j-1)-1 ≤ k ≤ i-3), calculate the RSS from k+1 to i-1, and then update DP[i, j] according to the following recursive formula as the state transition equation:

[0053]

[0054] In the formula, RSS[a, b] represents the residual sum of squares of linear fitting of data points a to b; The function represents traversing all k satisfying the condition, calculating the sum of "the minimum residual of the first k segments (DP[k, j-1])" and "the residual of the last segment (RSS[k+1, i-1])" for each k, and finally taking the minimum value in these sums.

[0055] The calculation formula of RSS is as follows:

[0056]

[0057] In the formula, l is the sequence length, y o,k And y s,k Are the kth measured value and the simulation value in the sequence respectively.

[0058] S5, according to the minimum residual sum of squares DP[i, j] obtained, calculate the average fitting goodness R of the segmentation result corresponding to all k ∈ [3, k_max] 2 , the average R 2 The maximum is selected as the optimal split point k value, and the segment number k and the segment position are obtained. After obtaining the segment number k and the segment position, the slope of each segment curve obtained by linear fitting of the scattered points between the segment positions is the attenuation coefficient a. The segment number, the segment position and the attenuation coefficient a are taken as the output segment result.

[0059] R 2 The calculation formula is as follows:

[0060]

[0061] In the formula, is the mean value of the measured value sequence.

[0062] In this embodiment, in step S5, the output of the segment result includes the following steps:

[0063] S51, the reverse order backtracking method is adopted, that is, starting from the end point currentIndex (the index of the element currently being processed) = n and currentSegments (the current segment number) = numSegments (the total segment number), the split points in the record matrix bestSplit matrix are backtracked and the segment list List1 is generated <segment>, List1 <segment>Each segment in the segment list List1 comprises: a runoff attenuation coefficient a of each segment curve reflecting the runoff change trend in unit time, and a goodness of fit R corresponding to the segment curve for evaluating the explanatory ability of the linear model to the actual data 2 value.

[0064] In this embodiment, the split points in the record matrix bestSplit matrix are traced back in segments in a reverse order tracing back manner, and a segment list List1 is generated <segment>In the step of the above-mentioned method, the first segment marker splitIndex = -1 needs to be handled separately, and the first segment marker is reset to the starting point 0. In addition, in generating the segment list List1 <segment>after the step of List1 <segment>segment validity verification, ensuring that each segment length currentIndex- splitIndex≥ 2, thus satisfying the linear fitting minimum data requirement.

[0065] S52, List1 obtained by backtracking <segment> <segment>, List2 <segment>Each segment in the list contains the following key parameters: the runoff decay coefficient a for each segment curve and the goodness of fit R corresponding to the segment 2 value;

[0066] S53, finally output the segment list List2 <segment>As a segmentation result.

[0067] The following specific embodiment illustrates the method for optimizing the segmented attenuation coefficient of karst watershed runoff based on dynamic programming according to this invention.

[0068] The DS River basin, controlled by the DM hydrological station in Guangxi Zhuang Autonomous Region, was selected as a typical study area. Karst landforms account for 100% of the basin, and hydrological processes are influenced by both rapid conduit flow and slow matrix flow, resulting in a multi-stage decay characteristic in the receding water curve. Specifically, hourly flow observation data from the DM hydrological station from 2019 to 2022 were collected, and five representative flood events were selected based on the flow process curve characteristics for subsequent analysis and calculations. The maximum number of iterations in the outer loop controlling the number of segments was defined as k_max = 7.

[0069] After processing by the dynamic programming-based karst watershed runoff segmental attenuation coefficient optimization method of this embodiment, the average R of the segmented results for all k∈[3,7] is calculated. 2 Choose the average R 2 The largest value is taken as the optimal k value. The segments obtained from backtracking are then arranged in reverse order to generate a list of segments. <segment>each segment contains the attenuation coefficient (a) and the goodness of fit (R 2 ) of the linear regression model. Specifically, the analysis results of different segment numbers of 5 flood events in the DS watershed are shown in Table 1. From Table 1, it can be seen that when the segment number of 5 flood events is 4, the average R 2 is the largest, and the fitting results of the attenuation coefficient (a) of each segment are shown in Figure 1, in which a Figures 2 to 6 , a S , a I , and a GC are the attenuation coefficients of surface runoff, soil water runoff, fast pipe flow, and slow matrix flow, respectively. GM

[0070] Table 1 Average R 2 of different segment numbers of 5 flood events in the DS watershed

[0071] Flood episode numSegments = 3 numSegments = 4 numSegments = 5 numSegments = 6 numSegments = 7 20190707 0.93 0.99 0.97 0.96 0.96 20190714 0.94 0.95 0.91 0.92 0.93 20200811 0.94 0.98 0.96 0.96 0.95 20210618 0.97 0.98 0.97 0.95 0.94 20220527 0.97 0.98 0.95 0.96 0.94

[0072] The above description is a detailed description of the preferred embodiments of the present application, but the embodiments are not intended to limit the scope of the patent application of the present application. Any equivalent changes or modifications made under the technical spirit of the present application should be included in the scope of the patent application of the present application.​< / segment> < / segment> < / segment> < / segment> ​< / segment> < / segment> < / segment> < / segment> < / segment> < / segment> < / segment> < / segment> < / segment> < / segment> < / segment> < / segment> ​< / segment> < / segment> < / segment>

Claims

1. A method for optimizing the number of runoff segments and attenuation coefficients in karst watersheds based on dynamic programming, characterized in that, Includes the following steps: S1, define measured runoff time series data, which includes an abscissa array t and a ordinate array Q of length n, wherein the time index is used as the abscissa array t and the measured flow observation value of the drainage curve is used as the ordinate array Q; define the maximum number of iterations k_max of the outer loop that controls the number of segments, k_max corresponds to the maximum number of media types, and use the iteration k as the number of segments of the drainage curve; S2. Establish a two-dimensional DP matrix DP[n+1,numSegments+1]. The element DP[i,j] of the two-dimensional DP matrix DP[n+1,numSegments+1] represents the minimum residual sum of squares when the first i ordinate data points of the measured runoff time series data are divided into j segments. S3, initialize the boundary conditions of the two-dimensional DP matrix DP[n+1,numSegments+1]; S4, through the outer loop, traverse all segment counts k, k∈[3,k_max], k=3 correspond to surface runoff, soil runoff, and subsurface runoff respectively. In each iteration of the outer loop, a multi-stage recursive optimization method is used to update DP[i,j]. S5. Based on the obtained minimum residual sum of squares DP[i,j], calculate the average goodness of fit R of the piecewise results corresponding to all k∈[3,k_max]. 2 Choose the average goodness of fit R 2 The largest value is taken as the optimal segmentation point k, and the number of segments k and the segment positions are obtained. After obtaining the number of segments k and the segment positions, the slope of each segment curve obtained by linear fitting between the scattered points of each segment position is the attenuation coefficient α. The number of segments, the segment positions and the attenuation coefficient α are used as the output segmentation results.

2. The method for optimizing the number of runoff segments and attenuation coefficient in karst watersheds based on dynamic programming as described in claim 1, characterized in that, During the boundary condition initialization process in step S3, DP[i,0], i.e., segment 0, is divided, and the minimum residual sum of squares is set to infinity to represent an invalid state. Only DP[0,0] = 0 is a reasonable initial state. All DP[i,j] with j≥1 are initialized to positive infinity for subsequent recursive optimization.

3. The method for optimizing the number of runoff segments and attenuation coefficient in karst watersheds based on dynamic programming as described in claim 1, characterized in that, Step S4, which uses a multi-stage recursive optimization method to update DP[i,j], includes the following steps: The DP[i,j] corresponding to the first curve segment j=1 is processed separately: the first curve segment from the starting point to the current point i is directly calculated, linearly fitted, and RSS is calculated to update DP[i,1]. For the j-th curve segment j≥2, a multi-stage recursive optimization method is used to update DP[i,j], including the following steps: Set constraints: ensure that each segment contains at least 2 data points, and the recursive range satisfies i≥2; State transition: Traverse all possible preorder segment points k (2*(j-1)-1≤k≤i-3), calculate the RSS from k+1 to i-1, and update DP[i,j] using the following recursive formula as the state transition equation: DP[i,j]=m k in(DP[i,j-1]+RSS[k+1,i-1]); In the formula, RSS[a,b] represents the sum of squared residuals from the linear fit of data points a to b; m k The in(*) function iterates through all k that meet the conditions, calculates the sum of "the minimum residual of the first k segments (DP[k,j-1])" and "the residual of the last segment (RSS[k+1,i-1])" for each k, and finally takes the minimum value among these sums.

4. The method for optimizing the number of runoff segments and attenuation coefficient in karst watersheds based on dynamic programming as described in claim 1, characterized in that, In step S1, it is also necessary to verify the consistency of the array lengths of t and Q. If they are not equal, an exception is thrown to ensure the rationality of the data.

5. The method for optimizing the number of runoff segments and attenuation coefficient in karst watersheds based on dynamic programming as described in claim 3, characterized in that, In step S2, a record matrix bestSplit[n+1,numSegments+1] is also established to store the optimal split point k value and the position of the optimal split point k obtained in step S5.

6. The method for optimizing the number of runoff segments and attenuation coefficient in karst watersheds based on dynamic programming as described in claim 5, characterized in that, In step S5, the output of the segmentation results includes the following steps: The algorithm uses reverse backtracking, starting from the endpoints currentIndex = n and currentSegments = numSegments, and records the split points in the bestSplit matrix segment by segment, generating a segment list List1. <segment>,List1 <segment>Each segment includes: the runoff attenuation coefficient α, which reflects the trend of runoff change per unit time, and the goodness of fit R, which is used to evaluate the interpretability of the linear model for the actual data. 2 value;< / segment> < / segment> List1 obtained from backtracking <segment>Reverse the order of the segments in the array to generate a list of segments, List2. <segment>,List2 <segment>Each segment contains the following key parameters: the runoff attenuation coefficient α of each curve segment and the goodness of fit R corresponding to that curve segment. 2 value;< / segment> < / segment> < / segment> Finally, output the segment list List2. <segment> As a segmentation result.< / segment> 7. The method for optimizing the number of runoff segments and attenuation coefficient in karst watersheds based on dynamic programming as described in claim 6, characterized in that, During the separate processing of DP[i,j] corresponding to the first curve segment j=1, bestSplit[i,1]=-1 is marked as the first segment has no preceding segmentation; the segmentation points in the bestSplit matrix are recorded segment by segment using a reverse backtracking method, and a segmentation list List1 is generated. <segment> In the process, the first segment marker splitIndex = -1 needs to be processed separately, and the first segment marker is reset to the starting point 0.< / segment> 8. The method for optimizing the number of runoff segments and attenuation coefficient in karst watersheds based on dynamic programming as described in claim 6, characterized in that, In generating the segmented list List1 <segment>After the steps, also List1 <segment> Validation of segmentation effectiveness: Ensure that the length of each segment, currentIndex-splitIndex, is ≥2, to meet the minimum data requirement for linear fitting.< / segment> < / segment>

Images