Abnormity detection method and device based on dynamic graph bidirectional transformation

Abstract

The invention discloses an anomaly detection method and device based on dynamic graph bidirectional transformation, and the method comprises the steps: analyzing a brand new feature representation learning algorithm through dynamic graph bidirectional edge transformation, and comprehensively capturing hidden spatial-temporal features which are beneficial to describing the data change of a normal operation state in industrial process data, thereby improving the anomaly detection effect; and the possibility of abnormal report omission is reduced. The anomaly detection method disclosed by the invention not only discloses dynamic graph bidirectional edge transformation analysis and an improved technical scheme thereof, but also provides three improved schemes for determining spatial neighbors, so that change characteristics of learning data in time sequence dynamics and spatial distribution are effectively represented. Besides, through establishment of hardware and software, an anomaly detection device based on dynamic graph bidirectional transformation is established, and the anomaly detection device comprises a parameter initialization module, a data acquisition module, a data processing module, a display module and an anomaly alarm so as to execute the implementation process of the disclosed anomaly detection method.

Description

Technical Field

[0001] The present invention relates to the technical field of data-driven anomaly detection, and in particular to an anomaly detection method and device based on dynamic graph bidirectional transformation. Background Art

[0002] Real-time monitoring of the operating status of industrial processes and timely detection of anomalies or failures in operation are the key to ensuring production safety and cost optimization control. Especially with the continuous expansion and increasing complexity of modern industrial processes, effective anomaly detection technology is increasingly needed to monitor production safety and product quality stability. Considering that complex industrial processes usually have a large number of sensors measuring process data in real time, as well as the complexity and uncertainty of large-scale modern industrial processes, it is almost impossible to establish an accurate mechanism model for anomaly detection. In contrast, it has become increasingly feasible to use a large amount of process data to mine useful information for detecting abnormal operating conditions, and it has also received attention from the industry.

[0003] An English paper titled "Density-Based Structure Preserving Projections Process Monitoring Model for Fused Magnesia Smelting Process" published in the journal "IEEE Transactions on Industrial Informatics" in 2023 proposed a process monitoring method based on a density structure preserving projection algorithm. By learning the representation of the density information of the spatial distribution of process data, it realized the anomaly detection of the magnesium smelting process. Another English paper titled "Multineighbor Ensemble Preserving Embedding for Process Monitoring" published in the journal "IEEE Transactions on Instrument and Measurement" in 2024 disclosed a process monitoring method based on a multi-neighbor ensemble preserving embedding algorithm. By learning the representation of the multi-faceted neighbor structure information of the process data in the spatial distribution, the sensitivity of anomaly detection in the magnesium smelting process was further improved. However, these methods do not take into account the temporal change information of the process data in the sampling time, that is, the dynamic nature of the process data, and the effectiveness of its anomaly detection needs to be further improved.

[0004] The Chinese invention patent with application number 202310063661.8 discloses a multi-scale process fault detection method based on a full-structure dynamic autoregressive latent variable model. This method implements representation learning of time series change information from the perspective of the full structure of process data, establishes a corresponding autoregressive latent variable model, and thus implements abnormality detection of chemical processes. The Chinese invention patent with application number 202410177636.7 discloses an industrial process fault detection method and device with enhanced spatiotemporal neighborhood features of data. It uses spatiotemporal-stacked sparse autoencoder-canonical variable analysis to implement representation learning of spatial distribution and time series change information of process data, solving the problem that traditional process monitoring methods ignore the spatiotemporal neighborhood of data and do not consider the dynamic nature of data, resulting in insufficient fault detection efficiency. In addition, the Chinese invention patent with application number 202410878009.6 discloses a fault detection method based on a dynamic latent variable neighborhood preservation embedding algorithm. By sequentially implementing dynamic analysis and neighbor preservation embedding analysis on process data, representation learning of the change characteristics of process data is achieved, and the sensitivity of abnormality detection for semi-intermittent production processes is improved.

[0005] It can be seen that, on the basis of learning the representation of spatial distribution features of process data, further considering the dynamics of process data is conducive to improving the anomaly detection effect. However, the technical solutions disclosed in these three invention patents all perform feature representation learning in sequence in a serial manner, either first for spatial distribution features and then for dynamic features in time series; or vice versa, first performing dynamic analysis and then performing spatial distribution feature analysis. It is impossible to coordinately consider the dynamics of process data in time series and the neighbor relationship in spatial distribution when learning feature representation. From this perspective, these methods cannot comprehensively and effectively implement feature representation learning of process data, which may lead to anomalies or insufficient fault detection rates in certain faults. Summary of the invention

[0006] The main technical problem to be solved by the present invention is: how to collaboratively represent the changing characteristics of learning process data in terms of temporal dynamics and spatial distribution, and comprehensively capture the characteristics hidden in industrial process data that are beneficial to describing the changes in normal operating status data, thereby improving the effect of anomaly detection and reducing the possibility of missed anomalies.

[0007] The technical solution adopted by the present invention to solve the above-mentioned problem is: an anomaly detection method and device based on bidirectional transformation of dynamic graphs; wherein, the implementation process of an anomaly detection method based on bidirectional transformation of dynamic graphs disclosed by the present invention includes steps 1 to 4 shown below.

[0008] Step 1: When the industrial process is running normally, obtain process data at N consecutive sampling times and arrange them into N column vectors x1, x1, …, x2 in the order of sampling time.N , and then according to the formula X i =[x i+d ,x i+d-1 ,…,x i ] Build n dynamic graph matrices X1, X2, ..., X n After that, each dynamic graph matrix is ​​standardized; among them, X i represents the i-th dynamic graph matrix, N represents an integer not less than 500, the dynamic order d is an integer not less than 2 and not greater than 10, n=Nd, and number i=1,2,…,n.

[0009] Step 2: Determine the k nearest neighbors corresponding to each of the normalized dynamic graph matrices, and calculate the average value of these k nearest neighbors, which are recorded as Here, k is an integer not less than 5 and not more than 10.

[0010] Step 3: Perform bidirectional transformation analysis on the dynamic graph to obtain J left transformation vectors u1, u2, ..., u J and J right transformation vectors v1,v2,…,v J Then, through the formula Calculate n score vectors s1,s2,…,s n ; Among them, s i represents the i-th score vector, and the parameter J is a smaller than the dynamic graph matrix X i The positive integer is the total number of row vectors in . The superscript T indicates the transpose of the matrix or vector.

[0011] Step 4: According to the formula After calculating the abnormal index corresponding to each score vector and determining the threshold of the abnormal index, the following steps 4-1 to 4-3 are repeatedly performed according to the sampling time interval; wherein, ζ i Indicates i The corresponding abnormal index, Λ represents the covariance matrix of the score vector, that is: Λ=[s1,s2,…,s J ][s1,s2,…,s J ] T / (n-1).

[0012] Step 4-1: Obtain the process data of the industrial process at the latest sampling time and form the corresponding column vector x t , and then according to the formula X t =[x t ,x t-1 ,x t-2 ,…,x t-d ] Build the dynamic graph matrix X of the latest sampling time t t After that, X tPerform the same normalization as in step 1; where x t-1 ,x t-2 ,…,x t-d Column vector representing the d consecutive sampling times before t.

[0013] Step 4-2: Calculate the score vector s in the same way as in step 3 t Then, according to the formula Calculate the abnormal index ζ at the latest sampling time t , and determine whether it is greater than the threshold determined in step 4; if so, execute step 4-3 to decide whether to trigger an abnormal alarm; if not, return to step 4-1.

[0014] Step 4-3, repeat steps 4-1 to 4-2 at the sampling time interval until abnormal indicators of f consecutive latest sampling moments are obtained. If they are all greater than the threshold determined in step 4, an abnormal alarm is triggered; where f is a positive integer not less than 3.

[0015] In the above implementation steps, the core technical feature of the method of the present invention that distinguishes it from the prior art solution lies in the dynamic graph bidirectional transformation analysis in step 3, including steps 3-1 to 3-4 shown below.

[0016] Step 3-1: After determining the parameter J, set j=1, 2, ..., J in sequence, and execute the following step 3-2 to obtain the corresponding j-th left transformation vector u j and the jth right transformation vector v j .

[0017] Step 3-2: Randomly initialize u j , and then repeat the steps 3-3 to 3-4 shown below until u j convergence.

[0018] Step 3-3: Solve the generalized eigenvalue problem The smallest eigenvalue λ j After the corresponding eigenvector v, Calculate the jth right transformation vector v j .

[0019] Step 3-4: If j is equal to 1, solve the generalized eigenvalue problem Gu = λ j The smallest eigenvalue λ in Lu j The corresponding eigenvector u; if j is greater than 1, solve the eigenvalue problem The eigenvector u corresponding to the minimum eigenvalue in Update j ;in, I represents the identity matrix, Uj-1 =[u1,u2,…,u j-1 ],

[0020] The dynamic graph matrix X involved in the method of the present invention i It is to integrate the d+1 column vectors with continuous sampling time, and then analyze them through the dynamic graph bidirectional transformation, that is, through the left and right bidirectional transformation of u and v, that is: u T X i v, to mine each dynamic graph matrix X1,X2,…,X n The spatial neighbor relationship characteristics between them can be described by mathematical formula:

[0021] in, It can be seen that the constrained minimization problem defined in the above formula ① aims to minimize u T X i v and The distance between Represents X i Therefore, the above formula ① aims to minimize the distance between each dynamic graph matrix and its neighbor center through the left and right bidirectional transformation of u and v, so as to maintain the original spatial neighbor distribution relationship of the dynamic graph matrix in the feature extraction analysis.

[0022] The above equation ① can be solved by the Lagrange multiplier method. First, the Lagrange function l as shown below must be constructed:

[0023] When the partial differential equations of l with respect to u, v and λ are all equal to 0, the objective function in the above formula ① takes the minimum value, that is:

[0024] Formula ⑤ is equivalent to the constraint condition in formula ①. Formula ④ can directly lead to the generalized eigenvalue problem in step 3-3, while formula ③ can directly lead to the generalized eigenvalue problem Gu=λLu in step 3-4. Since formula ① is for obtaining the minimum value, the minimum eigenvalue is calculated in both steps 3-3 and 3-4.

[0025] When j = 1, the first left transformation vector u1 and the first right transformation vector v1 are obtained by repeatedly solving the generalized eigenvalue problem in steps 3-3 and 3-4 until convergence; when more left transformation vectors and right transformation vectors need to be solved, that is, when j>1, a constraint condition can be added to the above formula ①, and then the above formula ① is transformed into the minimization problem as shown below:

[0026] Similarly, use the Lagrange multiplier method to first construct the Lagrange function l j :

[0027] Among them, η j-1 ,…,η1 is expressed as a constraint condition The Lagrange multiplier introduced. Then, let l j Relative to u j The partial differential equation is equal to 0, and we can get:

[0028] Multiply both sides of the equal sign of the above formula ⑧ on the left Then we can get j-1 sets of equations, namely:

[0029] The above formula ⑨ can be equivalently written in matrix form: Where η=[η1,…,η j-1 ] T Multiply the above formula ⑧ by L on the left -1 Then substitute You can infer the eigenvalue problem that needs to be solved in step 3-4

[0030] It should be emphasized that the dynamic graph bidirectional transformation analysis involved in the present invention is a brand-new data feature analysis algorithm, and the idea and reasoning process of its algorithm design have been shown in formulas ① to ⑨. From the above reasoning process, it can be seen that the algorithm takes into account the temporal dynamics of process data by constructing a dynamic graph matrix, and then performs left and right bidirectional transformations on the dynamic graph matrix to mine and analyze the neighbor relationship characteristics of the dynamic graph matrix in spatial distribution. Therefore, the dynamic graph bidirectional transformation analysis involved in the present invention can collaboratively represent the changing characteristics of learning process data in temporal dynamics and spatial distribution, thereby comprehensively capturing the characteristics hidden in industrial process data that are conducive to describing the changes in normal operating status data.

[0031] As an improvement of the method of the present invention, the above step 3-1 is preferably replaced by: after setting the parameter J, set j=1, 2, ..., M in sequence, and execute the above step 3-2 to obtain the corresponding j-th left transformation vector u j , the jth right transformation vector v j , and the minimum eigenvalue λ j ; Wherein, M is equal to the dimension of the column vector in step 1; and then after step 3-4, add step 3-5 as shown below.

[0032] Step 3-5: Replace λ1,λ2,…,λ M After arranging in ascending order, the left transformation vectors and right transformation vectors corresponding to the first J smallest eigenvalues ​​are used as the J left transformation vectors and J right transformation vectors obtained in the above step 3, respectively.

[0033] Although in the implementation, the minimum eigenvalues ​​λ1,λ2,…,λ M are obtained in sequence, but λ1,λ2,…,λ M The arrangement order is not necessarily strictly in ascending order. Therefore, the improved technical solution for the dynamic graph bidirectional transformation analysis is to ensure that the final J left transformation vectors and J right transformation vectors are selected to correspond to the J smallest eigenvalues ​​in terms of value, so as to ensure that they absolutely meet the minimization target in the above formula ① or ⑥, so as to ensure that the dynamic graph bidirectional transformation analysis in step 3 can play the best role in feature analysis and extraction.

[0034] As an improvement of the method of the present invention, step 2 is as follows: Steps A-1 to A-3 are as follows: i Determine its corresponding k nearest neighbors.

[0035] Step A-1: ​​According to the formula Calculate X i The distance between the remaining n-1 dynamic graph matrices; where X m Represents the mth dynamic graph matrix, numbered m=1,2,…,n and m≠i, d i,m Represents X i With X m The distance between Indicates calculation of X i -X m The sum of the squares of all elements in .

[0036] Step A-2: For any given number p=1, 2, ..., n not equal to i and m, determine d i,p +d p,m The minimum value that can be obtained, and then update d according to the minimum value i,m , the update rule is: if the minimum value is less than d i,m , then update di,m is equal to the minimum value; if not, keep d i,m unchanged; among them, d i,p Represents X i The distance between the pth dynamic graph matrix, d p,m Represents the pth dynamic graph matrix and X m The distance between.

[0037] Step A-3: According to the updated distance, i The k dynamic graph matrices with the smallest distance between them are determined as X i The corresponding k nearest neighbors.

[0038] Generally speaking, through The calculated distance can be directly determined as X i Determine the corresponding k nearest neighbors, which is also the most common way to determine the nearest neighbors in the technical field. The implementation process of the above steps A-1 to A-3 further updates the calculated distance on this basis. The updated distance is called geodesic distance in the technical field, that is, the shortest path distance. This is why d in step A-2 i,p +d p,m The purpose of taking the minimum value is to find the shortest path. The advantage of using geodesic distance is that it can better describe the positional relationship of data spatial distribution.

[0039] As another improvement of the method of the present invention, step 2 is as follows: step B-1 to step B-2 are as follows: i Determine its corresponding k nearest neighbors.

[0040] Step B-1, randomly initialize the feature vector u, and then repeatedly perform steps B-1.1 to B-1.2 as shown below until u converges, according to the formula Calculate X i With X m The difference distance between Among them, X m Represents the mth dynamic graph matrix, numbered m=1,2,…,n and m≠i.

[0041] Step B-1.1. Solve the generalized eigenvalue problem After the eigenvector v corresponding to the maximum eigenvalue η in Update v.

[0042] Step B-1.2: Solve the generalized eigenvalue problem After the eigenvector u corresponding to the maximum eigenvalue η in Update u.

[0043] Step B-2: Repeat step B-1 until you get X i After the difference distance between the other n-1 dynamic graph matrices, i The k dynamic graph matrices with the smallest difference distance are determined as X i The corresponding k nearest neighbors.

[0044] The solutions to the two generalized eigenvalue problems in the above steps B-1.1 and B-1.2 are derived from the following optimization problem:

[0045] The optimization problem aims to maximize ||u T X m v|| 2 while minimizing ||u T X i v|| 2 , and thus identify X through u and v i With X m The maximum difference between them is then used to determine X i The corresponding neighbors.

[0046] The solution of the above ⑩ can be achieved by the Lagrange multiplier method, which can be set v T X i T uu T X i v = 1, thus the above formula ⑩ is equivalently transformed into an optimization problem with constraints as shown below:

[0047] Then, construct the Lagrangian function And set The partial differential with respect to v is equal to 0, that is:

[0048] From the above formula It can be inferred that This is the generalized eigenvalue problem in step B-1.1. Similarly, the Lagrangian function can also be constructed Reset The partial differential with respect to u is equal to 0, and we get the generalized eigenvalue problem in step B-1.2.

[0049] The difference distance calculated by maximizing the difference between the two dynamic graph matrices can fully explore the difference between the two dynamic graph matrices, so that the k dynamic graph matrices with the smallest difference distance can be determined as X i The corresponding k nearest neighbors. In terms of the way of calculating distance, the improved technical solution of using difference distance to determine the nearest neighbors can fully explore the change differences between dynamic graph matrices, thereby better describing the spatial distribution relationship of data.

[0050] As a further improvement of the method of the present invention, step 2 is as follows: Steps C-1 to C-6 are as follows: i Determine its corresponding k nearest neighbors

[0051] Step C-1, according to the formula Calculate the distance between n dynamic graph matrices; where X m Represents the mth dynamic graph matrix, numbered m=1,2,…,n and m≠i, d i,m Represents X i and the mth dynamic graph matrix X m The distance between.

[0052] Step C-2: After setting two width parameters β1 and β2 and the skewness parameter ε, the basis function values ​​between each of the n dynamic graph matrices are calculated according to the formula shown below.

[0053] Among them, f i,m Represents X i With X m The basis function values ​​between β1 and β2 are both greater than zero, ε is not greater than 0.001 and greater than zero, Represents X i The k-order nearest neighbor distance, exp() represents an exponential function with a natural constant as the base.

[0054] Step C-3, according to the formula Calculate each element in the transfer matrix P; where P i,m represents the element in the i-th row and m-th column of P, and ζ represents the i The numbers of the k dynamic graph matrices with the smallest distance between them.

[0055] Step C-4: After randomly initializing the feature vector α, continuously update α according to the formula α=αP until it converges, then merge n repeated α to form the feature matrix A, and calculate the basic matrix B=(I n -P+A) -1 .

[0056] Step C-5: According to the formula Calculate X i With X m The similarity S between i,m , where C i,m =α m Bi,m +α i B m,i -α i α m -α i δ i,m , α i and α m denote the i-th and m-th elements in α, respectively, and B i,m represents the element in the i-th row and m-th column of B. When i=m, δ i,m =1, other casesδ i,m =0.

[0057] Step C-6, and X i The k dynamic graph matrices with the greatest similarity are determined as X i The corresponding k nearest neighbors.

[0058] The k-order nearest neighbor distance in step C-2 above refers to the distance between the remaining n-1 dynamic graph matrices and X i After the distances are arranged in ascending order, the kth smallest distance is defined as X i The k-order nearest neighbor distance.

[0059] The above steps C-1 to C-6 determine the nearest neighbors by means of similarity, and the method of calculating the similarity uses a Markov chain, which can further take into account the temporal relationship between dynamic graph matrices with continuous sampling time, and step C-1 starts with the calculation of spatial distance. Therefore, the similarity determines the nearest neighbors from both time and space perspectives, and more comprehensively considers the temporal and spatial distribution of the data, which is more conducive to subsequent feature extraction.

[0060] As a further improvement of the present invention, in the above step C-2, β1 can be set equal to the average value of the distances between the n dynamic graph matrices calculated in step C-1.

[0061] The above preferred technical solution for setting the width parameter β1 provides an objective parameter setting method, which can avoid underfitting and overfitting problems caused by setting β1 too large or too small.

[0062] In addition, considering the average value of k nearest neighbors It is the key to describe the spatial distribution of data. The present invention preferentially uses the weighted average of k nearest neighbors and provides the following two methods for calculating the weighted coefficients.

[0063] First, the weighting coefficient is determined by the softmax function, and the corresponding weighted average calculation process includes steps D-1 to D-2 shown below.

[0064] Step D-1: X i The k nearest neighbors of Then, according to the formula Calculate X i and its cth nearest neighbor The weighting coefficient between ; Among them, the width parameter σ i Equal to X i The average distance between a value and its k nearest neighbors, where the nearest neighbors are numbered c = 1, 2, …, k.

[0065] Step D-2: According to the formula calculate

[0066] Secondly, the weighting coefficient is determined by the least squares regression algorithm, and the corresponding weighted average calculation process includes the following steps E-1 to E-2:

[0067] Step E-1: X by column i Expand into column vector z i , and then in the same way Expand into corresponding column vectors, recorded as in, Represents X i The k nearest neighbors of .

[0068] Step E-2: Combine to form the input matrix Then, through the formula Calculate the weight coefficient vector w i .

[0069] Step E-3: According to the formula calculate in, Indicates w i The first to kth elements in are the weighting coefficients used to calculate the weighted average.

[0070] The previous method of calculating the weighted average further takes into account X i and its k nearest neighbors The difference in distance between neighbors, the larger the distance, the smaller the weight coefficient. The latter method of calculating the weighted average value uses the least squares regression fitting technique to determine the weight coefficient, which is determined from the perspective of neighbor reconstruction. Both methods can convert X i The spatial analysis differences between the k nearest neighbors are taken into account, which is beneficial to the dynamic graph bidirectional transformation analysis in the subsequent step 3.

[0071] Based on the same inventive concept, the present invention also provides an anomaly detection device based on bidirectional transformation of dynamic graphs, including a parameter initialization module, a data acquisition module, a data processing module, a display module, and an anomaly alarm; the functions of each component are as follows.

[0072] The parameter initialization module provides the data processing module with the mean and standard deviation required for standardization, as well as the J left-transformed vectors u1, u2, …, u required for calculating the anomaly detection index. J , J right transformation vectors v1,v2,…,v J and covariance matrix Λ, and provide the display module with a threshold for judging whether it is abnormal.

[0073] The data acquisition module dynamically stores the latest process data of d consecutive sampling moments, and can obtain the latest process data of the sampling moment in real time according to the sampling time interval, and send it together with the dynamically stored process data of d consecutive sampling moments to the data processing module in real time.

[0074] After the data processing module receives the process data of d+1 sampling moments, it constructs the dynamic graph matrix X t , and then use the mean value and standard deviation provided by the parameter initialization module to standardize it, and then use the J left transformation vectors u1,u2,…,u provided by the parameter initialization module J , J right transformation vectors v1,v2,…,v J and covariance matrix Λ, calculate the corresponding anomaly detection index ζ t , and finally t The specific value is sent to the display module.

[0075] The display module receives the abnormal detection indicator ζ sent by the data processing module t After obtaining the specific value of the anomaly detection indicator, the corresponding scatter points are plotted with the specific value of the anomaly detection indicator as the ordinate and the latest sampling time as the abscissa. At the same time, the threshold provided by the parameter initialization module is displayed as a horizontal line. Then, it is determined whether the anomaly detection indicators of the consecutive f latest sampling times are greater than the threshold; if so, the anomaly alarm is turned on; if not, the anomaly alarm is turned off; wherein f is a positive integer not less than 3.

[0076] The core of the anomaly detection device is that the parameter initialization module will immediately obtain the process data of N consecutive sampling moments of normal operation from the historical database of the industrial process according to the initialization instruction, firstly form n dynamic graph matrices, and calculate the average value and standard deviation of these n dynamic graph matrices; secondly, after standardizing each dynamic graph matrix using the calculated average value and standard deviation, determine its corresponding k nearest neighbors; then, implement the dynamic graph bidirectional transformation analysis according to steps 3-1 to 3-4, and obtain J left transformation vectors u1, u2, ..., u J , J right transformation vectors v1,v2,…,v J , and the covariance matrix Λ; Finally, calculate the abnormal index corresponding to each score vector and determine the threshold of the abnormal index.

[0077] It can be seen that the parameter initialization module executes the implementation process of steps 1 to 4 in the above-mentioned anomaly detection method, while the data processing module mainly executes the implementation process of steps 4-1 to 4-4; other modules involve functions such as real-time data acquisition, abnormal indicator display and triggering of abnormal alarms.

[0078] Therefore, the anomaly detection device disclosed in the present invention and the anomaly detection method disclosed previously are both based on the same inventive concept of bidirectional transformation of dynamic graphs.

[0079] In addition, the present invention also provides an improved anomaly detection device, in which the parameter initialization module replaces step 3-1 with setting parameter J, sets j=1, 2, ..., M in sequence, and executes step 3-2 to obtain the corresponding jth left transformation vector u j , the jth right transformation vector v j , and the minimum eigenvalue λ j ; and add the implementation process of step 3-5 after step 3-4.

[0080] This improved anomaly detection device implements an improved technical solution for the dynamic graph bidirectional transformation analysis involved in the above-mentioned anomaly detection method to ensure that the J left transformation vectors and J right transformation vectors obtained by parameter initialization can play the best role in feature analysis and extraction. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] Figure 1 This is the process flow chart of the sewage treatment process.

[0082] Figure 2 The present invention is a flowchart for implementing the bidirectional transformation analysis of dynamic graphs.

[0083] Figure 3 The figure is a schematic diagram of online detection of the abnormality detection method involved in the present invention.

[0084] Figure 4 The present invention is a schematic diagram for implementing the improved solution for the bidirectional transformation analysis of dynamic graphs.

[0085] Figure 5 Schematic diagram of online detection of anomaly detection based on improved dynamic graph bidirectional transformation analysis.

[0086] Figure 6 Schematic diagram of three improved neighbor determination process implementations of the anomaly detection method involved in the present invention.

[0087] Figure 7 The present invention is a schematic diagram of the structure of an abnormality detection device. DETAILED DESCRIPTION

[0088] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0089] In order to clearly and specifically describe the specific implementation of the abnormality detection method according to the present invention, Figure 1 The industrial objects shown are used for anomaly detection. Figure 1 The industrial process shown in the figure is a sewage treatment process, which consists of five interconnected biochemical reaction tanks and a secondary sedimentation tank; the first two reaction tanks operate under anoxic conditions, and the last three reaction tanks operate in the aeration zone. PI represents a proportional-integral controller, which is used to control the process operation.

[0090] Figure 1 The sewage treatment process shown can measure the process data of 7 variables at each sampling moment. Using the process data of N = 672 consecutive sampling moments under normal operating conditions, a corresponding anomaly detection model is first established, and then anomaly detection is performed on the process data at the 2016 latest sampling moments, including steps 1 to 4 shown below.

[0091] Step 1: When the sewage treatment process is running normally, obtain the process data of N = 672 consecutive sampling moments and arrange them into N column vectors x1, x1, ..., x2 in the order of sampling time. N , and then according to the formula X i =[x i+d ,x i+d-1 ,…,x i ] Build n dynamic graph matrices X1, X2, ..., X n Afterwards, each dynamic graph matrix is ​​standardized, wherein the dynamic order d is set to 2.

[0092] Step 2: Determine the k nearest neighbors corresponding to each of the normalized dynamic graph matrices, and calculate the average value of these k nearest neighbors, which are recorded as Here, k is equal to 5.

[0093] Step 3: Follow the steps below: Figure 2 The process shown in the figure implements the dynamic graph bidirectional transformation analysis and obtains J left transformation vectors u1,u2,…,u J and J right transformation vectors v1,v2,…,v J Then, through the formula Calculate n score vectors s1,s2,…,s n , and then calculate the covariance matrix Λ=[s1,s2,…,s J ][s1,s2,…,s J ] T / (n-1); where Figure 2 The process shown in includes steps 3-1 to 3-4 shown below.

[0094] Step 3-1, after setting the parameter J = 4, set j = 1, 2, ..., J in sequence, and execute the following step 3-2 to obtain the corresponding j-th left transformation vector u j and the jth right transformation vector v j .

[0095] Step 3-2: Randomly initialize u j , and then repeat the steps 3-3 to 3-4 shown below until u j convergence.

[0096] Step 3-3: Solve the generalized eigenvalue problem The smallest eigenvalue λ j After the corresponding eigenvector v, Calculate the jth right transformation vector v j .

[0097] Step 3-4: If j is equal to 1, solve the generalized eigenvalue problem Gu = λ j The smallest eigenvalue λ in Lu j The corresponding eigenvector u; in other cases, solve the eigenvalue problem The eigenvector u corresponding to the minimum eigenvalue in Update j .

[0098] Step 4: According to the formula After calculating the abnormal index corresponding to each score vector and determining the threshold of the abnormal index, steps 4-1 to 4-3 shown below are repeatedly performed according to the sampling time interval.

[0099] Step 4-1: Obtain the process data of the industrial process at the latest sampling time and form the corresponding column vector x t , and then according to the formula X t=[x t ,x t-1 ,x t-2 ,…,x t-d ] Build the dynamic graph matrix X at the latest sampling time t t After that, X t The same normalization process as in step 1 was performed.

[0100] Step 4-2, follow Calculate the score vector s t Then, according to the formula Calculate the abnormal index ζ at the latest sampling time t , and determine whether it is greater than the threshold determined in step 4; if so, execute step 4-3 to decide whether to trigger an abnormal alarm; if not, return to step 4-1.

[0101] Step 4-3: Repeat steps 4-1 to 4-2 at sampling time intervals until abnormal indicators at consecutive f=3 latest sampling moments are obtained. If they are all greater than the threshold determined in step 4, an abnormal alarm is triggered.

[0102] Plot all the abnormal indicators of the latest sampling time in 2016 on Figure 3 At the same time, the threshold value determined in step 4 is also drawn as a horizontal line on Figure 3 in; from Figure 3 It can be found that from the 673 latest sampling moments, the abnormal indicators are all greater than the threshold, which can trigger an abnormal alarm.

[0103] In addition, in this embodiment, when the threshold is determined in step 4, the probability density of the n=Nd=670 abnormal indicators obtained in step 4 is estimated by using the kernel density estimation method, and then the value under the condition that the confidence limit is equal to 99% is taken as the determined threshold. In addition, the maximum value of the 670 abnormal indicators can also be directly determined as the threshold, thereby maximizing the reduction of abnormal false alarms.

[0104] The second embodiment of the present invention utilizes the improved dynamic graph bidirectional transformation analysis in step 3, and its implementation process is as follows: Figure 4 As shown, it includes steps 3A to 3E shown below.

[0105] Step 3A, after setting the parameter J = 4, set j = 1, 2, ..., 7 in sequence, and execute step 3B to obtain the corresponding j-th left transformation vector u j , the jth right transformation vector v j , and the minimum eigenvalue λ j .

[0106] Step 3B: Randomly initialize u j, and then repeat steps 3C to 3D as shown below until u j convergence.

[0107] Step 3C: Solve the generalized eigenvalue problem The smallest eigenvalue λ j After the corresponding eigenvector v, Calculate the jth right transformation vector v j .

[0108] Step 3D: If j is equal to 1, solve the generalized eigenvalue problem Gu = λ j The smallest eigenvalue λ in Lu j The corresponding eigenvector u; in other cases, solve the eigenvalue problem The eigenvector u corresponding to the minimum eigenvalue in Update j .

[0109] Step 3E: Replace λ1,λ2,…,λ M Arrange in ascending order, and then use the left transformation vectors and right transformation vectors corresponding to the first J smallest eigenvalues ​​as the J left transformation vectors and J right transformation vectors obtained in step 3, respectively.

[0110] The anomaly detection schematic diagram corresponding to the second embodiment of the present invention is as follows Figure 5 shown; Figure 3 Compared with the above, after the anomaly occurs, its anomaly detection index becomes more sensitive because Figure 5 The difference between the anomaly detection indicator in and the threshold represented by the horizontal line is more obvious since the 673 latest sampling moments.

[0111] Figure 5 In relative to Figure 3 The advantage of the proposed method is mainly due to the improved bidirectional transformation analysis of dynamic graphs. The left transformation vector and the right transformation vector corresponding to the four smallest eigenvalues ​​are selected for feature analysis and extraction.

[0112] The third embodiment of the present invention improves the neighbor determination method based on the first two embodiments, involving three different improved methods, respectively: Figure 6 As shown in the three sub-graphs numbered (a), (b) and (c), (a) determines the nearest neighbors according to the geodesic distance between the dynamic graph matrices, (b) determines the nearest neighbors according to the maximized difference distance, and (c) determines the nearest neighbors through the similarity between the dynamic graph matrices based on the idea of ​​Markov chain.

[0113] If you use Figure 3 The nearest neighbor determination method shown in (c) is a key parameter β1 that can be set equal to Figure 3The average distance between the n dynamic graph matrices calculated in step C-1.

[0114] In addition, through the construction of hardware and software, an anomaly detection device based on dynamic graph bidirectional transformation is established, and its structure is as follows Figure 7 As shown, it includes a parameter initialization module, a data acquisition module, a data processing module, a display module, and an abnormal alarm. The functions of each module are introduced as follows.

[0115] The parameter initialization module provides the data processing module with the mean and standard deviation required for standardization, as well as the J left-transformed vectors u1, u2, …, u required for calculating the anomaly detection index. J , J right transformation vectors v1,v2,…,v J and covariance matrix Λ, and provide the display module with a threshold for judging whether it is abnormal.

[0116] The data acquisition module dynamically stores the latest process data of d consecutive sampling moments, and can obtain the latest process data of the sampling moment in real time according to the sampling time interval, and send it together with the dynamically stored process data of d consecutive sampling moments to the data processing module in real time.

[0117] After the data processing module receives the process data of d+1 sampling moments, it constructs the dynamic graph matrix X t , and then use the mean value and standard deviation provided by the parameter initialization module to standardize it, and then use the J left transformation vectors u1,u2,…,u provided by the parameter initialization module J , J right transformation vectors v1,v2,…,v J and covariance matrix Λ, calculate the corresponding anomaly detection index ζ t , and finally t The specific value is sent to the display module.

[0118] The display module receives the abnormal detection indicator ζ sent by the data processing module t After the specific value of the anomaly detection indicator is obtained, the corresponding scatter points are plotted with the specific value of the anomaly detection indicator as the ordinate and the latest sampling time as the abscissa. At the same time, the threshold provided by the parameter initialization module is displayed as a horizontal line. The corresponding displayed image is Figure 3 and Figure 5 The same; then determine whether the anomaly detection indicators at the consecutive f latest sampling moments are all greater than the threshold; if so, turn on the anomaly alarm; if not, turn off the anomaly alarm.

[0119] The above is the hardware and its functional implementation, while the key to the software part is to realize the dynamic graph bidirectional transformation analysis of the anomaly detection method provided by the present invention; Figure 7As shown in FIG. 1 , according to the initialization instruction, the parameter initialization module will immediately obtain the process data of N consecutive sampling moments of normal operation from the historical database of the industrial process, and firstly form n dynamic graph matrices and calculate the average value and standard deviation of these n dynamic graph matrices; secondly, after standardizing each dynamic graph matrix using the calculated average value and standard deviation, the k nearest neighbors corresponding to each of them are determined; then, according to steps 3-1 to 3-4, the dynamic graph bidirectional transformation analysis is performed to obtain J left transformation vectors u1, u2, …, u J , J right transformation vectors v1,v2,…,v J , and the covariance matrix Λ; Finally, calculate the abnormal index corresponding to each score vector and determine the threshold of the abnormal index.

[0120] It can be seen that the mean value and standard deviation, the left and right transformation vectors, and the threshold value provided by the parameter initialization module to the data processing module, as well as the threshold value provided to the display module, all come from the process executed after the initialization instruction.

[0121] Furthermore, when designing and implementing the execution program of the parameter initialization module, the improved dynamic graph bidirectional transformation analysis used in the second embodiment is executed, that is, the dynamic graph bidirectional transformation analysis is implemented according to steps 3A to 3E.

[0123] The above explanation of the technical solution disclosed in the present invention in combination with the accompanying drawings is illustrative rather than restrictive. The technical solution disclosed in the present invention is not limited to the above-mentioned implementation cases. Any modifications made without departing from the purpose of the technical solution disclosed in the present invention and the scope of protection of the claims shall fall within the protection scope of the present invention.

Claims

1. An anomaly detection method based on bidirectional transformation of dynamic graphs, comprising the following steps 1 to 4: Step 1: When the industrial process is running normally, obtain process data at N consecutive sampling times and arrange them into N column vectors x1, x1, …, x2 in the order of sampling time. N , and then according to the formula X i =[x i+d ,x i+d-1 ,…,x i ] Build n dynamic graph matrices X1, X2, ..., X n After that, each dynamic graph matrix is ​​standardized; N represents an integer not less than 500, the dynamic order d is an integer not less than 2 and not greater than 10, n = Nd, number i = 1, 2, ..., n; Step 2: Determine the k nearest neighbors corresponding to each of the normalized dynamic graph matrices, and calculate the average value of these k nearest neighbors, which are recorded as Wherein k is an integer not less than 5 and not greater than 10; Step 3: Perform bidirectional transformation analysis on the dynamic graph to obtain J left transformation vectors u1, u2, ..., u J and J right transformation vectors v1,v2,…,v J Then, through the formula Calculate n score vectors s1,s2,…,s n , and then calculate the covariance matrix Λ=[s1,s2,…,s J ][s1,s2,…,s J ] T / (n-1); where s i represents the i-th score vector, and the parameter J is a smaller than the dynamic graph matrix X i The positive integer of the total number of row vectors in the matrix. The superscript T indicates the transpose of the matrix or vector. Step 4: According to the formula After calculating the abnormal index corresponding to each score vector and determining the threshold of the abnormal index, the following steps 4-1 to 4-3 are repeatedly performed according to the sampling time interval; wherein, ζ i Indicates i Corresponding abnormal indicators; Step 4-1: Obtain the process data of the industrial process at the latest sampling time and form the corresponding column vector x t , and then according to the formula X t =[x t ,x t-1 ,x t-2 ,…,x t-d ] Build the dynamic graph matrix X at the latest sampling time t t After that, X t Perform the same normalization as in step 1; where x t-1 ,x t-2 ,…,x t-d represents the column vector of d consecutive sampling times before t; Step 4-2: Calculate the score vector s in the same way as in step 3 t Then, according to the formula Calculate the abnormal index ζ at the latest sampling time t , and determine whether it is greater than the threshold determined in step 4; if so, execute step 4-3 to decide whether to trigger an abnormal alarm; if not, return to step 4-1; Step 4-3, repeating steps 4-1 to 4-2 at the sampling time interval until abnormal indicators at consecutive f latest sampling moments are obtained, and if they are all greater than the threshold determined in step 4, an abnormal alarm is triggered; wherein f is a positive integer not less than 3; Its characteristics are: the implementation process of the dynamic graph bidirectional transformation analysis includes the following steps 3-1 to 3-4; Step 3-1: After setting the parameter J, set j=1, 2, ..., J in sequence, and execute the following step 3-2 to obtain the corresponding j-th left transformation vector u j and the jth right transformation vector v j ; Step 3-2: Randomly initialize u j , and then repeat the steps 3-3 to 3-4 shown below until u j convergence; Step 3-3: Solve the generalized eigenvalue problem The smallest eigenvalue λ j After the corresponding eigenvector v, Calculate the jth right transformation vector v j ; Step 3-4: If j is equal to 1, solve the generalized eigenvalue problem Gu = λ j The smallest eigenvalue λ in Lu j The corresponding eigenvector u; if j>1, solve the eigenvalue problem The eigenvector u corresponding to the minimum eigenvalue in Update j ;in, I represents the identity matrix, U j-1 =[u1,u2,…,u j-1 ], 2. The anomaly detection method according to claim 1, characterized in that: The step 3-1 of the dynamic graph bidirectional transformation analysis is replaced by: after setting the parameter J, set j=1, 2, ..., M in sequence, and execute the step 3-2 to obtain the corresponding j-th left transformation vector u j , the jth right transformation vector v j , and the minimum eigenvalue λ j ; Wherein, M is equal to the dimension of the column vector in step 1; and then add the following step 3-5 after step 3-4: Step 3-5: Substitute λ1,λ2,Ω,λ M After arranging in ascending order, the left transformation vectors and the right transformation vectors corresponding to the first J smallest eigenvalues ​​are used as the J left transformation vectors and the J right transformation vectors obtained in step 3, respectively.

3. The anomaly detection method according to claim 1 or 2, characterized in that: The step 2 is as follows: Step A-1 to Step A-3 are as follows: i Determine its corresponding k nearest neighbors: Step A-1: ​​According to the formula Calculate X i The distance between the remaining n-1 dynamic graph matrices; where X m Represents the mth dynamic graph matrix, numbered m=1,2,…,n and m≠i, d i,m Represents X i With X m The distance between Indicates calculation of X i -X m The sum of the squares of all elements in ; Step A-2: For any given number p=1, 2, ..., n not equal to i and m, determine d i,p +d p,m The minimum value that can be obtained, and then update d according to the minimum value i,m , the update rule is: if the minimum value is less than d i,m , then update d i,m is equal to the minimum value; if not, keep d i,m unchanged; among them, d i,p Represents X i The distance between the pth dynamic graph matrix, d p,m Represents the pth dynamic graph matrix and X m The distance between Step A-3: According to the updated distance, i The k dynamic graph matrices with the smallest distance between them are determined as X i The corresponding k nearest neighbors.

4. An abnormality detection method according to claim 1 or 2, characterized in that: The step 2 is as follows: Steps B-1 to B-3 are as follows: i Determine its corresponding k nearest neighbors: Step B-1, randomly initialize the feature vector u, and then repeatedly perform steps B-1.1 to B-1.2 as shown below until u converges, according to the formula Calculate X i With X m The difference distance between Among them, X m Represents the mth dynamic graph matrix, numbered m=1,2,…,n and m≠i; Step B-1.

1. Solve the generalized eigenvalue problem After the eigenvector v corresponding to the largest eigenvalue η in Update v; Step B-1.2: Solve the generalized eigenvalue problem After the eigenvector u corresponding to the maximum eigenvalue η in Update u; Step B-2: Repeat step B-1 until you get X i After the difference distance between the other n-1 dynamic graph matrices, i The k dynamic graph matrices with the smallest difference distance are determined as X i The corresponding k nearest neighbors.

5. The anomaly detection method according to claim 1 or 2, characterized in that: The step 2 is as follows: Steps C-1 to C-6 are as follows: i Determine its corresponding k nearest neighbors: Step C-1, according to the formula Calculate the distance between n dynamic graph matrices; where X m Represents the mth dynamic graph matrix, numbered m=1,2,…,n and m≠i, d i,m Represents X i and the mth dynamic graph matrix X m The distance between Step C-2, after setting two width parameters β1 and β2, and the skewness parameter ε, calculate the basis function values ​​between n dynamic graph matrices according to the formula shown below; Among them, f i,m Represents X i With X m The basis function values ​​between β1 and β2 are both greater than zero, ε is not greater than 0.001 and greater than zero, Represents X i The k-order nearest neighbor distance; Step C-3, according to the formula Calculate each element in the transfer matrix P; where P i,m represents the element in the i-th row and m-th column of P, Indicates that X i The numbers of the k dynamic graph matrices with the smallest distance between them; Step C-4: After randomly initializing the feature vector α, continuously update α according to the formula α=αP until it converges, then merge n repeated α to form the feature matrix A, and calculate the basic matrix B=(I n -P+A) -1 ; Step C-5: According to the formula Calculate X i With X m The similarity S between i,m , where C i,m =α m B i,m +α i B m,i -α i α m -α i δ i,m , α i and α m denote the i-th and m-th elements in α, respectively, and B i,m represents the element in the i-th row and m-th column of B. When i=m, δ i,m =1, other casesδ i,m =0; Step C-6, and X i The k dynamic graph matrices with the greatest similarity are determined as X i The corresponding k nearest neighbors.

6. The anomaly detection method according to claim 5, characterized in that: In step C-2, β1 is set equal to the average value of the distances between the n dynamic graph matrices calculated in step C-1.

7. An anomaly detection method according to any one of claims 1, 2, or 6, characterized in that: The calculation process of the average value of the k nearest neighbors in step 2 is as follows: Step D-1: X i The k nearest neighbors of Then, according to the formula Calculate X i and its cth nearest neighbor The weighting coefficient between Among them, the width parameter σ i Equal to X i The average distance between a node and its k nearest neighbors, where the nearest neighbors are numbered c = 1, 2, ..., k; Step D-2: According to the formula calculate 8. The anomaly detection method according to any one of claims 1, 2, or 6, characterized in that: The calculation process of the average value of the k nearest neighbors in step 2 is as follows: Step E-1: X by column i Expand into column vector z i , and then in the same way Expand into corresponding column vectors, recorded as in, Represents X i The k nearest neighbors of ; Step E-2: Combine to form the input matrix Then, through the formula Calculate the weight coefficient vector w i ; Step E-3: According to the formula calculate in, Indicates w i The first to the kth elements in .

9. An abnormality detection device based on dynamic graph bidirectional transformation, comprising a parameter initialization module, a data acquisition module, a data processing module, a display module, and an abnormality alarm; wherein, The parameter initialization module provides the data processing module with the mean and standard deviation required for standardization, as well as the J left-transformed vectors u1, u2, …, u required for calculating the anomaly detection index. J , J right transformation vectors v1,v2,…,v J and covariance matrix Λ, and provide a threshold for judging whether it is abnormal for the display module; The data acquisition module dynamically stores the latest process data of d consecutive sampling moments, and can acquire the latest process data of the sampling moment in real time according to the sampling time interval, and send it together with the dynamically stored process data of d consecutive sampling moments to the data processing module in real time; After the data processing module receives the process data of d+1 sampling moments, it constructs the dynamic graph matrix X t , and then use the mean value and standard deviation provided by the parameter initialization module to standardize it, and then use the J left transformation vectors u1,u2,…,u provided by the parameter initialization module J , J right transformation vectors v1,v2,…,v J and covariance matrix Λ, calculate the corresponding anomaly detection index ζ t , and finally t The specific value is sent to the display module; the display module receives the abnormal detection index ζ sent by the data processing module t After obtaining the specific value of the anomaly detection index, the corresponding scatter points are plotted with the specific value of the anomaly detection index as the ordinate and the latest sampling time as the abscissa. At the same time, the threshold provided by the parameter initialization module is displayed as a horizontal line. Then, it is determined whether the anomaly detection index of the consecutive f latest sampling times is greater than the threshold; if so, the anomaly alarm is turned on; If not, disconnect the abnormal alarm; Its characteristic is that the parameter initialization module will immediately obtain the process data of N consecutive sampling moments of normal operation from the historical database of the industrial process according to the initialization instruction, firstly form n dynamic graph matrices, and calculate the average value and standard deviation of these n dynamic graph matrices; secondly, after standardizing each dynamic graph matrix using the calculated average value and standard deviation, determine its corresponding k nearest neighbors; then, implement the dynamic graph bidirectional transformation analysis according to steps 3-1 to 3-4, and obtain J left transformation vectors u1, u2, ..., u J , J right transformation vectors v1,v2,…,v J , and the covariance matrix Λ; Finally, calculate the abnormal index corresponding to each score vector and determine the threshold of the abnormal index.

10. The abnormality detection device according to claim 9, characterized in that: The parameter initialization module implements the dynamic graph bidirectional transformation analysis described in claim 2.

Images