Electric power industrial control system abnormal behavior detection method, system, device and medium

Abstract

The invention provides an electric power industrial control system abnormal behavior detection method, system, device and medium, and relates to the technical field of electric power system security, and the method comprises the steps: collecting a communication message, a device state and an operation instruction, extracting a message fingerprint vector to construct a feature matrix, calculating a state transition feature to form a feature correlation tensor, and obtaining a feature correlation tensor; and performing time domain fragmentation on the tensor to construct a graph structure, extracting time sequence graph features by using a graph convolutional network, training a twinning neural network to construct a discrimination boundary, and realizing real-time detection of abnormal behaviors. According to the method, the abnormal behaviors in the complex electric power industrial control system can be accurately identified, the detection accuracy is improved, and the false alarm rate is reduced.

Description

Technical Field

[0001] This invention relates to the field of power system safety technology, and in particular to a method, system, equipment and medium for detecting abnormal behavior in power industrial control systems. Background Technology

[0002] Power control systems are the core control systems of critical infrastructure in the power industry, responsible for the automated control and monitoring of power production, transmission, and distribution processes. With the advancement of digital transformation and smart grid construction, power control systems are increasingly interconnected with the internet and enterprise networks. While improving operational efficiency, this also exposes them to growing cybersecurity threats. Attacks on power control systems can lead to equipment failures, system crashes, and even large-scale power outages, severely impacting the national economy and social life. Therefore, real-time and effective detection of abnormal behavior in power control systems is crucial for ensuring the safe and stable operation of the power system.

[0003] Currently, methods for detecting abnormal behavior in power control systems mainly include rule-based detection methods, statistical analysis methods, and traditional machine learning methods. Rule-based detection methods determine whether system behavior is abnormal by pre-setting security policies and rules, but they struggle to cope with unknown attacks and complex intrusion scenarios, and rule maintenance is costly and updates are delayed. Statistical analysis methods build normal behavior models based on historical data and identify anomalies through deviation metrics, but in complex and ever-changing industrial control environments, these methods often have high false alarm rates and lack adaptability to evolving attack methods. While traditional machine learning methods improve detection accuracy to some extent, they cannot effectively integrate multi-dimensional features such as communication messages, device status, and operation commands when faced with multi-source heterogeneous data from power control systems, resulting in limited ability to identify complex attack scenarios.

[0004] Existing technologies lack effective mechanisms for fusing multi-source heterogeneous data, making it difficult to comprehensively utilize multi-dimensional data such as communication messages, equipment status, and operation instructions in power industrial control systems for abnormal behavior detection. Traditional methods are insufficient in modeling the temporal and state transition characteristics of industrial control systems, failing to accurately capture subtle abnormal changes caused by attack behaviors. Existing detection methods generally lack adaptive learning capabilities, making it difficult to cope with continuously changing network traffic patterns and evolving attack techniques in the power industrial control environment, resulting in a significant decline in detection effectiveness over time. Summary of the Invention

[0005] This invention provides a method, system, device, and medium for detecting abnormal behavior in a power industrial control system, which can solve the problems in the prior art.

[0006] A first aspect of this invention provides a method for detecting abnormal behavior in a power industrial control system based on deep learning, comprising: Collect communication messages, equipment status, and operation commands from the power industrial control system, and perform timing alignment and standardization processing; Extract byte sequence features from communication messages, map each byte to a high-dimensional space, compress the message fingerprint vector using an autoencoder, and construct a communication behavior feature matrix based on the message fingerprint vector; The time difference sequence of device state is calculated, state transition features are identified based on recursive decomposition algorithm, and a feature association tensor is constructed by combining the state transition features with the communication behavior feature matrix. The hybrid feature tensor is obtained by tensor decomposition. Temporal segmentation is performed on the hybrid feature tensor, and feature correlation scores are calculated within each time segment. The feature correlation scores are constructed into a graph structure, and temporal graph features are extracted through a graph convolutional network. The Siamese neural network is trained using the temporal graph features of historical anomaly samples as input, and the discrimination boundary is constructed through comparative learning to obtain the anomaly detector; The time-series graph features generated from real-time data are input into the anomaly detector to calculate the anomaly score; When the abnormal score exceeds the adaptive threshold, an alarm message is sent to the power control system and the abnormal behavior characteristics are recorded.

[0007] In one optional embodiment, byte sequence features are extracted from the communication message, each byte is mapped to a high-dimensional space, and a message fingerprint vector is obtained by compression using an autoencoder. A communication behavior feature matrix is ​​then constructed based on the message fingerprint vector, including: The communication message is segmented according to the functional field and the data field, the byte sequence is extracted, the byte sequence is sampled in segments, and the frequency domain features are obtained through Fourier transform; By sequentially traversing each byte in the byte sequence, a byte co-occurrence matrix is ​​constructed to record the occurrence frequency and transition probability of adjacent bytes, thereby obtaining the byte sequence features; An overlapping sliding window is set on the byte sequence, and multiple local byte sequences are obtained by moving the overlapping sliding window. The entropy feature of each local byte sequence is calculated, and the entropy feature, the frequency domain feature and the byte sequence feature are fused to obtain the sequence statistical feature. Correlation analysis is performed on the sequence statistical features, and a subset of features is selected based on mutual information values ​​and mapped to a high-dimensional space. Location information and protocol information are then added to the high-dimensional space. A context window is set on the byte sequence, the association relationship between bytes is calculated based on the context window, and the association relationship is added to the mapping result in the high-dimensional space; Using a pre-trained autoencoder network, the mapping result of the high-dimensional space is compressed to obtain an initial fingerprint vector, and the feature dimension with high variance contribution rate is selected to construct the message fingerprint vector. Multiple message fingerprint vectors are collected based on temporal relationships, and the similarity features of adjacent vectors are calculated to construct a communication behavior feature matrix.

[0008] In one optional embodiment, the time difference sequence of device states is calculated, state transition features are identified based on a recursive decomposition algorithm, and a feature association tensor is constructed by combining the state transition features with the communication behavior feature matrix. The resulting hybrid feature tensor is obtained through tensor decomposition, including: A multidimensional state matrix is ​​constructed, wherein the row vectors of the multidimensional state matrix correspond to the device state parameters and the column vectors correspond to the time series. The time difference sequence is obtained by calculating the difference between adjacent time points in the multidimensional state matrix. The recursive decomposition algorithm is applied to the time difference sequence to obtain the intrinsic mode function set. The Hilbert transform of the intrinsic mode function set is calculated to obtain the instantaneous frequency features. The power spectral density distribution based on the instantaneous frequency features is used to identify state transition points. The state transition points are used as nodes to construct a state transition graph. The weights of the edges between nodes are calculated based on the changes in physical quantities during the state transition process. The degree distribution and centrality features of the state transition graph are extracted to obtain the state transition features. Calculate the adjacency relationship between nodes in the state transition graph and the conditional mutual information between the communication behavior feature matrix, and construct a feature association tensor. The three dimensions of the feature association tensor correspond to the state transition features, communication behavior features, and time series, respectively. The feature association tensor is decomposed into the product of multiple tensor kernels using a tensor network, and a hybrid feature tensor is obtained based on the association between the tensor kernels.

[0009] In one optional embodiment, the feature association tensor is decomposed into a product of multiple tensor kernels using a tensor network, and a hybrid feature tensor is calculated based on the association between the tensor kernels, including: The feature-associated tensor is optimized through tensor linear decomposition to construct an optimization objective function containing tensor dimension constraints. Multiple tensor kernels are obtained through iterative calculation based on the index shrinkage rule. The connection dimension between the tensor kernels is calculated. The connection topology of the tensor kernels is determined according to the connection dimension and the optimization objective function to form a tensor network. The shrinkage complexity matrix is ​​calculated for the tensor kernel in the tensor network. Based on the shrinkage complexity matrix, the optimal shrinkage path is calculated using a dynamic programming algorithm. Based on the optimal shrinkage path, sub-tensor blocks are divided. Parallel shrinkage calculation is performed on the sub-tensor blocks to obtain intermediate tensor results. Construct the density matrix corresponding to the intermediate tensor result, calculate the reduced density matrix of the density matrix, solve the eigenvalue spectrum of the reduced density matrix, and calculate the entanglement entropy between tensor kernels based on the eigenvalue spectrum; The feature weight coefficients are determined based on the entanglement entropy, and the feature weight coefficients are combined with the tensor kernel in a weighted combination operation to generate a hybrid feature tensor. Tensor dimensionality reduction analysis is performed on the hybrid feature tensor to extract the principal feature components. Based on the principal feature components, the information retention rate and feature discriminant are calculated to generate an optimized hybrid feature tensor.

[0010] In one optional embodiment, temporal segmentation is performed on the hybrid feature tensor, feature correlation scores are calculated within each time segment, the feature correlation scores are constructed into a graph structure, and temporal graph features are extracted through a graph convolutional network, including: Temporal segmentation is performed on the hybrid feature tensor. The size of the dynamic time window is determined based on the variance variation trend of the feature vectors in the hybrid feature tensor. The hybrid feature tensor is divided into multiple time segments according to the size of the dynamic time window, and there is a preset overlap interval between adjacent time segments. For each time segment, a set of feature pairs is constructed, the conditional mutual information and nonlinear correlation coefficient between features in the set of feature pairs are calculated, and the results are fused to obtain a feature correlation score matrix; The feature correlation score matrix is ​​constructed into a graph structure. The nodes in the graph structure are pre-screened based on the locality sensitive hashing algorithm. The graph structure is then optimized using the spectral clustering method to obtain an optimized graph structure containing multiple communities. Based on the community division of the optimized graph structure, node information is propagated in the graph structure through a graph convolutional network with multi-head attention mechanism and a message passing function. A corresponding receptive field range is set for each community. Within the receptive field range, the node features within the community are fused, and the node features between communities are fused based on the connection relationship between communities. Combined with the temporal dependency relationship, the temporal graph features are obtained.

[0011] In one optional embodiment, the feature correlation score matrix is ​​constructed into a graph structure, nodes in the graph structure are pre-screened based on the locality-sensitive hashing algorithm, and the graph structure is optimized using the spectral clustering method to obtain an optimized graph structure containing multiple communities, including: Multiple random mapping transformations are performed on the feature vectors of nodes in the graph structure, the mapping results are quantized into binary codes, and the node similarity matrix is ​​obtained by calculating the Hamming distance of the binary codes. Calculate the Laplacian matrix eigenvalues ​​of the node similarity matrix, determine the optimal number of communities based on the difference of the Laplacian matrix eigenvalues, and select the eigenvectors corresponding to the optimal number of communities to construct a low-dimensional node representation; The neighborhood density and relative distance of each node are calculated using the low-dimensional node representation. The community center node is identified through the density clustering effect, and the initial center point of the community is obtained. Other nodes are assigned to the nearest initial community center point, and the structural distance and attribute distance between the nodes and the initial community center point are calculated to obtain the initial community division results. Perform an edge pruning operation on the initial community division result, remove edges whose internal connection weights are lower than the preset weight threshold, determine the remaining connection edges, and merge communities whose number of remaining connection edges is less than the preset connection threshold to obtain an optimized community structure. Based on the connection relationships of the reconstructed graph structure according to the optimized community structure, an optimized graph structure containing multiple communities is obtained.

[0012] In one optional embodiment, a Siamese neural network is trained using the temporal graph features of historical anomaly samples as input, and a discrimination boundary is constructed through contrastive learning to obtain an anomaly detector comprising: The temporal graph features of historical anomaly samples are input into two sub-networks of the Siamese neural network. The two sub-networks have the same network structure and shared weight parameters. Each sub-network maps the input temporal graph features into a contrast feature vector. Construct a contrastive loss function, which is used to maximize the mutual information between the contrastive feature vectors of similar abnormal samples and minimize the mutual information between the contrastive feature vectors of different types of abnormal samples; Based on the contrastive loss function, a Siamese neural network is trained to obtain a feature extractor for extracting discriminative features; The feature extractor is used to extract the temporal spectrum features of the abnormal samples, and a support vector machine classifier is constructed in the feature space. The classification hyperplane of the support vector machine classifier is the anomaly discrimination boundary. The feature extractor and the anomaly detection boundary are integrated into an anomaly detector.

[0013] A second aspect of this invention provides a deep learning-based abnormal behavior detection system for power industrial control systems, comprising: The first unit is used to collect communication messages, equipment status and operation instructions of the power industrial control system, and to perform timing alignment and standardization processing. The second unit is used to extract byte sequence features from communication messages, map each byte to a high-dimensional space, compress it through an autoencoder to obtain a message fingerprint vector, and construct a communication behavior feature matrix based on the message fingerprint vector; The third unit is used to calculate the time difference sequence of device status, identify state transition features based on recursive decomposition algorithm, construct a feature association tensor by combining state transition features with communication behavior feature matrix, and obtain a hybrid feature tensor by tensor decomposition. The fourth unit is used to perform temporal segmentation on the hybrid feature tensor, calculate the feature correlation score in each time segment, construct the feature correlation score into a graph structure, and extract the temporal graph features through a graph convolutional network; The fifth unit is used to train a Siamese neural network with the temporal graph features of historical anomaly samples as input, and to construct a discrimination boundary through contrastive learning to obtain an anomaly detector; The sixth unit is used to input the time-series graph features generated from real-time data into the anomaly detector to calculate the anomaly score; The seventh unit is used to send an alarm message to the power control system and record the abnormal behavior characteristics when the abnormal score exceeds the adaptive threshold.

[0014] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0015] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0016] In this embodiment of the invention, by extracting message fingerprint vectors and constructing feature association tensors, a deep fusion of communication behavior and equipment status in the power industrial control system is achieved. This effectively captures subtle features of abnormal behavior, improving the accuracy and reliability of detection. Based on the contrastive learning mechanism of time-series graph features and Siamese neural networks, this method can adapt to the dynamically changing working modes of the power industrial control system, reducing the false alarm rate. At the same time, it has strong generalization ability against unknown types of attacks, improving the real-time performance and robustness of anomaly detection. By adopting an adaptive threshold mechanism for anomaly judgment, combined with alarm and feature recording functions, a comprehensive security protection system is provided for the power industrial control system, enhancing the system's ability to cope with complex network attacks and ensuring the safe and stable operation of the power industrial control system. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating the abnormal behavior detection method for power industrial control systems based on deep learning, according to an embodiment of the present invention. Figure 2 Optimize the flowchart for tensor networks; Figure 3 This is a graph showing the impact of the number of random mapping transformations on detection accuracy and computation time. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0020] Figure 1 This is a flowchart illustrating the abnormal behavior detection method for power industrial control systems based on deep learning, as described in an embodiment of the present invention. Figure 1 As shown, the method includes: Collect communication messages, equipment status, and operation commands from the power industrial control system, and perform timing alignment and standardization processing; Extract byte sequence features from communication messages, map each byte to a high-dimensional space, compress the message fingerprint vector using an autoencoder, and construct a communication behavior feature matrix based on the message fingerprint vector; The time difference sequence of device state is calculated, state transition features are identified based on recursive decomposition algorithm, and a feature association tensor is constructed by combining the state transition features with the communication behavior feature matrix. The hybrid feature tensor is obtained by tensor decomposition. Temporal segmentation is performed on the hybrid feature tensor, and feature correlation scores are calculated within each time segment. The feature correlation scores are constructed into a graph structure, and temporal graph features are extracted through a graph convolutional network. The Siamese neural network is trained using the temporal graph features of historical anomaly samples as input, and the discrimination boundary is constructed through comparative learning to obtain the anomaly detector; The time-series graph features generated from real-time data are input into the anomaly detector to calculate the anomaly score; When the abnormal score exceeds the adaptive threshold, an alarm message is sent to the power control system and the abnormal behavior characteristics are recorded.

[0021] In one optional implementation, byte sequence features are extracted from the communication message, each byte is mapped to a high-dimensional space, and a message fingerprint vector is obtained by compression using an autoencoder. A communication behavior feature matrix is ​​then constructed based on the message fingerprint vector, including: The communication message is segmented according to the functional field and the data field, the byte sequence is extracted, the byte sequence is sampled in segments, and the frequency domain features are obtained through Fourier transform; By sequentially traversing each byte in the byte sequence, a byte co-occurrence matrix is ​​constructed to record the occurrence frequency and transition probability of adjacent bytes, thereby obtaining the byte sequence features; An overlapping sliding window is set on the byte sequence, and multiple local byte sequences are obtained by moving the overlapping sliding window. The entropy feature of each local byte sequence is calculated, and the entropy feature, the frequency domain feature and the byte sequence feature are fused to obtain the sequence statistical feature. Correlation analysis is performed on the sequence statistical features, and a subset of features is selected based on mutual information values ​​and mapped to a high-dimensional space. Location information and protocol information are then added to the high-dimensional space. A context window is set on the byte sequence, the association relationship between bytes is calculated based on the context window, and the association relationship is added to the mapping result in the high-dimensional space; Using a pre-trained autoencoder network, the mapping result of the high-dimensional space is compressed to obtain an initial fingerprint vector, and the feature dimension with high variance contribution rate is selected to construct the message fingerprint vector. Multiple message fingerprint vectors are collected based on temporal relationships, and the similarity features of adjacent vectors are calculated to construct a communication behavior feature matrix.

[0022] In one specific implementation, during the communication message processing stage, the received communication message is segmented into two parts: a function field and a data field. Taking a certain industrial control protocol as an example, the total message length is 256 bytes, of which the first 24 bytes are the function field, containing information such as the source address, destination address, and function code, and the last 232 bytes are the data field, containing specific control instructions or status information. A complete byte sequence is extracted from the segmented message, for example, the extracted byte sequence is [0x45, 0x00, 0x00, 0x28, 0xab, 0xcd, 0x40, 0x00...].

[0023] The extracted byte sequence is segmented and sampled at a rate of 128 Hz with a sampling window size of 64 bytes and a window overlap rate of 50%. For each byte sequence within a sampling window, a Fast Fourier Transform (FFT) algorithm is applied to convert it to the frequency domain to obtain frequency domain features. For example, after transforming the 64 bytes of data within a certain window, 32 frequency components are obtained, with energy concentrated in the low-frequency range of 2-8 Hz. The spectral feature vector is denoted as [0.82, 0.65, 0.43, 0.37, 0.22…].

[0024] Each byte in the byte sequence is traversed sequentially to construct a byte co-occurrence matrix. This matrix is ​​a 256×256 two-dimensional array, where the row and column indices represent the values ​​of the previous and current byte, respectively, and the matrix elements represent the number of times each byte pair co-occurs in the sequence. By normalizing the matrix elements, the byte transition probabilities can be obtained. For example, the probability of byte 0x45 being followed by byte 0x00 is 0.75, indicating that these two bytes are highly correlated. Based on the co-occurrence matrix, row and column distribution features, diagonal features, and eigenvalue decomposition results are extracted as byte sequence features.

[0025] An overlapping sliding window is set on the byte sequence, with a window size of 16 bytes and a step size of 4 bytes. For each local byte sequence, its entropy feature is calculated. Taking information entropy as an example, the information entropy value of 16 bytes within a certain local window is 3.45, indicating that the byte values ​​within this window are relatively evenly distributed. Simultaneously, conditional entropy, relative entropy, and other multi-dimensional entropy features are calculated to obtain an entropy feature vector [3.45, 2.87, 1.23, ...]. These entropy features are then concatenated and fused with the previously obtained frequency domain features and byte sequence features to form a comprehensive sequence statistical feature vector with a dimension of 128.

[0026] Correlation analysis was performed on the sequence statistical features to calculate the mutual information value between each pair of features. A mutual information threshold of 0.6 was set, and a subset of features with mutual information values ​​greater than the threshold was selected. The selected feature subset was then mapped to a 256-dimensional high-dimensional space using a nonlinear mapping function. The mapping function employed a radial basis function with a kernel width parameter of 0.8. In the high-dimensional space, the relative position information of bytes within the message was added, and position encoding was used to generate a 20-dimensional position vector using sine and cosine functions. Simultaneously, protocol information was added based on the message's protocol type, represented using one-hot encoding, supporting 10 common protocol types and generating a 10-dimensional protocol vector. The final mapping result in the high-dimensional space has a dimension of 286.

[0027] A context window with a size of 5 is set on the byte sequence, considering the context information of the two bytes before and after the current byte. Based on the context window, the association relationships between bytes are calculated using cosine similarity as a metric to generate association feature vectors. For example, for the byte 0xab in the sequence, its association degree with other bytes in the context is [0.92, 0.87, 0.65, 0.42]. These association features are added to the mapping result in a high-dimensional space, expanding the dimension to 350.

[0028] A pre-trained autoencoder network is used to compress high-dimensional features. This autoencoder consists of 5 layers: an input layer with 350 nodes, an encoder with three hidden layers of 256, 128, and 64 nodes, a bottleneck layer of 32 nodes, and a decoder that is a mirror image of the encoder. The autoencoder compresses the 350-dimensional high-dimensional features into a 32-dimensional initial fingerprint vector. Further analysis of the variance contribution rate of each dimension is conducted, and the top 16 dimensions with a cumulative contribution rate exceeding 95% are selected to construct the final message fingerprint vector. For example, the fingerprint vector of a certain message might be [0.72, -0.56, 0.34, -0.28, 0.21...].

[0029] Multiple message fingerprint vectors from a continuous communication process are collected based on temporal relationships, typically 20-50 consecutive messages from each session. Similarity features between adjacent fingerprint vectors are calculated, including Euclidean distance, cosine similarity, and Hamming distance. For example, the Euclidean distance between two adjacent message fingerprint vectors is 0.23, the cosine similarity is 0.92, and the Hamming distance is 4. These similarity features are combined with the message fingerprint vectors to construct a communication behavior feature matrix. The final generated feature matrix has a dimension of N×20, where N is the number of messages, and 20 is the feature dimension of each message (16-dimensional fingerprint vector + 4-dimensional similarity features). This feature matrix can be used for subsequent communication behavior analysis, anomaly detection, and other tasks.

[0030] In one optional implementation, the time difference sequence of device states is calculated, state transition features are identified based on a recursive decomposition algorithm, and a feature association tensor is constructed by combining the state transition features with a communication behavior feature matrix. Tensor decomposition is then used to obtain a hybrid feature tensor, including: A multidimensional state matrix is ​​constructed, wherein the row vectors of the multidimensional state matrix correspond to the device state parameters and the column vectors correspond to the time series. The time difference sequence is obtained by calculating the difference between adjacent time points in the multidimensional state matrix. The recursive decomposition algorithm is applied to the time difference sequence to obtain the intrinsic mode function set. The Hilbert transform of the intrinsic mode function set is calculated to obtain the instantaneous frequency features. The power spectral density distribution based on the instantaneous frequency features is used to identify state transition points. The state transition points are used as nodes to construct a state transition graph. The weights of the edges between nodes are calculated based on the changes in physical quantities during the state transition process. The degree distribution and centrality features of the state transition graph are extracted to obtain the state transition features. Calculate the adjacency relationship between nodes in the state transition graph and the conditional mutual information between the communication behavior feature matrix, and construct a feature association tensor. The three dimensions of the feature association tensor correspond to the state transition features, communication behavior features, and time series, respectively. The feature association tensor is decomposed into the product of multiple tensor kernels using a tensor network, and a hybrid feature tensor is obtained based on the association between the tensor kernels.

[0031] In one specific implementation, power control systems generate a large amount of equipment status data and communication behavior data. For the equipment status data, constructing a multidimensional state matrix is ​​fundamental to analyzing system behavior. The row vectors of this multidimensional state matrix correspond to various equipment status parameters, such as transformer temperature, voltage, current, and power factor; the column vectors correspond to time series, i.e., data points collected in chronological order. For example, for a system containing 50 status parameters, if data is collected every minute and monitored continuously for 24 hours, a 50×1440 state matrix is ​​formed.

[0032] The time difference sequence is obtained by calculating the difference between adjacent time points in the multidimensional state matrix. Specifically, for each row of data in the state matrix, the difference between two adjacent time points is calculated to form a difference sequence. This step helps to highlight state changes and eliminate static baseline interference. Taking transformer temperature as an example, if the temperature sequence is [45.2, 45.3, 45.8, 46.5, 47.2]℃, the corresponding difference sequence is [0.1, 0.5, 0.7, 0.7]℃.

[0033] A recursive decomposition algorithm is applied to a time-difference sequence to obtain a set of intrinsic mode functions (EMFs). During the recursive decomposition, local extrema of the sequence are identified. Upper and lower envelopes are constructed using these maxima and minima, and the mean curve is calculated. The mean curve is subtracted from the original signal to obtain candidate EMFs. This process is repeated until a candidate EMF satisfies the definition of an EMF. The EMF is then subtracted from the original signal. The above steps are repeated for the residual signal until it becomes a monotonic function or a constant. Taking voltage fluctuations as an example, the decomposition may yield 5-10 EMFs reflecting different frequency characteristics.

[0034] The instantaneous frequency characteristics are obtained by calculating the Hilbert transform of the set of eigenmode functions. The Hilbert transform is then applied to each eigenmode function to obtain an analytic signal, from which the instantaneous amplitude and instantaneous phase are calculated. The instantaneous frequency is equal to the time derivative of the instantaneous phase. For example, calculating the eigenmode functions of a current fluctuation yields characteristics reflecting different frequency components within the range of 0.1Hz to 50Hz.

[0035] State transition points are identified based on the power spectral density distribution of instantaneous frequency characteristics. The power spectral density is calculated within a time window at the instantaneous frequency, and its distribution characteristics are analyzed. When the system state changes, the power spectral density distribution exhibits significant changes. A threshold is set to judge the magnitude of the change; time points exceeding the threshold are marked as state transition points. In practical applications, a sudden change in power spectral density exceeding three times the average value can be identified as a state transition point.

[0036] A state transition graph is constructed using state transition points as nodes. Each state transition point represents the moment when the system transitions from one stable state to another. Directed edges are formed between adjacent state transition points, representing state transition paths. For example, during a 24-hour monitoring period, 15-20 state transition points may be identified, forming the nodes of the state transition graph.

[0037] The weights of edges between nodes are calculated based on the changes in physical quantities during state transitions. These changes can include the magnitude of changes in parameters such as voltage, current, and power before and after the state transition. Weights can be calculated using methods such as Euclidean distance, Mahalanobis distance, or correlation coefficients. For example, if the voltage changes by 0.5kV, the current changes by 2A, and the power changes by 0.8kW between two states, the normalized combined distance can be used as the edge weight.

[0038] State transition features are obtained by extracting the degree distribution and centrality features from the state transition graph. The degree distribution reflects the statistical count of the number of edges connected to each node, while the centrality features include degree centrality, proximity centrality, and betweenness centrality. Degree centrality represents the number of directly connected neighbors of a node, proximity centrality measures the average distance from a node to all other nodes, and betweenness centrality represents the frequency with which a node lies on the shortest path between other node pairs. These features reflect the patterns and critical states of state transitions. In practice, anomalous behavior often manifests as abrupt changes in centrality features.

[0039] The conditional mutual information between the adjacency relationships of nodes in the state transition graph and the communication behavior feature matrix is ​​calculated to construct a feature correlation tensor. The communication behavior feature matrix includes indicators such as communication traffic, protocol type, and communication frequency. The conditional mutual information measures the correlation between state transition features and communication behavior features under specific time conditions. The constructed feature correlation tensor has three dimensions corresponding to state transition features, communication behavior features, and time series, respectively. For example, for 10 state transition features, 15 communication behavior features, and 1440 time points, a 10×15×1440 tensor is formed.

[0040] Tensor networks can be used to decompose feature association tensors into the product of multiple tensor kernels. Tensor network decomposition can employ methods such as Tucker decomposition, CP decomposition, or tensor column decomposition. Taking Tucker decomposition as an example, the original tensor is decomposed into the product of a kernel tensor and three factor matrices. In practical applications, a 10×15×1440 tensor can be decomposed into a 5×7×50 kernel tensor and its corresponding three factor matrices.

[0041] A hybrid feature tensor is obtained based on the correlation calculation between tensor kernels. The element values ​​and distribution characteristics of the tensor kernels are analyzed to identify the feature combinations corresponding to high-value elements in the kernel tensor, and the correlation strength of these feature combinations is calculated. Based on these feature combinations, a hybrid feature tensor is constructed, which comprehensively reflects the coupling relationship between state transition characteristics and communication behavior characteristics. Practice shows that abnormal behavior often manifests as abnormal coupling between specific state transition characteristics and communication behavior characteristics. By analyzing the variation patterns of the hybrid feature tensor, abnormal behavior in power control systems can be effectively detected.

[0042] In one optional implementation, the feature association tensor is decomposed into a product of multiple tensor kernels using a tensor network, and a hybrid feature tensor is calculated based on the association between the tensor kernels, including: The feature-associated tensor is optimized through tensor linear decomposition to construct an optimization objective function containing tensor dimension constraints. Multiple tensor kernels are obtained through iterative calculation based on the index shrinkage rule. The connection dimension between the tensor kernels is calculated. The connection topology of the tensor kernels is determined according to the connection dimension and the optimization objective function to form a tensor network. The shrinkage complexity matrix is ​​calculated for the tensor kernel in the tensor network. Based on the shrinkage complexity matrix, the optimal shrinkage path is calculated using a dynamic programming algorithm. Based on the optimal shrinkage path, sub-tensor blocks are divided. Parallel shrinkage calculation is performed on the sub-tensor blocks to obtain intermediate tensor results. Construct the density matrix corresponding to the intermediate tensor result, calculate the reduced density matrix of the density matrix, solve the eigenvalue spectrum of the reduced density matrix, and calculate the entanglement entropy between tensor kernels based on the eigenvalue spectrum; The feature weight coefficients are determined based on the entanglement entropy, and the feature weight coefficients are combined with the tensor kernel in a weighted combination operation to generate a hybrid feature tensor. Tensor dimensionality reduction analysis is performed on the hybrid feature tensor to extract the principal feature components. Based on the principal feature components, the information retention rate and feature discriminant are calculated to generate an optimized hybrid feature tensor.

[0043] In one specific implementation, raw feature data is acquired, which can come from multimodal data sources such as images, text, and audio. Preprocessing operations are performed on the feature data for different modalities, including data standardization, missing value imputation, and outlier handling. For example, image features can be standardized with a mean of 0 and a standard deviation of 1; text features can be weighted using TF-IDF. The preprocessed features are organized into a multidimensional tensor structure, forming a feature association tensor. Specifically, for a dataset containing 100 samples, each with 20 image features and 30 text features, a third-order feature association tensor with a shape of 100×20×30 can be constructed.

[0044] To effectively decompose the feature association tensor, a tensor network is used to decompose it into a product of multiple tensor kernels. Specifically, tensor linear decomposition optimization is applied to the feature association tensor, constructing an optimization objective function that includes tensor dimensionality constraints. This objective function aims to minimize the Frobenius norm difference between the original and reconstructed tensors, while considering the dimensionality constraints of the tensor kernels. Multiple tensor kernels are iteratively computed using alternating least squares. For example, a third-order tensor is decomposed into four connected low-order tensor kernels, with the rank of each kernel controlled to be within 30% of the original dimension.

[0045] After the tensor kernels are computed, the connection dimensions between each kernel are determined. Based on the rules of virtual and physical edges, the connection topology of the tensor kernels is determined. For example, for four tensor kernels A, B, C, and D, a star topology may be formed where A connects to B, B connects to C and D, and the connection dimensions are 10, 8, and 12, respectively.

[0046] To efficiently compute the shrinkage result of a tensor network, a shrinkage complexity matrix is ​​calculated for the tensor kernels in the tensor network. This matrix records the computational complexity required to shrink any two tensor kernels. For example, in a 4×4 shrinkage complexity matrix, each element represents the number of floating-point operations required to shrink the corresponding two tensor kernels. Based on this matrix, a dynamic programming algorithm is used to calculate the optimal shrinkage path that minimizes the overall computational complexity. In practical implementation, for a network containing 6 tensor kernels, the possible optimal shrinkage path is to first shrink kernels 1 and 2, then shrink the result and kernel 3, and so on, reducing the total computational complexity from the original 10. 12 Reduced to 10 8 Floating-point arithmetic.

[0047] Based on the optimal shrinkage path, sub-tensor blocks are divided for parallel computation. For example, if a shrinkage task is divided into 3 sub-blocks, each responsible for computation within a specific dimension range, parallel shrinkage is performed through multi-threaded or distributed computation to obtain intermediate tensor results.

[0048] To assess the correlation between tensor kernels, a density matrix corresponding to the intermediate tensor results is constructed. Taking a fourth-order intermediate tensor as an example, it is reshaped into a two-dimensional matrix, and its covariance is calculated as the density matrix. Subsequently, the reduced density matrix is ​​calculated, specifically by performing partial trace operations on some dimensions. For a density matrix of size 16×16, an 8×8 reduced density matrix may be obtained. The eigenvalue spectrum of this reduced density matrix is ​​solved, yielding an eigenvalue distribution such as [0.5, 0.25, 0.125, 0.0625, 0.0312, 0.0156, 0.0078, 0.0039]. Based on the eigenvalue spectrum, the entanglement entropy between tensor kernels is calculated; a higher entanglement entropy indicates a more complex correlation between features.

[0049] Based on the entanglement entropy value, feature weight coefficients are assigned to different tensor kernels. For example, a kernel with an entanglement entropy of 1.2 might receive a weight of 0.8, while a kernel with an entanglement entropy of 0.5 might receive a weight of 0.4. The feature weight coefficients are then weighted and combined with the tensor kernels to generate a hybrid feature tensor. Specifically, each tensor kernel is multiplied by its corresponding weight, and the kernels are merged through a shrinking operation to generate the hybrid feature tensor.

[0050] To optimize the hybrid feature tensor, tensor dimensionality reduction analysis is performed. High-order singular value decomposition (SVD) is used to extract principal feature components, for example, reducing the dimensionality from the original 50-dimensional feature space to a 15-dimensional principal feature space. The system calculates the information retention rate based on the principal feature components; a retention rate of 95% or higher indicates that the reduced features retain the main information of the original. Simultaneously, feature discriminancy is calculated to evaluate the distinguishing ability of the features in classification or clustering tasks. Combining the information retention rate and feature discriminancy, the system generates an optimized hybrid feature tensor.

[0051] Existing technologies mainly employ traditional tensor decomposition methods such as Tucker decomposition and CP decomposition for feature extraction and dimensionality reduction. However, these methods suffer from high computational complexity and insufficient feature correlation capture when dealing with high-dimensional, multimodal feature data. In particular, when processing heterogeneous data sources, traditional methods struggle to effectively represent the complex correlation structures between features of different modalities.

[0052] The method in this embodiment addresses the efficiency and accuracy issues of existing feature fusion techniques when processing high-dimensional multimodal data, particularly the insufficient expressive power for complex correlation structures between heterogeneous data sources. By introducing the concepts of tensor networks and entanglement entropy from quantum information theory, a theoretical model of the correlation between features is established, achieving more accurate feature extraction and fusion.

[0053] like Figure 2The diagram illustrates the temporal interaction between components in a tensor network optimization system, depicting the complete workflow from data input to final output. Five main participants are shown: the data scientist, the tensor optimization engine, the shrinking path optimizer, the entanglement entropy analyzer, and the feature dimensionality reduction analyzer. They are represented by lifelines and message passing to show the chronological order of system operation. The interaction begins with the data scientist submitting a feature association tensor to the tensor optimization engine, which processes and constructs the tensor network, returning a confirmation message upon completion. The tensor optimization engine then passes the constructed tensor network to the shrinking path optimizer for calculating the shrinking path. After calculation, the shrinking path optimizer returns the optimal shrinking path result and passes the intermediate tensor results to the entanglement entropy analyzer. Upon receiving the intermediate tensor results, the entanglement entropy analyzer performs density matrix construction and entanglement entropy calculation, returning the entanglement entropy calculation result and providing the feature weight coefficients to the feature dimensionality reduction analyzer. The feature dimensionality reduction analyzer uses these weight coefficients to perform tensor dimensionality reduction analysis, extracting the principal feature components and ultimately generating an optimized hybrid feature tensor. The result is then returned to the data scientist, completing the entire tensor network optimization process.

[0054] In one optional implementation, temporal segmentation is performed on the hybrid feature tensor, feature correlation scores are calculated within each time segment, and the feature correlation scores are constructed into a graph structure. Temporal graph features are extracted through a graph convolutional network, including: Temporal segmentation is performed on the hybrid feature tensor. The size of the dynamic time window is determined based on the variance variation trend of the feature vectors in the hybrid feature tensor. The hybrid feature tensor is divided into multiple time segments according to the size of the dynamic time window, and there is a preset overlap interval between adjacent time segments. For each time segment, a set of feature pairs is constructed, the conditional mutual information and nonlinear correlation coefficient between features in the set of feature pairs are calculated, and the results are fused to obtain a feature correlation score matrix; The feature correlation score matrix is ​​constructed into a graph structure. The nodes in the graph structure are pre-screened based on the locality sensitive hashing algorithm. The graph structure is then optimized using the spectral clustering method to obtain an optimized graph structure containing multiple communities. Based on the community division of the optimized graph structure, node information is propagated in the graph structure through a graph convolutional network with multi-head attention mechanism and a message passing function. A corresponding receptive field range is set for each community. Within the receptive field range, the node features within the community are fused, and the node features between communities are fused based on the connection relationship between communities. Combined with the temporal dependency relationship, the temporal graph features are obtained.

[0055] In one specific implementation, when performing time-domain slicing on the hybrid feature tensor, raw data from the power industrial control system is acquired, including heterogeneous data from multiple sources such as equipment operating status, network traffic, and operation logs. This data is preprocessed, including data cleaning, normalization, and feature extraction, to obtain the hybrid feature tensor. The hybrid feature tensor has dimensions of N×M×T, where N represents the number of samples, M represents the feature dimension, and T represents the time step.

[0056] The size of the dynamic time window is determined based on the variance variation trend of the feature vectors in the mixed feature tensor. Specifically, the variance sequence of the features in the time dimension, Var = {var1, var2, ..., var...}, is calculated. t}, where var t Let represent the variance of the feature at time t. The rate of change of variance is calculated using the sliding window method. When the rate of change of variance exceeds a preset threshold (e.g., 0.15), the feature is considered to have changed significantly. The frequency of significant feature changes under different window sizes is statistically analyzed, and the smallest window size that can capture more than 80% of the feature changes is selected as the dynamic time window size W. For example, in the monitoring data of a power substation, analysis shows that when the window size is 12 hours, approximately 85% of the significant feature changes can be captured; therefore, the dynamic time window size is set to 12 hours.

[0057] The hybrid feature tensor is divided into multiple time segments based on a dynamic time window size W, with a preset overlap interval between adjacent time segments. The size of the overlap interval is set to 25% of the time window size, i.e., 0.25W. For example, when W is 12 hours, the overlap interval is 3 hours. This division ensures a smooth transition between time segments and avoids information loss.

[0058] For each time segment, a feature pair set is constructed. The feature pair set contains all possible combinations of feature pairs within that time segment, totaling M×(M-1) / 2 feature pairs. For a certain power industrial control system, if the feature dimension M is 64, then the number of feature pairs is 2016.

[0059] The conditional mutual information between features in the feature pair set is calculated. For a feature pair (X, Y), the conditional mutual information represents the degree of dependence between X and Y given other features Z. The conditional mutual information is calculated using a nonparametric method based on kernel density estimation. For continuous features, a Gaussian kernel function is used; for discrete features, a histogram is used to estimate the probability distribution. To improve computational efficiency, an approximate algorithm is used; when the sample size exceeds 10,000, 5,000 samples are randomly sampled for calculation.

[0060] The nonlinear correlation coefficient between features in the feature pair set is calculated. This nonlinear correlation coefficient is calculated using the distance correlation method, which detects nonlinear dependencies between features. For a feature pair (X, Y), the distance matrix between the features is calculated, the distance matrix is ​​centered, and the distance correlation coefficient is calculated. In practical applications, when dealing with high-dimensional features, dimensionality reduction techniques are used to reduce computational complexity, such as reducing the original 64-dimensional features to 32 dimensions using principal component analysis before calculating the nonlinear correlation coefficient.

[0061] The feature correlation score matrix is ​​obtained by fusing conditional mutual information and nonlinear correlation coefficient. A weighted average is used for fusion, with the weights dynamically adjusted based on the historical performance of the two indicators. For example, in an anomaly detection task of a power control system, the empirical weights are set to 0.6 for conditional mutual information and 0.4 for the nonlinear correlation coefficient. The fusion yields an M×M dimensional feature correlation score matrix S.

[0062] The feature relevance score matrix is ​​constructed as a graph structure G=(V,E,W), where V represents the set of nodes (corresponding to features), E represents the set of edges (corresponding to feature pairs), and W represents the edge weights (corresponding to relevance scores). To reduce the influence of noise, a threshold θ (e.g., 0.3) is set; when the relevance score is less than θ, the corresponding edge is pruned.

[0063] Locality-Sensitive Hash (LSH) is used to pre-filter nodes in the graph structure. Specifically, the MinHash algorithm is employed, generating k hash signatures (k = 10) for each node. Nodes with similar hash signatures are grouped together. This step initially clusters the original 64 nodes into approximately 15 groups, significantly reducing the computational load in subsequent processing.

[0064] Spectral clustering is used to optimize the graph structure, resulting in an optimized graph structure containing multiple communities. The Laplacian matrix of the graph is calculated using spectral clustering, and the first c feature vectors are extracted (c is the preset number of communities, typically set to 4-8 depending on business characteristics). K-means clustering is then applied in the feature space to complete the clustering. In a power substation monitoring system, 64 feature nodes are divided into 6 communities, corresponding to functional groups such as voltage monitoring, current monitoring, equipment status, network communication, operation control, and environmental parameters.

[0065] Based on the community partitioning of the optimized graph structure, node information is propagated through a graph convolutional network with a multi-head attention mechanism. The multi-head attention mechanism contains h attention heads (h is set to 8), each independently learning the relationships between nodes and generating an attention score. The graph convolutional operation updates the representation of the central node by aggregating information from neighboring nodes. For node i, its feature update formula includes a weighted sum of its own features and neighboring features, with the weights determined by the attention score.

[0066] Message passing functions propagate node information within a graph structure. The message passing process consists of three steps: message generation, message aggregation, and node update. In the message generation phase, each node sends transformed feature information to its neighbors; in the message aggregation phase, each node collects messages from its neighbors and performs a weighted sum; in the node update phase, the node representation is updated by combining the original features and the aggregated messages. Message passing can be performed in multiple rounds (e.g., 3 rounds) to capture a wider range of dependencies.

[0067] A corresponding sensing field is set for each community. The sensing field is defined as the farthest node distance that a node can sense. Different sensing field ranges are set for different communities based on their functional characteristics. For example, for the voltage monitoring community, the sensing field range is set to 2, which can sense all nodes within a 2-hop distance; for the network communication community, the sensing field range is set to 3 to capture a wider range of network dependencies.

[0068] The features of nodes within the receptive field are fused. The fusion method employs attention-weighted summation, with the attention weight dynamically calculated based on the node's centrality within the community and the feature importance of the current time segment. For example, in the current monitoring community, the centrality of the main loop current feature is high, and its attention weight is approximately 0.25, while the attention weight of the auxiliary loop current feature is approximately 0.1.

[0069] The feature fusion between communities is based on the connectivity relationships between them. These connections are defined by the edge density of nodes between communities; higher edge density indicates stronger inter-community relationships. Feature fusion employs a gating mechanism, learning a gating parameter to control the degree of information flow. For example, in an anomaly detection task, a gating parameter of 0.7 between the voltage monitoring community and the equipment status community indicates strong information exchange between the two communities.

[0070] Temporal graph features are obtained by combining temporal dependencies. The graph features of different time segments are arranged chronologically, and temporal dependencies are captured through a bidirectional long short-term memory (LSTM) network. Specifically, the hidden layer dimension of the LSTM network is set to 128, and the time step is equal to the number of time segments. After temporal modeling, temporal graph features containing spatiotemporal information are obtained. These features have a dimension of N×F, where N is the number of samples and F is the feature dimension (set to 256).

[0071] Based on the extracted temporal graph features, classifiers (such as random forests or deep neural networks) are used to detect abnormal behavior.

[0072] In one optional implementation, the feature correlation score matrix is ​​constructed into a graph structure, nodes in the graph structure are pre-screened based on the locality-sensitive hashing algorithm, and the graph structure is optimized using spectral clustering to obtain an optimized graph structure containing multiple communities, including: Multiple random mapping transformations are performed on the feature vectors of nodes in the graph structure, the mapping results are quantized into binary codes, and the node similarity matrix is ​​obtained by calculating the Hamming distance of the binary codes. Calculate the Laplacian matrix eigenvalues ​​of the node similarity matrix, determine the optimal number of communities based on the difference of the Laplacian matrix eigenvalues, and select the eigenvectors corresponding to the optimal number of communities to construct a low-dimensional node representation; The neighborhood density and relative distance of each node are calculated using the low-dimensional node representation. The community center node is identified through the density clustering effect, and the initial center point of the community is obtained. Other nodes are assigned to the nearest initial community center point, and the structural distance and attribute distance between the nodes and the initial community center point are calculated to obtain the initial community division results. Perform an edge pruning operation on the initial community division result, remove edges whose internal connection weights are lower than the preset weight threshold, determine the remaining connection edges, and merge communities whose number of remaining connection edges is less than the preset connection threshold to obtain an optimized community structure. Based on the connection relationships of the reconstructed graph structure according to the optimized community structure, an optimized graph structure containing multiple communities is obtained.

[0073] In one specific implementation, the feature correlation score matrix is ​​constructed as a graph structure, where each node in the graph structure represents a feature, and the edge weights between nodes represent the correlation scores between features. The graph structure can be represented as a weighted undirected graph G=(V,E,W), where V represents the set of nodes, E represents the set of edges, and W represents the edge weight matrix.

[0074] When pre-screening nodes in a graph structure, the Locality Sensitive Hash (LSH) algorithm is used for initial clustering. Specifically, multiple random mapping transformations are performed on the feature vector of each node in the graph structure. In practice, 10-20 random mapping transformations can be performed, with each transformation using a different random projection vector. Assuming the dimension of the node feature vector is d, k random projection vectors r_i (i=1, 2, ..., k) can be generated, each with a dimension of d. For each node's feature vector v, its dot product with the random projection vector is calculated. If the dot product is greater than 0, it is encoded as 1; otherwise, it is encoded as 0, resulting in a binary code of length k. For example, for a node's feature vector [0.5, 0.3, -0.2, 0.8], using four random projection vectors yields dot product values ​​of 0.6, -0.3, 0.4, and -0.1, respectively, then the binary code for this node is 1010.

[0075] A node similarity matrix S is constructed by calculating the Hamming distance between the binary codes of different nodes. The Hamming distance refers to the number of different characters at corresponding positions in two strings of equal length. For example, if node A's code is 1010 and node B's code is 1110, their Hamming distance is 1. Similarity can be defined as 1 minus the normalized Hamming distance, i.e., similarity = 1 - (Hamming distance / code length). In the example above, the similarity between nodes A and B is 1 - (1 / 4) = 0.75.

[0076] After constructing the node similarity matrix, calculate its Laplacian matrix. The Laplacian matrix L can be obtained by subtracting the similarity matrix S from the degree matrix D, i.e., L = DS, where D is a diagonal matrix and the diagonal elements are the degrees of the corresponding nodes. Calculate the eigenvalues ​​of the Laplacian matrix and sort them in ascending order. Observe the difference between adjacent eigenvalues. When there is a significant jump in the difference between eigenvalues, the number of eigenvalues ​​corresponding to that position is the optimal number of communities c. For example, assuming the sorted eigenvalues ​​are [0.01, 0.03, 0.08, 0.42, 0.67, 0.89], the difference between adjacent eigenvalues ​​is [0.02, 0.05, 0.34, 0.25, 0.22]. The jump from 0.08 to 0.42 is the most significant, with a difference of 0.34. In this case, the optimal number of communities is determined to be 3.

[0077] We select the eigenvectors corresponding to the first c eigenvalues ​​of the Laplacian matrix to construct a low-dimensional node representation matrix. Each node is represented as a c-dimensional vector in this low-dimensional space. Based on these low-dimensional representations, we calculate the neighborhood density and relative distance of each node. The neighborhood density can be calculated by counting the number of neighbors of a node within a given radius, and the relative distance is the minimum distance from a node to a node with a higher density. Specifically, we choose a suitable cutoff distance dc, and for each node i, we calculate its neighborhood density ρi, which is the number of nodes whose distance to node i is less than dc. Then, for each node i, we calculate its relative distance δi, which is the distance to the nearest node i among all nodes with a higher density.

[0078] Community center nodes are identified by analyzing the distribution of neighborhood density and relative distance among nodes. Community center nodes typically exhibit high neighborhood density and high relative distance. A density-distance decision graph can be plotted to select nodes with both high ρ and high δ values ​​as initial community center points. In practical applications, thresholds can be set, such as selecting nodes where ρ > ρ threshold and δ > δ threshold as community centers. For example, if the ρ threshold is set to 10 and the δ threshold to 0.5, nodes satisfying ρ > 10 and δ > 0.5 will be identified as community centers.

[0079] After determining the initial community center, other nodes are assigned to the nearest community center. Distance calculation considers both structural distance and attribute distance. Structural distance is measured by the shortest path length between nodes in the graph, while attribute distance is based on the Euclidean distance of the node's feature vectors. The combined distance can be expressed as a weighted sum of structural and attribute distances, and the weights can be adjusted according to the specific application scenario. For example, the structural distance weight can be set to 0.6, and the attribute distance weight to 0.4.

[0080] After obtaining the initial community partitioning results, an edge pruning operation is performed to optimize the community structure. A weight threshold, w, is set as the lower limit, and edges with connection weights within a community that are lower than this threshold are removed. For example, if w is set to 0.3, all edges with weights less than 0.3 will be removed. Simultaneously, a connection threshold, c, is set. If the number of remaining connection edges in a community is less than c, that community is merged with its closest community. The similarity between communities can be measured by the average connection strength between nodes in two communities.

[0081] Through the edge pruning and community merging operations described above, an optimized community structure is obtained. Based on this optimized community structure, the connection relationships of the graph structure are reconstructed, resulting in an optimized graph structure containing multiple communities. In the optimized graph structure, nodes within the same community are more tightly connected, while connections between different communities are relatively sparse, thus achieving effective optimization of the graph structure.

[0082] like Figure 3 As shown in the figure, the blue curve represents the NMI score (left Y-axis), the red curve represents the computation time (right Y-axis), and the X-axis represents the number of random mapping transformations. It is clear from the figure that as the number of random mapping transformations increases, the NMI score exhibits a trend of "rapid increase - gradual stabilization," while the computation time increases approximately linearly. Specifically, when the number of random mapping transformations is 5, the NMI score is only 0.62, and the computation time is 12.3 seconds; when the number increases to 20, the NMI score significantly improves to 0.79, and the computation time increases to 22.4 seconds; when the number reaches 50, the NMI score reaches 0.87, and the computation time is 43.5 seconds; thereafter, the improvement in the NMI score becomes slow, and when the number reaches 100, the NMI score only slightly improves to 0.89, but the computation time increases to 85.2 seconds. This result indicates that there is an optimal range for the number of random mapping transformations, approximately between 40 and 50, at which point high detection accuracy can be maintained while avoiding excessive increases in computation time. This discovery has important guiding significance for parameter optimization in practical applications, especially when dealing with large-scale networks, where choosing an appropriate number of random mapping transformations can achieve a good balance between accuracy and efficiency.

[0083] In one optional implementation, a Siamese neural network is trained using the temporal graph features of historical anomaly samples as input, and a discrimination boundary is constructed through contrastive learning to obtain an anomaly detector comprising: The temporal graph features of historical anomaly samples are input into two sub-networks of the Siamese neural network. The two sub-networks have the same network structure and shared weight parameters. Each sub-network maps the input temporal graph features into a contrast feature vector. Construct a contrastive loss function, which is used to maximize the mutual information between the contrastive feature vectors of similar abnormal samples and minimize the mutual information between the contrastive feature vectors of different types of abnormal samples; Based on the contrastive loss function, a Siamese neural network is trained to obtain a feature extractor for extracting discriminative features; The feature extractor is used to extract the temporal spectrum features of the abnormal samples, and a support vector machine classifier is constructed in the feature space. The classification hyperplane of the support vector machine classifier is the anomaly discrimination boundary. The feature extractor and the anomaly detection boundary are integrated into an anomaly detector.

[0084] In one specific implementation, the time-series spectral feature extraction process includes segmenting the original time-series data using a sliding window, and extracting statistical features, frequency domain features, and time-frequency features for each segment. Specifically, statistical features include mean, standard deviation, kurtosis, skewness, etc.; frequency domain features are obtained through Fast Fourier Transform to extract the spectral energy distribution; and time-frequency features are obtained using wavelet transform to acquire the signal's characteristics at different time scales. These features are combined to form a multidimensional vector representing the behavioral characteristics of the time-series data in different dimensions. For example, for the vibration data of a certain device, 96-dimensional time-series spectral features can be extracted, including 32-dimensional statistical features, 32-dimensional frequency domain features, and 32-dimensional time-frequency features.

[0085] A Siamese neural network consists of two identical sub-networks that share the same weight parameters. Each sub-network contains multiple convolutional layers, pooling layers, and fully connected layers. In this implementation, the sub-network consists of three convolutional blocks, each containing a convolutional layer, a batch normalization layer, and a ReLU activation function, followed by a max-pooling layer. The number of filters in the convolutional layers is set to 32, 64, and 128, respectively, and the kernel size is 3×3. After the convolutional blocks, two fully connected layers are connected, with 256 and 128 nodes, respectively, ultimately outputting a 128-dimensional feature vector, i.e., the contrastive feature vector.

[0086] In the comparative learning training process, it is necessary to construct positive sample pairs and negative sample pairs. A positive sample pair refers to two different samples from the same anomaly type, while a negative sample pair refers to two samples from different anomaly types. For example, if there are three anomaly types: bearing failure, gear failure, and motor failure, then two bearing failure samples constitute a positive sample pair, and one bearing failure sample and one gear failure sample constitute a negative sample pair.

[0087] During training, a batch of samples is randomly selected from historical outlier data, with each batch containing 64 samples. For each sample, another sample from the same class is randomly selected to form a positive sample pair, and a sample from a different class is randomly selected to form a negative sample pair. These sample pairs are input into the Siamese neural network, and the loss is calculated using a contrastive loss function. The goal of the contrastive loss function is to minimize the distance between positive sample pairs in the feature space and maximize the distance between negative sample pairs. In the specific implementation, the InfoNCE loss function with a temperature parameter of 0.07 is used; the higher the similarity between two feature vectors, the smaller the loss value; the lower the similarity between two feature vectors, the larger the loss value.

[0088] The Adam optimizer was used during training, with an initial learning rate of 0.001. The learning rate was decayed to 0.1 times its original value every 30 epochs, for a total of 100 epochs. To prevent overfitting, a dropout ratio of 0.5 was added to the fully connected layers.

[0089] After training, a subnetwork of the Siamese neural network is used as a feature extractor. This feature extractor maps the input temporal graph features to a more discriminative feature space. In this feature space, features of similar anomalous samples cluster together, while features of different anomalous samples are more widely distributed.

[0090] Next, a Support Vector Machine (SVM) classifier is constructed in this feature space. The SVM uses a Radial Basis Function (RBF) as its kernel function, and the optimal parameters are determined through grid search, with the penalty parameter C set to 10 and the kernel parameter gamma set to 0.01. The decision boundary of the SVM is the anomaly detection boundary.

[0091] The final anomaly detector consists of a feature extractor and anomaly detection boundaries. In practical applications, when new data arrives, its time-series spectral features are first extracted, then mapped to the feature space by the feature extractor, and finally, a support vector machine classifier is used to determine whether it is an anomaly and what type of anomaly it belongs to.

[0092] The abnormal behavior detection system for power industrial control systems based on deep learning, as described in this embodiment of the invention, includes: The first unit is used to collect communication messages, equipment status and operation instructions of the power industrial control system, and to perform timing alignment and standardization processing. The second unit is used to extract byte sequence features from communication messages, map each byte to a high-dimensional space, compress it through an autoencoder to obtain a message fingerprint vector, and construct a communication behavior feature matrix based on the message fingerprint vector; The third unit is used to calculate the time difference sequence of device status, identify state transition features based on recursive decomposition algorithm, construct a feature association tensor by combining state transition features with communication behavior feature matrix, and obtain a hybrid feature tensor by tensor decomposition. The fourth unit is used to perform temporal segmentation on the hybrid feature tensor, calculate the feature correlation score in each time segment, construct the feature correlation score into a graph structure, and extract the temporal graph features through a graph convolutional network; The fifth unit is used to train a Siamese neural network with the temporal graph features of historical anomaly samples as input, and to construct a discrimination boundary through contrastive learning to obtain an anomaly detector; The sixth unit is used to input the time-series graph features generated from real-time data into the anomaly detector to calculate the anomaly score; The seventh unit is used to send an alarm message to the power control system and record the abnormal behavior characteristics when the abnormal score exceeds the adaptive threshold.

[0093] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0094] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0095] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0096] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for detecting abnormal behavior in power industrial control systems based on deep learning, characterized in that, include: Collect communication messages, equipment status, and operation commands from the power industrial control system, and perform timing alignment and standardization processing; Extract byte sequence features from communication messages, map each byte to a high-dimensional space, compress the message fingerprint vector using an autoencoder, and construct a communication behavior feature matrix based on the message fingerprint vector; The time difference sequence of device state is calculated, state transition features are identified based on recursive decomposition algorithm, and a feature association tensor is constructed by combining the state transition features with the communication behavior feature matrix. The hybrid feature tensor is obtained by tensor decomposition. Temporal segmentation is performed on the hybrid feature tensor, and feature correlation scores are calculated within each time segment. The feature correlation scores are constructed into a graph structure, and temporal graph features are extracted through a graph convolutional network. The Siamese neural network is trained using the temporal graph features of historical anomaly samples as input, and the discrimination boundary is constructed through comparative learning to obtain the anomaly detector; The time-series graph features generated from real-time data are input into the anomaly detector to calculate the anomaly score; When the abnormal score exceeds the adaptive threshold, an alarm message is sent to the power control system and the abnormal behavior characteristics are recorded.

2. The method according to claim 1, characterized in that, Byte sequence features are extracted from communication messages, each byte is mapped to a high-dimensional space, and a message fingerprint vector is obtained through autoencoder compression. A communication behavior feature matrix is ​​constructed based on the message fingerprint vector, including: The communication message is segmented according to the functional field and the data field, the byte sequence is extracted, the byte sequence is sampled in segments, and the frequency domain features are obtained through Fourier transform; By sequentially traversing each byte in the byte sequence, a byte co-occurrence matrix is ​​constructed to record the occurrence frequency and transition probability of adjacent bytes, thereby obtaining the byte sequence features; An overlapping sliding window is set on the byte sequence, and multiple local byte sequences are obtained by moving the overlapping sliding window. The entropy feature of each local byte sequence is calculated, and the entropy feature, the frequency domain feature and the byte sequence feature are fused to obtain the sequence statistical feature. Correlation analysis is performed on the sequence statistical features, and a subset of features is selected based on mutual information values ​​and mapped to a high-dimensional space. Location information and protocol information are then added to the high-dimensional space. A context window is set on the byte sequence, the association relationship between bytes is calculated based on the context window, and the association relationship is added to the mapping result in the high-dimensional space; Using a pre-trained autoencoder network, the mapping result of the high-dimensional space is compressed to obtain an initial fingerprint vector, and the feature dimension with high variance contribution rate is selected to construct the message fingerprint vector. Multiple message fingerprint vectors are collected based on temporal relationships, and the similarity features of adjacent vectors are calculated to construct a communication behavior feature matrix.

3. The method according to claim 1, characterized in that, The time difference sequence of device states is calculated, and state transition features are identified based on a recursive decomposition algorithm. These state transition features are then combined with the communication behavior feature matrix to construct a feature correlation tensor. Tensor decomposition yields a hybrid feature tensor, including: A multidimensional state matrix is ​​constructed, wherein the row vectors of the multidimensional state matrix correspond to the device state parameters and the column vectors correspond to the time series. The time difference sequence is obtained by calculating the difference between adjacent time points in the multidimensional state matrix. The recursive decomposition algorithm is applied to the time difference sequence to obtain the intrinsic mode function set. The Hilbert transform of the intrinsic mode function set is calculated to obtain the instantaneous frequency features. The power spectral density distribution based on the instantaneous frequency features is used to identify state transition points. The state transition points are used as nodes to construct a state transition graph. The weights of the edges between nodes are calculated based on the changes in physical quantities during the state transition process. The degree distribution and centrality features of the state transition graph are extracted to obtain the state transition features. Calculate the adjacency relationship between nodes in the state transition graph and the conditional mutual information between the communication behavior feature matrix, and construct a feature association tensor. The three dimensions of the feature association tensor correspond to the state transition features, communication behavior features, and time series, respectively. The feature association tensor is decomposed into the product of multiple tensor kernels using a tensor network, and a hybrid feature tensor is obtained based on the association between the tensor kernels.

4. The method according to claim 3, characterized in that, The feature association tensor is decomposed into a product of multiple tensor kernels using a tensor network. Based on the association between the tensor kernels, a hybrid feature tensor is calculated, including: The feature-associated tensor is optimized through tensor linear decomposition to construct an optimization objective function containing tensor dimension constraints. Multiple tensor kernels are obtained through iterative calculation based on the index shrinkage rule. The connection dimension between the tensor kernels is calculated. The connection topology of the tensor kernels is determined according to the connection dimension and the optimization objective function to form a tensor network. The shrinkage complexity matrix is ​​calculated for the tensor kernel in the tensor network. Based on the shrinkage complexity matrix, the optimal shrinkage path is calculated using a dynamic programming algorithm. Based on the optimal shrinkage path, sub-tensor blocks are divided. Parallel shrinkage calculation is performed on the sub-tensor blocks to obtain intermediate tensor results. Construct the density matrix corresponding to the intermediate tensor result, calculate the reduced density matrix of the density matrix, solve the eigenvalue spectrum of the reduced density matrix, and calculate the entanglement entropy between tensor kernels based on the eigenvalue spectrum; The feature weight coefficients are determined based on the entanglement entropy, and the feature weight coefficients are combined with the tensor kernel in a weighted combination operation to generate a hybrid feature tensor. Tensor dimensionality reduction analysis is performed on the hybrid feature tensor to extract the principal feature components. Based on the principal feature components, the information retention rate and feature discriminant are calculated to generate an optimized hybrid feature tensor.

5. The method according to claim 1, characterized in that, Temporal segmentation is performed on the hybrid feature tensor, and feature correlation scores are calculated within each time segment. These feature correlation scores are then constructed into a graph structure, and temporal graph features are extracted using a graph convolutional network, including: Temporal segmentation is performed on the hybrid feature tensor. The size of the dynamic time window is determined based on the variance variation trend of the feature vectors in the hybrid feature tensor. The hybrid feature tensor is divided into multiple time segments according to the size of the dynamic time window, and there is a preset overlap interval between adjacent time segments. For each time segment, a set of feature pairs is constructed, the conditional mutual information and nonlinear correlation coefficient between features in the set of feature pairs are calculated, and the results are fused to obtain a feature correlation score matrix; The feature correlation score matrix is ​​constructed into a graph structure. The nodes in the graph structure are pre-screened based on the locality sensitive hashing algorithm. The graph structure is then optimized using the spectral clustering method to obtain an optimized graph structure containing multiple communities. Based on the community division of the optimized graph structure, node information is propagated in the graph structure through a graph convolutional network with multi-head attention mechanism and a message passing function. A corresponding receptive field range is set for each community. Within the receptive field range, the node features within the community are fused, and the node features between communities are fused based on the connection relationship between communities. Combined with the temporal dependency relationship, the temporal graph features are obtained.

6. The method according to claim 5, characterized in that, The feature correlation score matrix is ​​constructed into a graph structure. Nodes in the graph structure are pre-screened using the Locality Sensitive Hashing algorithm. The graph structure is then optimized using spectral clustering to obtain an optimized graph structure containing multiple communities, including: Multiple random mapping transformations are performed on the feature vectors of nodes in the graph structure, the mapping results are quantized into binary codes, and the node similarity matrix is ​​obtained by calculating the Hamming distance of the binary codes. Calculate the Laplacian matrix eigenvalues ​​of the node similarity matrix, determine the optimal number of communities based on the difference of the Laplacian matrix eigenvalues, and select the eigenvectors corresponding to the optimal number of communities to construct a low-dimensional node representation; The neighborhood density and relative distance of each node are calculated using the low-dimensional node representation. The community center node is identified through the density clustering effect, and the initial center point of the community is obtained. Other nodes are assigned to the nearest initial community center point, and the structural distance and attribute distance between the nodes and the initial community center point are calculated to obtain the initial community division results. Perform an edge pruning operation on the initial community division result, remove edges whose internal connection weights are lower than the preset weight threshold, determine the remaining connection edges, and merge communities whose number of remaining connection edges is less than the preset connection threshold to obtain an optimized community structure. Based on the connection relationships of the reconstructed graph structure according to the optimized community structure, an optimized graph structure containing multiple communities is obtained.

7. The method according to claim 1, characterized in that, The Siamese neural network is trained using the temporal graph features of historical anomaly samples as input. A discrimination boundary is constructed through contrastive learning, resulting in anomaly detectors including: The temporal graph features of historical anomaly samples are input into two sub-networks of the Siamese neural network. The two sub-networks have the same network structure and shared weight parameters. Each sub-network maps the input temporal graph features into a contrast feature vector. Construct a contrastive loss function, which is used to maximize the mutual information between the contrastive feature vectors of similar abnormal samples and minimize the mutual information between the contrastive feature vectors of different types of abnormal samples; Based on the contrastive loss function, a Siamese neural network is trained to obtain a feature extractor for extracting discriminative features; The feature extractor is used to extract the temporal spectrum features of the abnormal samples, and a support vector machine classifier is constructed in the feature space. The classification hyperplane of the support vector machine classifier is the anomaly discrimination boundary. The feature extractor and the anomaly detection boundary are integrated into an anomaly detector.

8. A deep learning-based abnormal behavior detection system for power industrial control systems, used to implement the method of any one of claims 1-7, characterized in that, include: The first unit is used to collect communication messages, equipment status and operation instructions of the power industrial control system, and to perform timing alignment and standardization processing. The second unit is used to extract byte sequence features from communication messages, map each byte to a high-dimensional space, compress it through an autoencoder to obtain a message fingerprint vector, and construct a communication behavior feature matrix based on the message fingerprint vector; The third unit is used to calculate the time difference sequence of device status, identify state transition features based on recursive decomposition algorithm, construct a feature association tensor by combining state transition features with communication behavior feature matrix, and obtain a hybrid feature tensor by tensor decomposition. The fourth unit is used to perform temporal segmentation on the hybrid feature tensor, calculate the feature correlation score in each time segment, construct the feature correlation score into a graph structure, and extract the temporal graph features through a graph convolutional network; The fifth unit is used to train a Siamese neural network with the temporal graph features of historical anomaly samples as input, and to construct a discrimination boundary through contrastive learning to obtain an anomaly detector; The sixth unit is used to input the time-series graph features generated from real-time data into the anomaly detector to calculate the anomaly score; The seventh unit is used to send an alarm message to the power control system and record the abnormal behavior characteristics when the abnormal score exceeds the adaptive threshold.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images