Semi-automatic method for extraction, analysis and functional characterization of karst hydrological time series

Abstract

The application discloses a kind of karst hydrology time series extraction, analysis and semi-automatic method of functional characterization, which comprises: obtaining karst hydrology time series, and pretreatment;Extract flood sample time series from the karst hydrology time series after pretreatment;From flood sample time series, determine flood peak position, extract attenuation process according to flood peak position;Extract functional characterization from flood sample time series, including: Gaussian fitting is carried out to attenuation process, and attenuation coefficient is extracted in segments to extract effective functional characterization information.The application obtains the flood peak process of long time series of karst spring groundwater in the same region quickly by semi-automatic method, and then can carry out other frequency or time domain statistical analysis work on the basis, and obtain more dynamic characteristics of karst groundwater.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of karst groundwater dynamic characteristic extraction, and particularly relates to a semi-automatic method for karst hydrological time series extraction, analysis and functional characterization. BACKGROUND

[0002] Karst groundwater resources are high-quality freshwater resources available for human use and are the main source of water for living and production in eight provinces in Southwest China. In the process of karst groundwater development and utilization, the complexity and high connectivity of karst water-bearing media result in poor risk resistance and strong vulnerability of groundwater. Therefore, it has become increasingly necessary to optimize karst mountainous groundwater resource management, understand the recharge-runoff process of groundwater resources from the perspective of hydrogeology, and develop adequate resource protection strategies. Karst springs are the main discharge points of karst water systems. From the perspective of system theory, the groundwater of karst springs is the output of the whole karst water system and reflects the comprehensive regulation and storage function of the karst water system to rainfall input. Therefore, the hydrological time series of karst spring groundwater are often used to analyze the medium structure characteristics, regulation and storage capacity, and hysteresis of the karst water system. However, the hydrodynamic characteristics of karst groundwater are complex, and compared with the dynamics of surface rivers, karst groundwater has the characteristics of fast response time, long decay process, large variation range, and much noise. Under the control of different karst water systems, the flood hydrograph shape is specific and often has multiple superimposed peaks. Short-term and severe changes require that the monitoring time series of karst springs are generally at a frequency of 30 minutes, 60 minutes, or higher. Monitoring at a frequency of day scale or month scale often fails to reflect the change rule, but high-frequency monitoring also brings more challenges to the extraction of flood samples and the calculation of a large number of data samples. Meanwhile, the extraction of time series samples is inconsistent for the characteristics of different karst water systems. For example, the analysis of water-bearing medium characteristics often only targets the decay process of the underground flood peak, while the hysteresis analysis needs to analyze a single complete flood process, and the correlation analysis with rainfall needs the entire time series. The extraction of flood samples, the extraction of decay processes, and the segmentation of decay processes are the basis for all time series analysis. The original segmentation of rainfall-runoff time series cannot segment multiple flood peak processes and cannot automatically extract satisfactory flow decay processes, so as to automatically extract and characterize the function of the karst water system implied by the flood samples. Therefore, with the development of the Internet of Things technology and the exponential increase in collected data, it is urgent to develop an automatic data extraction method for such complex time series characteristics, so as to quickly obtain the characteristics and system function of different karst groundwater. SUMMARY

[0003] The technical problem solved by the present application is to provide an automatic karst hydrological data extraction method, so as to quickly analyze the characteristics and system functions of different karst underground water.

[0004] To solve the technical problem, the technical solution adopted by the present application is to provide a semi-automatic method for karst hydrological time series extraction, analysis and functional characterization, comprising the following steps:

[0005] S1: obtaining karst hydrological time series and performing pretreatment;

[0006] S2: extracting flood sample time series from the pretreated karst hydrological time series;

[0007] S3: determining the flood peak position from the flood sample time series, and extracting the decay process according to the flood peak position;

[0008] S4: extracting functional characterization from the flood sample time series, including: performing Gaussian fitting on the decay process, and extracting decay coefficients in segments to extract effective functional characterization information.

[0009] Further, step S1 specifically comprises:

[0010] Performing pressure influence elimination on the karst hydrological time series, and the specific calculation formula of the pressure influence elimination is as follows:

[0011] H Rt = H Jt -P t -ΔH

[0012] In the formula, H Rt , H Jt respectively represent the actual water level and the monitored water level of the karst water at time t; P t represents the air pressure of the water surface at time t; and ΔH represents the system error.

[0013] Performing sliding window smoothing processing on the karst hydrological time series after pressure influence elimination, and the selection of the window width depends on the frequency characteristics of the time series to be extracted;

[0014] During the sliding window smoothing processing, the detection and elimination of abnormal values are simultaneously performed, and the outliers are eliminated and interpolated according to the data before and after.

[0015] Further, step S2 specifically comprises:

[0016] For the pretreated karst hydrological time series, the base flow part of the flow process is obtained by using digital filtering method, and the specific operation is as follows:

[0017]

[0018] wherein b t is the base flow at time t, b t-1 is the base flow at time t-1, Q' t , Q' t-1 are the flow at time t, t-1 respectively, and β is a filter coefficient.

[0019] The filtering effect is controlled by adjusting the filtering times and the filter coefficient. After obtaining the base flow sequence, the following formula is used to separate the karst hydrological time series {Q'1, Q'2, …, Q' t}, thereby obtaining a new karst hydrological time sequence {Q1, Q2, …, Q t}.

[0020] Q t = Q' t -b t

[0021] The separated karst hydrological time sequence is extracted by using a three-point traversal method of flood process starting point, process determination point, and flood process ending point.

[0022] Further, the step of extracting the flood sample time sequence from the separated karst hydrological time sequence by using the three-point traversal method of flood process starting point, process determination point, and flood process ending point specifically includes the following steps.

[0023] The starting point of the flood sample is found, and the specific steps are as follows.

[0024] The first-order difference sequence {dif1, dif2, …, dif t} and the second-order difference sequence {diff1, diff2, …, diff t} of the karst hydrological time sequence are calculated, and the traversal starts at the starting point of the time sequence to preliminarily find the starting point Q(s) of the flood sample. The determination conditions of the starting point are as follows.

[0025] ① The flow at the next time is greater than the flow at the point Q(s), and the flow at the point Q(s) is greater than the threshold value of the basin, i.e., Q(s+1) > Q(s) and Q(s) > Q th .

[0026] ② The first-order difference at the point Q(s) is greater than the threshold value dif th , and the second-order difference is greater than 0, i.e., dif(s) > dif th and diff th > 0, indicating that the process of increasing flow starts.

[0027] ③ The flow at a certain time Δt after the starting point is greater than the flow at the starting point, and the first-order difference at the point Q(s+Δt) is also greater than 0, i.e., Q(s+Δt) > Q(s) and dif(s+Δt) > 0.

[0028] After the initial determination of the starting point of the flood process, verification is performed:

[0029] The first case is that the end point of the previous peak and the starting point of the next peak differ by more than a set value:

[0030] Perform reverse traversal to find the minimum value within the preset time step as the starting point of the real flood process;

[0031] The second case is that the two peaks are closely connected, and the end point of the previous peak is the starting point of the next peak:

[0032] The initial value of the process determination point Q(p) is set as the starting point position, and the first-order difference is greater than the threshold dif th and the second-order difference is greater than 0, and this point is also considered as the starting point Q(s) of the flood process.

[0033] Further, the step of extracting the flood sample time series from the separated karst hydrological time series using the flood process starting point, process determination point, and flood process end point three-point traversal method, specifically includes:

[0034] Finding the end point of the flood sample, specifically including the following steps:

[0035] The initial value of the process determination point Q(p) is set as the starting point position, and the first-order difference is greater than the threshold dif th and the second-order difference is greater than 0, and if the distance from the starting point position to this point position is greater than the flood sample width threshold W th , then this point is considered as the end point Q(e), and if the calculated flood sample width is lower than the width threshold, then the position of the process determination point is moved to this point, and from this point, continue to traverse backward to find the end point Q(e) that meets the requirements within the length of the time series.

[0036] The selection of the flood sample end point Q(e) considers two cases, one is that the flow naturally decays to the base flow or below, and the other is that the flow changes to stable before the next peak, and the latter considers the hydrological dynamic characteristics of multiple peaks, which specifically needs to meet the following two standards:

[0037] ① The flow or water level is lower than the set threshold, that is, the flow state meets the base flow, and it is considered that the flood process is over;

[0038] ② Before the next peak, the flow changes to stable, that is, the average value of the first-order difference and the second-order difference of the 2 positions before and after the point Q(e) is less than the corresponding threshold, that is:

[0039] mean(dif(e-2)+dif(e-1)+dif(e)+dif(e+1)+dif(e+2))<difth mean(diff(e-2)+diff(e-1)+diff(e)+diff(e+1)+diff(e+2))<diff th

[0040] wherein dif(e), diff(e) are the first-order and second-order difference values of the flood sample respectively, dif th , diff th are the set first-order and second-order difference thresholds respectively.

[0041] Further, the step S3 specifically comprises:

[0042] The first-order difference and neighborhood comparison method is adopted to find the peak point of the one-dimensional array, and the specific calculation steps are as follows: the first-order difference of the flood sample is calculated:

[0043] Δx f =x i+1 -x i

[0044] wherein x i represents the i-th data point of the flood sample;

[0045] The zero-crossing point in the first-order difference signal is detected, and if Δx i changes from positive to negative, it indicates that there is a peak value at this position; in addition, the minimum height, minimum distance, minimum width and other screening conditions of the peak value are set to further screen the peak value meeting the requirements, and the position sequence P idx of the peak value is obtained.

[0046] {Q(P idx (end)).....Q(j), P idx (end)<j<n}

[0047] wherein P idx (end) is the last value in the position sequence of the peak value, j is the coordinate of any point of the flood sample, and n is the length of the flood sample.

[0048] Further, the step S4 comprises:

[0049] The Gaussian fitting method is used to fit the attenuation process to obtain the fitting curve.

[0050] For the attenuation process with a length of m, the optimal segmentation points are P1 and P2, the enumeration method is adopted to search the attenuation process, and a self-defined statistical quantity MFE of multi-segment fitting error is constructed, which is defined as follows:

[0051]

[0052] wherein x i , x j , x k are the positions of the abscissa of three different attenuation stages, wherein 1 < x i < x j < x k < m, y0, y m are the ordinate values of the origin, P1, P2 and m point, respectively, y i , y j , y k are the ordinate values corresponding to x i , x j , x k in the three different attenuation periods.

[0053] The positions of the two points P1 and P2 satisfying the minimum MFE are the optimal segmentation points;

[0054] Then the water release rate of different water-containing media and the water release amount of different water-containing media in a single flood process are calculated according to the following formula, so as to study the differences of different karst springs:

[0055]

[0056]

[0057] wherein Q t represents the flow at time t after the attenuation starts; Q i represents the flow at the initial time of the i+1 attenuation period, Q0 is the flow at the initial time of the first attenuation period, Q1 is the flow at the initial time of the second attenuation period, and Q2 is the flow at the initial time of the third attenuation period; alpha i represents the flow attenuation coefficient of the i attenuation period; t i represents the end time of the i attenuation period, and t1, t2 and m are the end time of the first attenuation period, the end time of the second attenuation period and the length of the attenuation period, respectively.

[0058] The technical scheme provided by the present application has the following beneficial effects:

[0059] 1. The present application uses the method of reverse traversal and backward traversal within a certain range to prevent the turning point of the multi-peak flood process from being misjudged as the starting point or the ending point, thereby improving the extraction accuracy of the multi-peak flood sample, and the continuous flood sample can also be segmented.

[0060] 2. The attenuation process of the single flood sample can be automatically identified, and the time sequence similar to the attenuation process can be obtained by the method of Gaussian fitting distribution, thereby facilitating the subsequent feature extraction.

[0061] 3. The optimal segmentation point is determined by calculating the statistical quantity MFE through enumeration, and different decay processes are automatically divided, eliminating the influence of human subjective factors on the segmentation effect.

[0062] 4. By using a semi-automated method, the peak flow process of karst spring groundwater in the same area was quickly obtained over a long period of time. Based on this, other frequency or time-domain statistical analysis can be performed to obtain more dynamic characteristics of karst groundwater. Attached Figure Description

[0063] The present invention will be further described below with reference to the accompanying drawings and embodiments. In the accompanying drawings:

[0064] Figure 1 This is a flowchart illustrating the overall process of a semi-automatic method for extracting, analyzing, and functionally characterizing karst hydrological time series data according to the present invention.

[0065] Figure 2 This is a schematic diagram of the air pressure correction result of the present invention;

[0066] Figure 3 This is a schematic diagram of the results of flood sample extraction from a karst spring over a period of time according to the present invention;

[0067] Figure 4 This is a schematic diagram illustrating the flow attenuation process segmentation of the present invention. Detailed Implementation

[0068] To provide a clearer understanding of the technical features, objectives, and effects of the present invention, specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0069] like Figure 1 As shown in this embodiment, the semi-automatic method for karst hydrological time series extraction, analysis, and functional characterization includes the following steps:

[0070] S1: Obtain the time series of karst hydrology and perform preprocessing;

[0071] S2: Extract flood sample time series from preprocessed karst hydrological time series;

[0072] S3: Determine the location of the flood peak from the flood sample time series, and extract the attenuation process based on the location of the flood peak;

[0073] S4: Extract functional characteristics from the flood sample time series, including: performing Gaussian fitting on the attenuation process and extracting the attenuation coefficient in segments to extract effective functional characterization information.

[0074] Based on, but not limited to, this method, in step S1, the preprocessing of dynamic data of spring water in karst hydrological time series mainly includes three major steps: pressure influence removal, sliding window smoothing, and outlier detection and removal.

[0075] (1) Pressure influence elimination

[0076] In the actual work in the field, the pressure water level sensor is often used to monitor the water level change of karst groundwater. The monitored value includes the height of water surface to the monitoring position at a certain time and the size of the water surface pressure at that time. Therefore, the real water level change can be represented as formula (1), and the pressure correction result is as shown in Figure 2

[0077] H Rt = H Jt - P t - ΔH (1)

[0078] In the formula, H Rt , H Jt represent the actual water level and the monitored water level of karst water at t; P t represents the water surface pressure at t; and ΔH represents the system error, such as the error between the probe placement position and the actual flow section.

[0079] (2) Sliding window smoothing

[0080] The karst water time series after pressure influence elimination is subjected to sliding window smoothing. Sliding window smoothing is a low-pass filtering technique that can eliminate high-frequency noise and extract the long-term trend of the time series. This method obtains the smoothed time series by constructing a fixed-length sliding window and averaging the data in the window. The selection of the window width depends on the frequency characteristics of the time series to be extracted. Generally, a wider window can remove higher frequency noise. Considering that karst water has high frequency characteristics, it is recommended that the window be controlled within 4 hours in this embodiment.

[0081] (3) Abnormal value detection and elimination

[0082] Abnormal value detection and elimination are performed simultaneously in the sliding window smoothing. The purpose of abnormal value detection is to find the extreme points in the time series and judge whether they are real signals or noise. The outliers are removed, and interpolation is performed according to the data before and after them to prevent them from misleading the subsequent time series analysis.

[0083] Based on but not limited to the above method, in step S2, the preprocessed karst hydrological time series is first subjected to base flow component elimination, and then flood sample extraction. In view of the complex dynamic characteristics of karst water, a three-point traversal flood sample time series extraction method is first proposed.

[0084] (1) Base flow segmentation

[0085] ​For pretreated karst hydrological time series, the base flow part and flood part in the flow sequence are divided. The base flow part of the flow process is obtained by using digital filtering method, as shown in formula (2):

[0086]

[0087] In the formula, b t is the base flow at time t, b t-1 is the base flow at time t-1, Q' t , Q' t-1 are the flow at time t and t-1 respectively, and β is the filtering coefficient.

[0088] The filtering effect is controlled by adjusting the filtering times and filtering coefficient. After obtaining the base flow sequence, the karst hydrological time series {Q t} is separated by using formula (3), so as to obtain the new karst hydrological time series {Q t};

[0089] Q t =Q' t -b t (3)

[0090] For the separated karst hydrological time series, the flood sample time series is extracted by using the method of traversing the three points of flood process starting point (s), process determination point (p) and flood process ending point (e).

[0091] (1) The starting point of the flood sample is found, and the specific steps are as follows:

[0092] The first order difference sequence {dif1, dif2,..., dif t} and the second order difference sequence {diff1, diff2,..., diff t} of the karst hydrological time sequence are calculated, and the starting point Q(s) of the flood sample is found by traversing from the starting point of the time sequence. The determination conditions of the starting point are as follows:

[0093] 1) The flow at the next time is greater than the flow at the point Q(s), and the flow at the point Q(s) is greater than the threshold of the basin, that is: Q(s+1) > Q(s) and Q(s) > Q th ;

[0094] 2) The first order difference of the point Q(s) is greater than the threshold dif th , and the second order difference is greater than 0, that is: dif(s) > dif th and diff th > 0, which indicates that the flow begins to increase;

[0095] 3) the flow at a certain time Δt after the starting point is greater than the flow at the starting point, and the first-order difference of the point Q(s+Δt) is also greater than zero, i.e., Q(s+Δt)>Q(s) and dif(s+Δt)>0;

[0096] After the starting point of the flood process is preliminarily determined, verification is performed:

[0097] The first case is that the end point of the previous peak and the starting point of the next peak differ by more than a set value:

[0098] Reverse traversal is performed to find the minimum value within a certain time step (usually 10) as the starting point of the actual flood process;

[0099] The second case is that the two peaks are closely connected, and the end point of the previous peak is the starting point of the next peak:

[0100] The initial value of the process determination point Q(p) is set as the starting point position, and the first-order difference is greater than the threshold dif th and the second-order difference is greater than 0, the point is considered to be the starting point Q(s) of the flood process.

[0101] (2) Finding the end point of the flood sample, which includes the following steps:

[0102] The initial value of the process determination point Q(p) is set as the starting point position, and the first-order difference is greater than the threshold dif th and the second-order difference is greater than 0, if the distance from the starting point position to the point position is greater than the flood sample width threshold W th , then the point is considered to be the end point Q(e), if the calculated flood sample width is lower than the width threshold, then the position of the process determination point is moved to the point, and the search for the end point Q(e) that meets the requirements is continued from the point.

[0103] The selection of the flood sample end point Q(e) considers two cases, one is that the flow naturally decays to the base flow or below, and the other is that the flow changes to be stable before the next peak, the latter considers the hydrodynamic characteristics of multi-peak, which needs to meet the following two standards:

[0104] 1) the flow or water level is lower than the set threshold, i.e., the flow state meets the base flow, then the flood process is considered to be over;

[0105] 2) before the next peak, the flow changes to be stable, i.e., the average values of the first-order difference and the second-order difference of the point Q(e) before and after the point are less than the corresponding thresholds, i.e.,

[0106] mean(dif(e-2)+dif(e-1)+dif(e)+dif(e+1)+dif(e+2)) < dif th

[0107] mean(diff(e-2)+diff(e-1)+diff(e)+diff(e+1)+diff(e+2)) < diff th

[0108] wherein dif(e), diff(e) are first-order and second-order difference values of the flood sample respectively, dif th , diff th are set first-order and second-order difference thresholds respectively.

[0109] The above steps are repeated. During the traversal search process, the starting point Q(s), peak value and ending point Q(e) of each flood peak process are recorded, and the sequence of the length is extracted from the original data as the time sequence of the single flood process according to the starting point Q(s) and the ending point Q(e). The extraction effect of the method is shown in Figure 3 The automatic extraction result meets the requirements of subsequent refined data analysis.

[0110] For the extracted flood sample time sequence data, the flow decay process is also extracted, which is used for subsequent karst hydrology time sequence analysis.

[0111] The flow decay process is the minimum value of the flow after the maximum value of the flow in the flood sample with a length of n. For a single-peak flood sample, it is relatively simple to extract the decay process. For a multi-peak process, the flow of the first few peaks has not completely decayed, therefore, the application is more inclined to extract the decay process after the last peak, but the last peak process may not be the maximum value of the flood sample, so the application first identifies the inflection point, i.e. the peak point, in the flood sample.

[0112] Based on but not limited to the above method, step S3 judges the flood peak position from the flood sample time sequence, and extracts the decay process according to the flood peak position, and the specific implementation process is as follows:

[0113] The first-order difference and neighborhood comparison method is used to find the peak point of the one-dimensional array, and the specific calculation steps are as follows: the first-order difference of the flood sample is calculated:

[0114] Δx i =x i+1 -x i (4)

[0115] wherein x i represents the i-th data point of the flood sample;

[0116] Next, the zero-crossing points in the first-order differential signal are detected, if Δx i changes from positive to negative, it indicates that there is a peak value at the position; in addition, by setting screening conditions such as minimum height, minimum distance, minimum width of the peak value, the required peak value is further screened, and a position sequence P of the peak value is obtained idx , then the flood decay process can be expressed as:

[0117] {Q(P idx (end)).....Q(j),P idx (end)<j<n}

[0118] wherein P idx (end) is the last value in the position sequence of the peak value, j is the coordinate of any point of the flood sample, and n is the length of the flood sample.

[0119] Based on but not limited to the above method, in step S4, the functional characterization is extracted from the flood sample time sequence, including: performing Gaussian fitting on the decay process, and segmenting to extract the decay coefficient to extract effective functional characterization information;

[0120] The specific implementation process is as follows:

[0121] (1) Fitting of flow decay process

[0122] After extracting the flow decay process, the flow decay curve needs to be segmented and fitted to extract effective information. When automatically segmenting the flood decay process, first, in order to avoid the influence of local fluctuations on the overall decay process, the Gaussian fitting method is used to fit the decay process. As shown in Figure 4 , the fitting curve can completely conform to the overall trend of the decay process, and the data is more continuous and the curve is smoother, so that the change of the first-order difference in the lnQ-t curve is more continuous and more convenient to extract the position of the segmentation point.

[0123] (2) Automatic segmentation of flow decay process

[0124] For a decay process with a length of m, the optimal segmentation points are P1 and P2, the present application uses the enumeration method to search the decay process, and constructs a statistical quantity MFE (Multiple Fitting Error) of multi-segment fitting error, which is defined as follows:

[0125]

[0126] In the formula, x i , x j , and x k are the positions of the abscissa of the three different decay stages, wherein 1 < x i < xj <x k <m,y0、 y m The y-coordinates of the origin, P1, P2, and m are respectively. i y j y k The corresponding x in three different decay periods i x j x k The ordinate value.

[0127] The two points P1 and P2 that satisfy the minimum MFE are the optimal split points.

[0128] Then, according to formulas (5) and (6), the rate of water release from different water-bearing media and the amount of water released from different water-bearing media during a single flood can be calculated, so as to study the differences between different karst springs.

[0129]

[0130]

[0131] In the formula, Q t Q represents the flow rate at time t after the decay begins; i α represents the flow rate at the initial time of the (i+1)th decay period; i t represents the flow attenuation coefficient during the i-th attenuation period; i Let t1 represent the end time of the i-th decay period, and t2 represent the end time of the first decay period, the end time of the second decay period, and m represent the length of the decay period.

[0132] The key points of the embodiments of the present invention are as follows:

[0133] 1. Using the first and second differences of karst groundwater time series as judgment conditions, the starting and ending points of flood samples can be identified.

[0134] 2. The peak point of the flood is identified by first-order difference and neighborhood comparison, thereby realizing the extraction of the flood sample flow decay process.

[0135] 3. By using Gaussian fitting to fit the decay process, the impact of local fluctuations in the data on the segmentation effect is effectively eliminated.

[0136] 4. Use enumeration to search for the attenuation process and construct the multi-segment fitting error statistic MFE to find the optimal segmentation point, realize the automatic segmentation of the attenuation process, and obtain the attenuation coefficient of the segment.

[0137] The beneficial effects of implementing this invention are as follows:

[0138] 1. The application prevents the turning point of multi-peak flood process from being misjudged as the starting point or the end point by using the method of reverse traversal and a certain range of backward traversal, improves the extraction accuracy of multi-peak flood sample, and can also separate the continuous flood sample.

[0139] 2. The application can automatically identify the attenuation process of a single flood sample, and obtain a time sequence similar to the attenuation process by the method of Gaussian fitting distribution, so as to facilitate subsequent feature extraction.

[0140] 3. The application calculates the statistical quantity MFE by the enumeration method to determine the optimal segmentation point, automatically divides different attenuation processes, and eliminates the influence of human subjective factors on the segmentation effect.

[0141] 4. The application quickly obtains the long time sequence flood peak process of karst spring groundwater in the same region by the semi-automatic method, and can further perform other frequency or time domain statistical analysis work on the basis to obtain more dynamic characteristics of karst groundwater.

[0142] It should be noted that in this article, the terms "include", "contain" or any other variant thereof are intended to cover non-exclusive inclusion, so that the process, method, article or system including a series of elements not only includes those elements, but also includes other elements not explicitly listed, or includes elements inherent to such process, method, article or system. Without more limitations, the element defined by the statement "including a" does not exclude the presence of another identical element in the process, method, article or system including the element.

[0143] The above-mentioned application example serial numbers are only for description, and do not represent the advantages and disadvantages of the embodiments. In the unit claims of several devices, several of these devices can be embodied by the same hardware item. The use of the words first, second, and third does not represent any order, and these words can be interpreted as identifiers.

[0144] The above is only the preferred embodiment of the application, and does not limit the patent scope of the application. Any equivalent structure or equivalent process transformation using the content of the application specification and drawings, or direct or indirect application in other related technical fields, is also included in the patent protection scope of the application.

Claims

1. A semi-automatic method for karst hydrological time series extraction, analysis and functional characterization, characterized by, The method comprises the following steps: S1: obtaining a karst hydrological time sequence and preprocessing; S2: extracting a flood sample time sequence from the preprocessed karst hydrological time sequence; S3: determining a flood peak position from the flood sample time sequence and extracting a decay process according to the flood peak position; S4: extracting a functional representation from the flood sample time sequence, comprising: performing Gaussian fitting on the decay process, and extracting an effective functional representation information by segmenting and extracting decay coefficients; Step S2 specifically comprises: For the preprocessed karst hydrological time sequence, a digital filtering method is used to obtain the base flow part of the flow process, and the specific operation is as shown in the formula: wherein b t is the base flow at time t , b t-1 is the base flow at time t- 1, , is the flow at time t , t- 1, β is the filter coefficient; The filtering effect is controlled by adjusting the filtering times and filtering coefficients. After obtaining the base stream sequence, the karst hydrological time is separated by using the following formula , thereby obtaining a new karst hydrological time sequence . The separated karst hydrological time sequence is extracted by using the three-point traversal method of the flood process starting point, process determination point and flood process ending point; The step of extracting the flood sample time sequence from the separated karst hydrological time sequence by using the three-point traversal method of the flood process starting point, process determination point and flood process ending point specifically comprises: The starting point of the flood sample is found, and the specific steps are as follows: First-order difference sequence of karst hydrology time series and second-order difference sequence , starting from the beginning of the time series, initially searching for the starting point of the flood sample Q(s) The determination conditions of the starting point are as follows: ① the flow of the next time is greater than the flow of the point Q(s), and the flow of the point Q(s) is greater than the threshold of the basin, that is: and ; (ii) the first-order difference of the point Q(s) is greater than a threshold dith and the second-order difference is greater than 0, i.e., represents the process in which the flow rate starts to increase; ③ the flow rate at a certain time Δt after the starting point is greater than the flow rate at the starting point, and the first-order difference of the point Q(s+Δt) is also greater than zero, i.e. ; After the starting point of the flood process is preliminarily determined, verification is performed: The first case is that the difference between the ending point of the previous peak and the starting point of the next peak is greater than a set value: Reverse traversal is performed to find the minimum value within a preset time step as the starting point of the actual flood process; The second case is that the two peaks are closely connected, and the ending point of the previous peak is the starting point of the next peak: Process decision point Q ( p ) Initial value is set as the starting position, and then traverse to find a point where the first-order difference is greater than the threshold value dif th and the second-order difference is greater than 0, which is also considered as the starting point of the flood process Q ( s ); The step of extracting the flood sample time sequence from the separated karst hydrological time sequence by using the three-point traversal method of the flood process starting point, process determination point and flood process ending point specifically comprises: The ending point of the flood sample is found, and the specific steps are as follows: Setting the initial value as the start position, traversing backward to find a point where the first-order difference is greater than a threshold value and the second-order difference is greater than 0 Q p ①If the flow or water level is lower than the set threshold value, that is, the flow state meets the base flow, it is considered that the flood process is ended; th If the distance from the start position to the point is greater than a flood sample width threshold value W th , the point is considered as the end position Q e If the calculated flood sample width is lower than the width threshold value, the position of the process decision point is moved to the point, and the traversal is continued from the point to find a required end position in the time series length Q e ;​​​​ The end of the flood sample Q e The selection of the end of the flood sample takes into account two cases: one is that the flow naturally decays to the base flow or below, and the other is that the flow variation tends to be stable before the next peak, which takes into account the hydrological dynamic characteristics of multi-peak, and specifically meets the following two standards:​ ②Before the next peak, the flow changes tend to be stable, that is, the average values of the first-order difference and the second-order difference of the two positions before and after the point Q(e) are less than the corresponding threshold values, that is: dif wherein, diff ( e ), dif ( e ) are first-order and second-order difference values of the flood sample, respectively, diff th、 Step S3 specifically comprises: th are set first-order and second-order difference threshold values, respectively; The peak point of the one-dimensional array is found by using the first-order difference and neighborhood comparison method, and the specific calculation steps are as follows: the first-order difference of the flood sample is calculated: end wherein, represents the i-th data point of the flood sample; and i represents the i-th data point of the flood sample; and detecting zero-crossing points in the first-order differential signal, if from positive to negative, indicating that there is a peak value at this position; in addition, by setting screening conditions such as minimum height, minimum distance, and minimum width of the peak value, further screening of the required peak value is performed to obtain a position sequence of the peak value P idx , then the flood attenuation process can be represented as: wherein P idx Step S1 specifically comprises: ) is the last value in the sequence of positions of the peak, j is the coordinate of an arbitrary point of the flood sample, n is the length of the flood sample.​ 2. The semi-automatic method for karst hydrograph extraction, analysis and functional characterization according to claim 1, characterized in that, The karst hydrological time sequence is subjected to pressure influence elimination, and the specific calculation formula of the pressure influence elimination is as follows: The karst hydrological time sequence subjected to the pressure influence elimination is subjected to sliding window smoothing processing, and the window width is selected according to the frequency characteristics of the time sequence to be extracted; In the formula: H Rt , H Jt respectively represent t actual water level of karst water and monitored water level at the moment; P t represent t atmospheric pressure of water surface at the moment; Δ H represents system error; During the sliding window smoothing processing, the detection and elimination of abnormal values are simultaneously performed, and the outliers are eliminated and interpolated according to the data before and after the outliers. Step S4 comprises:

3. The semi-automatic method for karst hydrograph extraction, analysis and functional characterization according to claim 1, characterized in that, The decay process is fitted by using the Gaussian fitting method to obtain a fitting curve; Then, the rates of water release of different water-containing media and the water release amounts of different water-containing media in a single flood process are calculated according to the following formula, so as to study the differences of different karst springs: For the decay process with length m , the optimal segmentation points are P 1, P 2, respectively. The enumeration method is used to search the decay process, and a self-defined statistical quantity MFE of multi-segment fitting error is constructed, which is defined as follows: wherein x i , x j , x k are the positions of the abscissa for the three different decay phases, respectively, wherein 1 x i x j x k m , y 0, , , y m are the ordinate values of the points 0, P 1, P 2 and m , respectively, y i , y j , y k are the ordinate values of the points corresponding to x i , x j , x k in the three different decay phases, respectively.​​​ The two points that satisfy the minimum MFE P 1 and P 2, i.e. the position of the optimal split point; ​ wherein denotes the time of the beginning of the decay t flow rate at the time; denotes the first i+ 1 flow rate at the beginning of the first decay period, Q0 flow rate at the beginning of the second decay period, Q1 flow rate at the beginning of the third decay period, Q2 denotes the first i decay period flow rate decay coefficient; denotes the first i end of the decay period, respectively the end of the first decay period t 1, the end of the second decay period t 2, the length of the decay period m .

Images