A method for hydrological time series motif mining based on numerical feature clustering

Abstract

The present invention discloses a hydrological time series motif mining method based on numerical feature clustering, including a time series clustering process and a time series motif mining process; in the time series motif mining process, firstly, the time series in each class are fully connected and the global distance matrix is ​​calculated based on the class; secondly, a variable-size two-dimensional window is used to slide on the global matrix according to certain rules and constraints; then the shortest path is calculated for the local matrix in each window; finally, a candidate motif list is constructed according to the calculation results, and a suitable motif is selected for applications in different fields. Starting from the two perspectives of the numerical features of time series concerned in the hydrological field and the difficulty in determining the motif length in practical applications, the present invention introduces a clustering method based on the numerical features of time series and a DTW algorithm for solving the distance calculation of non-uniform time series into the existing motif mining method to achieve efficient and accurate mining effects.

Description

Technical Field

[0001] The present invention belongs to the technical field of hydrological data processing, and in particular relates to a hydrological time series motif mining method based on numerical feature clustering. Background Art

[0002] With the rapid development of Internet of Things technology, various sensor devices have recorded massive amounts of data that change over time, namely time series. An important research point in the field of time series analysis is motif mining, which is of great significance for discovering features in time series. Motifs refer to repeated sub-segments or frequently occurring trends in time series. Nowadays, motif mining has been widely used in many fields such as meteorology, earthquakes, and insect behavior analysis.

[0003] Time series motif mining methods can be divided into approximate motif mining and precise motif mining, and time series data can be one-dimensional data or multi-dimensional data. There are many classic methods in the field of time series motif mining. Patel et al. proposed an efficient fixed-length motif discovery algorithm. They first used a brute force algorithm to locate motif pairs. Shasha and Wang et al. optimized the brute force algorithm by using the approximate distance graph (ADM) and proposed the matrix approximation method (EMMA). Nevertheless, the motif discovery process in actual scenarios does not allow such high time calculation. Therefore, people have made many efforts to improve efficiency. Chiu et al. used symbolic aggregation approximation (SAX) to reduce the dimension of time series data and proposed the random projection (RP) algorithm. This method reduces the mining accuracy while compressing the data. Mueen's team pioneered the MatrixProfile algorithm, which reflects the distance between subsequences in the distance matrix, and completed the design of some pruning methods, such as excluding repeated subsequence pairs (trivial match). Since then, MP-based methods have been continuously improved, and the accuracy and efficiency have gradually improved. The most typical ones are the STOMP and STAMP methods proposed by Zhu et al.

[0004] However, in the field of hydrological applications, it is difficult for users to determine the length of the motif at the beginning. On the other hand, paired motifs are not necessarily two subsequences of strictly equal length. Subsequences of unequal lengths can also reflect similar hydrological characteristics. Therefore, variable-length motif discovery is particularly important. Linardi et al. proposed a variable-length motif discovery (VALMOD) method based on MP. This method searches for all possible motifs in a given range. Another feasible method is based on the dynamic time warping (DTW) algorithm, which can calculate the distance between sequences of unequal lengths. Alaee et al. proposed the SWAMP algorithm based on this. However, the above methods still do not solve the mining of dynamic-length motifs.

[0005] In addition, due to the domain specificity of hydrology, hydrological time series have strict interpretation significance, that is, they can be divided into units according to domain knowledge. A hydrological time series within a period of time can reflect an event or even describe a single scene. Most of the existing motif mining methods are based on long time series, which is difficult to handle time series with multiple features and a wide time span in the hydrological field. A feasible solution is to cluster the relevant numerical features of the hydrological time series in advance and mine motifs with hydrological interpretable features in the cluster set. Zhou et al. proposed the IF2CNN framework, which integrates the iterative filtering (IF) method and convolutional neural networks (CNNs) for automatic feature learning of time series. However, the CNN feature extraction process mainly serves specific prediction tasks and it is difficult to extract other important features besides temporal features. Tiano et al. proposed a feature-based semi-supervised clustering framework (FeatTS), which relies on time series labels to adaptively adjust the key features suitable for time series and realizes the clustering of time series based on existing clustering methods. However, this supervised method requires partial time series labels, which are difficult to meet in some applications. In addition, the community detection and correlation calculation steps are separated, and the excessive parameters involved can be considered for further integration. Summary of the invention

[0006] Purpose of the invention: To overcome the deficiencies of the above prior art, the present invention provides a hydrological time series motif mining method based on numerical feature clustering to achieve real-time, accurate and efficient motif mining.

[0007] Technical solution: The present invention provides a hydrological time series motif mining method based on numerical feature clustering, comprising the following steps:

[0008] (1) Extract the numerical features of multiple hydrological time series, select the features that are highly correlated with the time series labels according to the time series labels, and on this basis select the minimum number of features greater than the coverage threshold;

[0009] (2) Construct an edge-weighted graph with the minimum number of features as the basic unit, determine the retained distance threshold based on the weighted edge screening ratio, and perform initial community detection on the weighted graph of each feature;

[0010] (3) Calculate feature weights based on the number of communities in each weighted graph and the number of predefined clusters, and construct a weighted co-occurrence matrix; calculate the similarity between row vectors in the co-occurrence matrix to achieve time series clustering set;

[0011] (4) Taking each clustering set in step (3) as the basic unit, fully connect the time series in the set to form a complete time series, and compare it with itself to form a global distance matrix;

[0012] (5) Use a sliding window of variable size to slide on the global distance matrix, and calculate the shortest path distance of the local distance matrix within the window during each sliding process;

[0013] (6) A set of subsequence pair triplets is constructed based on the shortest path distance, and a list of candidate motifs is formed by screening, and a predefined number of motifs are selected according to the actual application.

[0014] Furthermore, the step (1) comprises the following steps:

[0015] (11) Form relevant numerical features based on the duration, amplitude and trend of the time series of interest in the hydrological field;

[0016] (12) Calculate the relevance of features based on the labels of the time series, i.e., the p-value, and select the top features after sorting them by p-value. f features; to reduce feature redundancy, select the minimum number of features that meets the coverage δ.

[0017] Furthermore, the step (2) comprises the following steps:

[0018] (21) Construct edge weighted graph G w :For a feature f i And the time series set TSS = {TS1, TS2, ..., TS m}, each time series TS p For edge weighted graph A node v of the vertex set V p ; Each edge e in the weighted edge set E p,q Corresponding to a weight Keep the edges whose distance is less than the threshold ζ;

[0019] (22) A community detection algorithm is applied to assign several nodes to different communities in each weighted graph. The community detection algorithm uses the connection structure of the graph to unify different representations of time series clustering.

[0020] Furthermore, the step (3) comprises the following steps:

[0021] (31) The weight of each feature is determined by calculating the ratio between the number of user-defined clusters and the number of communities. The specific calculation method is:

[0022]

[0023] Where C represents the number of predefined clusters, O i Representing a weighted graph the number of communities in

[0024] (32) Count the number of times each pair of time series appears in the same community and fill in the co-occurrence weighted matrix The specific method is:

[0025]

[0026] in, is the co-occurrence weighted matrix, Represents the time series TS i and time series TS j The weighted co-occurrence ratio between ; the calculation formula is:

[0027]

[0028] Among them, w k Represents the weight corresponding to any feature, w k' Indicates TS i and TS j The weight corresponding to any feature that appears in a community, m and m' represent the total number of features and TS respectively. i and TS j The total number of times they co-occur in a community;

[0029] (33) Yes The distance between rows is calculated, and the clustering algorithm is used to form the final clustering set.

[0030] Furthermore, the distance matrix in step (4) is:

[0031] In a single clustering set cs i In the above example, all time series in the set are arranged in chronological order to form a complete time series TS. glob ={ts1,ts2,…,ts m}; Complete time series TS glob Each element in TS glob Calculate the Euclidean distance for each element in and initialize the global distance matrix Mat dis_glob :

[0032]

[0033] in, It represents the distance between the sequence values ​​of the time series at time i and time j, and n represents the length of the time series.

[0034] Furthermore, the step (5) comprises the following steps:

[0035] (51) Use windows of different area sizes and different side length ratios to slide on the global distance matrix with a certain step size st;

[0036] (52) During a sliding process, the local distance matrix Mat within the sliding window is calculated dis_local The regular path distance, that is, the shortest path distance d sp :

[0037] d sp =Min{DisWarping(Mat dis_local )}

[0038] (53) gradually reducing the sliding window size according to the side length ratio and area setting, and repeating steps (51) and (52);

[0039] (54) For each window in step (52) and step (53), the local distance matrix Mat is constructed dis_local Calculate the shortest path distance d sp ; The shortest path distance uses the shortest path distance calculation optimization method, that is, the cumulative distance of the accumulated matrix elements in the reused previous large window and the local shortest path of the matrix in the reused previous neighbor window; the shortest path distance d sp The length needs to be normalized according to the subsequence.

[0040] Furthermore, the step (6) comprises the following steps:

[0041] (61) Construct subsequence pair triple SPT = <(idx1,len1),(idx2,len2),d sp >; (idx, len) represent the starting index and length of the subsequence respectively, d sp is the DTW distance between subsequence pairs;

[0042] (62) The subsequence pairs are screened according to the predefined candidate motif similarity threshold τ to form a candidate motif triple CMT, i.e., the normalized shortest path distance sim(d sp )≥τ; where sim is the similarity solution function, and the similarity is inversely proportional to the DTW distance;

[0043] (63) Select subsequence elements in the candidate motif triple CMT and add them to the candidate motif list CML;

[0044] (64) Each element of the candidate motif list CML (cmidx i ,cmlen i ) Associate a one-way linked list simlinklist i , each element in the linked list is the candidate motif triple CMT and (cmidx i ,cmlen i ) similar subsequences, that is, satisfying in, is any element of a singly linked list;

[0045] (65) Sort by the length of the one-way linked list, and according to the number of motifs K defined by the user, select the first K elements of the sorted candidate motif list CML, which is the K-motif.

[0046] Furthermore, the shortest path distance calculation optimization method described in step (54) comprises the following steps:

[0047] (541) According to the DTW algorithm shortest path solution process, the calculation domain Mat is calculated. dis_local The cumulative distance of all elements in [,:r+1];

[0048] (542) Multiplexing Field Mat dis_local The cumulative distance of all elements in [,N-st-1:N-st] is solved, and the previous order distance matrix is ​​reused from the calculation domain Mat dis_local The difference in the cumulative distance from [,:r+1] to the element, that is, in and They represent the previous order cumulative distance matrix and the current cumulative distance matrix respectively; (i, j) and (p, q) are located in the calculation domain and reuse domain of the previous order cumulative matrix respectively;

[0049] (543) According to the DTW shortest path solving process and the cumulative distance of all elements in the reuse domain known in step (542), the defined computational domain Mat is calculated. dis_local Cumulative distance of all elements in [,N-st:], Mat dis_local The cumulative distance of the elements in the Mth row and the Nth column is the shortest path distance d of the matrix. sp .

[0050] Beneficial effects: Compared with the prior art, the present invention has the following beneficial effects: the present invention utilizes a variable-size window to slide on the global distance matrix generated by DTW, realizes the similarity calculation of subsequence pairs of a one-dimensional time series, and maps them to the shortest path solution of the local matrix inside the two-dimensional sliding window; the local distance matrix calculation process adopts an optimization algorithm to simplify the solution of the shortest path, thereby realizing real-time, accurate and efficient model mining; at the same time, the method has practical application value in the field of hydrology. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 is a flow chart of the present invention;

[0052] Figure 2 It is a schematic diagram of the community detection results of various characteristics of the time series of the present invention;

[0053] Figure 3A schematic diagram of constructing a global distance matrix for the complete time series of the present invention;

[0054] Figure 4 It is a schematic diagram of dividing the reuse domain and calculation domain of the local distance matrix of the present invention;

[0055] Figure 5 Schematic diagram of the candidate motif list of the present invention. DETAILED DESCRIPTION

[0056] The present invention will be further described in detail below in conjunction with the accompanying drawings.

[0057] The present invention proposes a hydrological time series motif mining method based on numerical feature clustering, such as Figure 1 As shown, the specific steps include:

[0058] Step 1: Extract the numerical features of multiple hydrological time series, select the features that are highly correlated with the time series labels according to the time series labels, and on this basis select the minimum number of features greater than the coverage threshold.

[0059] In the time series set TSS = {TS1, TS2, ..., TS m}Use the tsfresh library to extract features, each time series TS i Corresponding to a eigenvector F i ={f1,f2,…,f n}; The feature vectors of the time series and the corresponding classification labels are combined into a two-dimensional matrix of m×(n+1), where m represents the total number of time series and n represents the total number of features. The Benjamini-Yekutieli procedure is used to quantify the importance of each feature according to the category label, i.e., the p-value. In the relevance sorting, the top f Features, top f are predefined parameters. The top f The feature screening is performed to ensure that the sum of the variances of the features after screening can cover most of the information of the remaining features. The coverage rate δ needs to be defined in advance. f After filtering, the number of features is top δ .

[0060] Step 2: Construct an edge-weighted graph with the minimum number of features as the basic unit, determine the retained distance threshold based on the weighted edge screening ratio, and perform initial community detection on the weighted graph of each feature.

[0061] Construct edge weighted graph G w , for a feature f i And the time series set TSS = {TS1, TS2, ..., TS m}, each time series TS p For edge weighted graph A node v of the vertex set V p . Each edge e in the weighted edge set E p,q Corresponding to a weight Edges with a distance less than the threshold ζ are retained. Community detection is performed on each weighted graph using the greedy modularity algorithm of the NetworkX library to search for combinations of nodes that are densely connected into communities, such as Figure 2 shown.

[0062] Step 3: Calculate feature weights based on the number of communities in each weighted graph and the predefined number of clusters, and construct a weighted co-occurrence matrix; calculate the similarity between row vectors in the co-occurrence matrix to achieve time series clustering set.

[0063] According to the community detection algorithm from the graph The number of communities derived from i And the predefined number of clusters C, for each feature f i Assign an approximate weight w i :

[0064]

[0065] Count the number of times each pair of time series appears in the same community and substitute it into the weight w i , construct a weighted co-occurrence matrix The specific method is:

[0066]

[0067] Among them, the element in row i and column j The calculation method is:

[0068]

[0069] Among them, w k Represents the weight corresponding to any feature, w k' Indicates TS i and TS j The weight corresponding to any feature that appears in a community, m and m' represent the total number of features and TS respectively. i and TS j The total number of times they co-occur in a community.

[0070] The weighted co-occurrence matrix The Euclidean distance is calculated between each row vector of . The K-Medoid algorithm is used to cluster similar time series and generate the final clustering set CS.

[0071] Step 4: Taking each cluster set in step 3 as the basic unit, fully connect the time series in the set to form a complete time series, and compare it with itself to form a global distance matrix, such as Figure 3 shown.

[0072] Each time series set c in the clustering set CS i ={TS1,TS2,…,TS n All time series in} are fully connected to form a complete time series TS glob ={ts1,ts2,…,ts m}.

[0073] The complete time series TS glob Each element in TS glob Calculate the Euclidean distance for each element in and construct the global distance matrix Mat dis_glob . Construct the global distance matrix Mat dis_glob The specific method is:

[0074]

[0075] in, It represents the Euclidean distance between the sequence values ​​of the time series at time i and time j, that is, d(i,j)=|ts i -ts j |; n represents the length of the time series.

[0076] Step 5: Use a sliding window of variable size to slide on the global distance matrix, and calculate the shortest path distance of the local distance matrix within the window during each sliding process.

[0077] (5.1) Use windows of different sizes and different side length ratios to slide on the global distance matrix. Take the upper left corner of the window as the reference point, the area is from large to small, and the direction is from left to right and from top to bottom. To improve time efficiency, slide according to a certain step length st;

[0078] (5.2) During a sliding process, the local distance matrix Mat within the sliding window is calculated dis_local The regular path distance, that is, the shortest path distance d sp The calculation method uses the shortest path solution method in the DTW algorithm:

[0079] d sp =Min{DisWarping(Mat dis_local )}

[0080] (5.3) Repeat the above steps (5.1) and (5.2). In the process of sliding between the left and right adjacent matrices, the distance matrix calculation optimization strategy is adopted. Figure 4 As shown, the local distance matrix It is divided into a reuse domain and a computational domain. The reuse domain is the (N-st)th column elements. The computational domain is the first (r+1)th column elements and the last st column elements defined by the fastDTW algorithm. The reuse and computation processes are performed based on the above demarcated areas.

[0081] The same applies to sliding between adjacent matrices;

[0082] (5.4) Repeat the above steps (5.1), (5.2) and (5.3). In the process of window size decreasing from large to small, the shortest path distance of the smaller local distance matrix can reuse the distance of the larger local distance matrix, that is, Among them, cum large (i,j) is the element of the i-th row and j-th column of the cumulative distance matrix when calculating the shortest path of the larger distance matrix. i and j are the number of rows and columns of the smaller local distance matrix.

[0083] The shortest path reuse and calculation process includes the following steps:

[0084] 1) Calculate the defined calculation Mat according to the shortest path solution process of the DTW algorithm dis_local The cumulative distance of all elements in [,:r+1];

[0085] 2) Multiplexing domain Mat dis_local The cumulative distance of all elements in [,N-st-1:N-st] is solved, and the previous order distance matrix is ​​reused from the calculation domain Mat dis_local The difference in the cumulative distance from [,:r+1] to the element, that is, in and Represent the previous order cumulative distance matrix and the current cumulative distance matrix respectively. (i,j) and (p,q) are located in the calculation domain and reuse domain of the previous order cumulative matrix respectively;

[0086] 3) According to the DTW shortest path solving process and the cumulative distance of all elements in the reuse domain known in step (532), the defined computational domain Mat is calculated dis_local The cumulative distance of all elements in [,N-st:]. Mat dis_local The cumulative distance of the elements in the Mth row and the Nth column is the shortest path distance d of the matrix. sp .

[0087] Step 6: Construct a set of subsequence pair triplets based on the shortest path distance, and filter them to form a list of candidate motifs, and select a predefined number of motifs based on actual applications.

[0088] (6.1) Construct subsequence pair triple SPT = <(idx1,len1),(idx2,len2),d sp >. Among them, (idx,len) represent the starting index and length of the subsequence respectively, d sp is the DTW distance between subsequence pairs;

[0089] (6.2) The subsequence pairs are screened according to the predefined candidate motif similarity threshold τ to form a candidate motif triple CMT, i.e., the normalized shortest path distance sim(d sp )≥τ. Where sim is the similarity solution function, and the similarity is inversely proportional to the DTW distance;

[0090] (6.3) Select subsequence elements in the candidate motif triple CMT and add them to the candidate motif list CML, such as Figure 5 As shown. Due to the binary symmetry of the candidate motif triple CMT, (idx1, len1) and (idx2, len2) need to be added to the CML;

[0091] (6.4) Each element (cmidx i ,cmlen i ) Associate a one-way linked list simlinklist i , each element in the linked list is the candidate motif triple CMT and (cmidx i ,cmlen i ) similar subsequences, that is, satisfying in, is any element of a singly linked list;

[0092] (6.5) Sort by the length of the one-way linked list. According to the number of motifs K defined by the user, select the first K elements of the sorted candidate motif list CML, which is the K-motif.

[0093] Taking rainfall, the dominant factor in flood forecasting applications in the field of hydrology, as an example, the details are as follows:

[0094] A total of 31 flood data sets from 1998 to 2010 in the Changhua Basin were selected, mainly including rainfall data and outlet flow values ​​at seven rain gauge monitoring points. The seven rain gauge monitoring points are: Daoshiwu, Taohua Village, Longmen Temple, Shuangshi, Lingxia, Yulingguan and Changhua. The monitoring frequency is one hour. The average rainfall in the Changhua Basin is calculated based on the spatial distribution ratio of the seven stations, and 31 average rainfall time series data sets TSS = {TS1, TS2, …, TS 31}. Among them, the surface average rainfall time series TS i Each data point The calculation formula is:

[0095]

[0096] in, Indicates the hourly rainfall measured at the monitoring station. rt i It indicates the ratio of the sub-basin area where the monitoring station is located to the total basin area.

[0097] Use tsfresh library to extract rainfall time series data set TSS = {TS1, TS2, ..., TS 31} feature set The p-value of the feature set F is calculated using the Benjamini-Yekutieli method. According to empirical research, the first 20 features sorted by relevance are sufficient to obtain high-quality clustering. Therefore, the first 20 features are selected. The feature set F is screened using the PFA algorithm. Based on various threshold experiments, the coverage rate is selected to be 90%. Therefore, the minimum number of features is selected, and the sum of their variances covers 90% of the information generated by the remaining features. The screened feature set The three features are quantile, trend_stderr, and trend_rvalue.

[0098] Constructing an edge-weighted graph Each edge weighted graph The vertex set V = {v1,v2,…,v 31} corresponds to the rainfall time series data set TSS = {TS1, TS2, …, TS 31}. Weighted edge set E = {e 1,2 ,e 1,3 ,…,e 30,31 The weight calculation method for each edge of} is: The weighted edges with a distance less than the threshold ζ are retained. In the empirical evaluation, the use of a threshold of 80% works well in practical scenarios. The set of weighted edges after screening is E = {e1, e2, …, e 93}; Community detection using greedy modularity algorithm. For edge-weighted graphs Community after testing

[0099] Record 31 rainfall time series in three edge weighted graphs The co-occurrence of the community. Among them:

[0100]

[0101]

[0102] Determine the weight w of each feature in the feature set Fi . According to the national standard GB / T 28592-2012, rainfall categories are mainly divided into: trace rainfall, light rain, moderate rain, heavy rain, rainstorm, heavy rainstorm and extremely heavy rainstorm. Combined with the knowledge related to floods, the main rainfall of floods is concentrated in the range of moderate rain to heavy rainstorm. Therefore, the number of clusters C is predefined as 4. According to the formula in the technical solution, w1=0.66, w2=1, w3=0.5. Construct a weighted co-occurrence matrix Take time series TS1 and TS3 as examples. Calculate the Euclidean distance between the row vectors of the weighted co-occurrence matrix, that is, The K-Medoid algorithm is used to cluster similar time series and generate the final clustering set CS = {cs1, cs2, cs3, cs4}.

[0103] The rainfall time series cluster cs1 = {TS1, TS3, TS8, TS 15 ,TS 23 ,TS 24} as an example, the six rainfall time series are fully connected, and the adjacent time series are labeled FL = {fl1, fl2, fl3, fl4, fl5} to form the complete rainfall time series TS of cluster cs1. glob ={ts1,ts2,…,ts 273 Generate a complete rainfall time series TS glob With TS glob The DTW distance matrix Each element of the matrix is ​​the Euclidean distance between the values.

[0104] According to relevant experiments of the regional generation network (RPN), the sliding window with a side length ratio of 1:1, 1:2 and 2:1 can detect the target position quickly. Referring to the national standard GB / T 28592-2012, the rainfall levels are divided into 12h rainfall and 24h rainfall. Therefore, the sliding window size is selected as {48×48, 48×24, 24×48, 24×24, 24×12, 12×24, 12×12, 12×6, 6×12}. The upper left corner of the window is used as the reference point, and the direction is from top to bottom and from left to right. According to the parameter sensitivity experiment, the sliding step size st is selected as 5. In one sliding process, the local matrix Mat in the window is calculated dis_local The shortest path distance d sp Taking the example of sliding a window of size 24×12 left and right, the current local distance matrix The reuse domain is the 7th column element The calculation domain is the overlapping area of ​​the first 28 columns and the last 5 columns of elements The same applies to the up and down sliding area allocation. The DTW algorithm is used to calculate the area The cumulative distance of the elements within. Reuse the previous local distance matrix The local shortest path of some areas is obtained The cumulative distance of the elements inside. Continue to use the DTW algorithm to calculate the area Taking the sliding window size from 48×48 to 48×24 as an example, for the sliding window in the same position (i.e. the upper left corner reference point is the same),

[0105] TS 27,48 and TS 91,24 For example, the normalized DTW distance between them is 1.266, and the subsequence pair triple is SPT = <(27,48), (91,24), 1.266>. According to the knowledge in the hydrological field, the similarity threshold τ is set to 0.75. The DTW distance of 1.266 is converted into a similarity of about 0.79. Therefore, this triple can be used as a motif candidate triple CMT. Taking the S601 candidate motif triple CMT = <(27,48), (91,24), 1.266> as an example, the candidate motif list CML needs to add (27,48) and (91,24) at the same time. And take (27,48) and (91,24) as the head nodes of the linked list respectively, and insert (91,24) and (27,48) at the tail node of the current linked list respectively. If a single linked list does not exist, create it directly and add it after the head node. According to all linked lists simlinklist i The length of the candidate motif list CML is sorted, and the first three candidate motifs in the sorted candidate motif list CML are selected as the final motifs. Therefore, there are 12 final motifs in the cluster set CS, and the 12 motifs have different hydrological characteristics. Here, the number of motifs K can be modified according to user needs.

Claims

1. A hydrological time series motif mining method based on numerical feature clustering, characterized in that: The following steps are involved: (1) Extract the numerical features of multiple hydrological time series, select the features that are highly correlated with the time series labels according to the time series labels, and on this basis select the minimum number of features greater than the coverage threshold; (2) Construct an edge-weighted graph with the minimum number of features as the basic unit, determine the retained distance threshold based on the weighted edge screening ratio, and perform initial community detection on the weighted graph of each feature; (3) Calculate feature weights based on the number of communities in each weighted graph and the number of predefined clusters, and construct a weighted co-occurrence matrix; calculate the similarity between row vectors in the co-occurrence matrix to achieve time series clustering set; (4) Taking each clustering set in step (3) as the basic unit, fully connect the time series in the set to form a complete time series, and compare it with itself to form a global distance matrix; (5) Use a sliding window of variable size to slide on the global distance matrix, and calculate the shortest path distance of the local distance matrix within the window during each sliding process; (6) constructing a set of subsequence pair triples based on the shortest path distance, and screening to form a list of candidate motifs, and selecting a predefined number of motifs according to actual applications; The step (5) comprises the following steps: (51) Use windows of different area sizes and different side length ratios to slide on the global distance matrix with a certain step size st; (52) During a sliding process, the local distance matrix Mat within the sliding window is calculated dis_local The regular path distance, that is, the shortest path distance d sp : the sp =Min{DisWarping(Mat dis_local )} (53) gradually reducing the sliding window size according to the side length ratio and area setting, and repeating steps (51) and (52); (54) For each window in step (52) and step (53), the local distance matrix Mat is constructed dis_local Calculate the shortest path distance d sp ; The shortest path distance uses the shortest path distance calculation optimization method, that is, the cumulative distance of the accumulated matrix elements in the reused previous large window and the local shortest path of the matrix in the reused previous neighbor window; the shortest path distance d sp The length needs to be normalized according to the subsequence; The shortest path distance calculation optimization method comprises the following steps: (541) According to the DTW algorithm shortest path solution process, the calculation domain Mat is calculated. dis_local The cumulative distance of all elements in [,:r+1]; (542) Multiplexing Field Mat dis_local The cumulative distance of all elements in [,N-st-1:N-st] is solved, and the previous order distance matrix is ​​reused from the calculation domain Mat dis_local The difference in the cumulative distance from [,:r+1] to the element, that is, in and They represent the previous order cumulative distance matrix and the current cumulative distance matrix respectively; (i, j) and (p, q) are located in the calculation domain and reuse domain of the previous order cumulative matrix respectively; (543) According to the DTW shortest path solving process and the cumulative distance of all elements in the reuse domain known in step (542), the defined computational domain Mat is calculated. dis_local Cumulative distance of all elements in [,N-st:], Mat dis_local The cumulative distance of the elements in the Mth row and the Nth column is the shortest path distance d of the matrix. sp ; The step (6) comprises the following steps: (61) Construct subsequence pair triple SPT = <(idx1,len1),(idx2,len2),d sp >; (idx, len) represents the starting index and length of the subsequence respectively, and dsp is the DTW distance between the subsequence pairs; (62) The subsequence pairs are screened according to the predefined candidate motif similarity threshold τ to form a candidate motif triple CMT, i.e., the normalized shortest path distance sim(d sp )≥τ; where sim is the similarity solution function, and the similarity is inversely proportional to the DTW distance; (63) Select subsequence elements in the candidate motif triple CMT and add them to the candidate motif list CML; (64) Each element of the candidate motif list CML (cmidx i ,cmlen i ) Associate a one-way linked list simlinklist i , each element in the linked list is the candidate motif triple CMT and (cmidx i ,cmlen i ) similar subsequences, that is, satisfying in, is any element of a singly linked list; (65) Sort by the length of the one-way linked list, and according to the number of motifs K defined by the user, select the first K elements of the sorted candidate motif list CML, which is the K-motif.

2. The hydrological time series motif mining method based on numerical feature clustering according to claim 1 is characterized in that: The step (1) comprises the following steps: (11) Form relevant numerical features based on the duration, amplitude and trend of the time series of interest in the hydrological field; (12) Calculate the relevance of features based on the labels of the time series, i.e., the p-value, and select the top features after sorting them by p-value. f features; to reduce feature redundancy, select the minimum number of features that meets the coverage δ.

3. The hydrological time series motif mining method based on numerical feature clustering according to claim 1 is characterized in that: The step (2) comprises the following steps: (21) Construct edge weighted graph G w :For a feature f i And the time series set TSS = {TS1, TS2, ..., TS m }, each time series TS p For edge weighted graph A node v of the vertex set V p ; Each edge e in the weighted edge set E p,q Corresponding to a weight Keep the edges whose distance is less than the threshold ζ; (22) A community detection algorithm is applied to assign several nodes to different communities in each weighted graph. The community detection algorithm uses the connection structure of the graph to unify different representations of time series clustering.

4. The hydrological time series motif mining method based on numerical feature clustering according to claim 1 is characterized in that: The step (3) comprises the following steps: (31) The weight of each feature is determined by calculating the ratio between the number of user-defined clusters and the number of communities. The specific calculation method is: Where C represents the number of predefined clusters, O i Representing a weighted graph the number of communities in (32) Count the number of times each pair of time series appears in the same community and fill in the co-occurrence weighted matrix The specific method is: in, is the co-occurrence weighted matrix, Represents the time series TS i and time series TD j The weighted co-occurrence ratio between ; the calculation formula is: Among them, w k Represents the weight corresponding to any feature, w k' Indicates TS i and TD j The weight corresponding to any feature that appears in a community, m and m' represent the total number of features and TS respectively. i and TD j The total number of times they co-occur in a community; (33) Yes The distance between rows is calculated, and the clustering algorithm is used to form the final clustering set.

5. The hydrological time series motif mining method based on numerical feature clustering according to claim 1 is characterized in that: The distance matrix in step (4) is: In a single clustering set cs i In the above example, all time series in the set are arranged in chronological order to form a complete time series TS. glob ={ts1,ts2,…,ts m }; Complete time series TS glob Each element in TS glob Calculate the Euclidean distance for each element in and initialize the global distance matrix Mat dis_glob : in, It represents the distance between the sequence values ​​of the time series at time i and time j, and n represents the length of the time series.

Images