Power grid equipment fault diagnosis and life prediction method and system

Abstract

This invention provides a method and system for fault diagnosis and lifespan prediction of power transmission network equipment, belonging to the field of fault diagnosis technology. The method includes acquiring time-series electrical characteristic data from multiple monitoring nodes, performing window segmentation and statistical feature extraction, establishing a dynamic correlation graph structure based on a spatiotemporal correlation constraint model, calculating the anomaly contribution of each node and marking candidate anomaly nodes, determining the fault propagation path through reverse tracing and path analysis, and finally outputting the fault location result. This invention can accurately identify anomaly nodes, precisely track fault propagation paths, and improve the accuracy and timeliness of power transmission network fault diagnosis.

Description

Technical Field

[0001] This invention relates to the field of fault diagnosis technology, and in particular to a method and system for fault diagnosis and life prediction of power transmission network equipment. Background Technology

[0002] The power transmission network is a crucial component of the power system, undertaking the task of long-distance power transmission. Its safe and stable operation is vital to the national economy and people's lives. With the rapid development and widespread application of smart grids, the scale of power transmission network equipment is constantly expanding and its structure is becoming increasingly complex. Fault diagnosis and life prediction of various equipment have become key technologies for ensuring the safe and stable operation of the power grid.

[0003] Currently, power transmission network equipment monitoring technology has gradually evolved from traditional periodic manual inspections to real-time online monitoring systems. Various sensors and monitoring devices are widely deployed at key nodes such as transmission lines and substations to collect electrical characteristic data such as voltage, current, temperature, and vibration in real time. This data contains rich information on equipment operating status, providing a data foundation for equipment fault diagnosis and lifespan prediction.

[0004] Traditional fault diagnosis methods for power transmission network equipment mainly rely on expert experience and simple threshold judgments. By setting safety thresholds for various electrical parameters, alarms are triggered when monitored data exceeds these thresholds. With the development of artificial intelligence and big data technologies, machine learning-based fault diagnosis methods are increasingly being applied to power systems. By establishing data-driven fault characteristic models, the accuracy and timeliness of fault diagnosis have been improved.

[0005] Existing technologies for processing time-series electrical characteristic data are relatively simple, mostly using fixed thresholds or simple statistical features for judgment, which cannot fully explore the deep information and changing trends in time series data, resulting in insufficient ability to identify gradual faults and early fault symptoms.

[0006] Existing methods often analyze each monitoring node as an independent entity, ignoring the topological correlation and fault propagation characteristics between various devices in the power transmission network. This makes it difficult to effectively capture the transmission relationship and chain reaction of faults between devices, and it is impossible to accurately locate the source of the fault, which can easily lead to misjudgment or missed judgment.

[0007] Existing fault diagnosis methods lack a quantitative assessment mechanism for abnormal states, making it difficult to distinguish the severity and scope of impact of different abnormal states. They also cannot accurately track the propagation path and development trend of faults, and cannot provide a scientific basis for equipment maintenance and replacement decisions. Summary of the Invention

[0008] The present invention provides a method and system for fault diagnosis and life prediction of power transmission network equipment, which can solve the problems in the prior art.

[0009] A first aspect of the present invention provides a method for fault diagnosis and lifespan prediction of power transmission network equipment, comprising:

[0010] Acquire time-series electrical characteristic data from multiple monitoring nodes of the power transmission network, wherein the time-series electrical characteristic data contains electrical physical quantities in multiple time-series dimensions;

[0011] The time-series electrical characteristic data is segmented into time windows, and statistical features are extracted from the electrical physical quantities within each time window to obtain a segmented statistical feature set.

[0012] Based on the spatiotemporal correlation constraint model, cross-port correlation analysis is performed on the segmented statistical feature set to establish a dynamic correlation graph structure between different monitoring nodes. The edge weights of the dynamic correlation graph structure reflect the synchronicity and causal transmission of changes in electrical physical quantities between adjacent monitoring nodes, thus obtaining spatiotemporal correlation features.

[0013] Based on the deviation between the local topological change and the global topological stability of each monitoring node in the spatiotemporal correlation features, the abnormal contribution of each monitoring node is calculated, and monitoring nodes whose abnormal contribution exceeds the preset judgment conditions are marked as candidate abnormal nodes.

[0014] Based on the positional relationship of the candidate abnormal nodes in the dynamic association graph structure and the changing direction of the edge weights, the fault propagation path is determined through reverse tracing and path analysis;

[0015] The fault location result is output based on the starting node position of the fault propagation path.

[0016] Based on the spatiotemporal correlation constraint model, cross-port correlation analysis is performed on the segmented statistical feature set to establish a dynamic correlation graph structure between different monitoring nodes, resulting in spatiotemporal correlation features including:

[0017] Calculate the similarity of the electrical physical quantity change trends of any two monitoring nodes in the segmented statistical feature set within the same time window, and the time series-based delayed correlation of electrical physical quantity changes between different monitoring nodes, to obtain the inter-node correlation metric;

[0018] An initial association graph is constructed based on the inter-node association metric. Edges in the initial association graph that do not satisfy physical connectivity are removed according to the topological connection constraints of the power transmission network. At the same time, the retained edges are assigned directional attributes according to the transmission directionality of electrical physical quantities between adjacent monitoring nodes, thus obtaining a topologically constrained directed association graph.

[0019] For each edge in the directed association graph after the topological constraints, the edge weight is calculated. The edge weight is obtained by weighted fusion of the synchronicity measure and the causal transitivity measure of the changes in electrical physical quantities between the two connected monitoring nodes.

[0020] The directed association graph with the topological constraints is combined with the edge weights to form the dynamic association graph structure, and the degree centrality, betweenness centrality and edge weight distribution characteristics of each monitoring node in the dynamic association graph structure are extracted as the spatiotemporal association characteristics.

[0021] The edge weight is obtained by weighted fusion of the synchronicity measure and the causal transitivity measure of the changes in electrical physical quantities between the two connected monitoring nodes, including:

[0022] For two connected monitoring nodes, the waveforms of changes in electrical physical quantities of the two monitoring nodes within the same time window are extracted. The morphological similarity and the time difference of the peak arrival time between the waveforms of changes in electrical physical quantities are calculated. The synchronization metric is obtained based on the normalized reciprocal relationship between the morphological similarity and the time difference.

[0023] For two connected monitoring nodes, the temporal sequence of the changes in electrical physical quantities between adjacent monitoring nodes is identified, and the transmission direction of the changes in electrical physical quantities is determined according to the temporal sequence. The response intensity ratio between the amplitudes of the changes in electrical physical quantities along the transmission direction of adjacent monitoring nodes is calculated, and the response intensity ratio is used as the causal transitivity metric.

[0024] The synchronization metric and the causal transitivity metric are respectively assigned adaptive weighting coefficients, which are dynamically adjusted according to the electrical distance and load transmission capacity of the two connected monitoring nodes in the power grid topology.

[0025] The edge weight is obtained by weighting and summing the synchronicity metric and the causal transitivity metric according to the adaptive weight coefficient.

[0026] Based on the deviation between the local topological changes and the global topological stability of each monitoring node in the spatiotemporal correlation features, the anomaly contribution of each monitoring node is calculated, including:

[0027] For each monitoring node in the spatiotemporal correlation features, a topological feature evolution trajectory of the monitoring node is constructed within multiple consecutive time windows. The topological feature evolution trajectory includes the degree centrality sequence, betweenness centrality sequence, and connection edge weight sequence of the monitoring node. The trajectory deviation distance of the topological feature evolution trajectory relative to the historical steady-state trajectory is calculated to obtain the local topological change amount.

[0028] A global topology feature matrix is ​​constructed based on the dynamic correlation graph structure. Singular value decomposition is performed on the global topology feature matrix to obtain the dominant topology pattern and the secondary topology pattern. The energy distribution ratio between the projection coefficient of each monitoring node in the dominant topology pattern and the projection coefficient in the secondary topology pattern is calculated. When the energy distribution ratio is reversed, the degree of deviation of the monitoring node from the global topology stability is determined.

[0029] The anomaly contribution of each monitoring node is obtained by coupling the trajectory deviation distance in the local topology change with the energy distribution ratio reversal amplitude in the deviation degree of the global topology stability.

[0030] Singular value decomposition is performed on the global topology feature matrix to obtain the dominant topology pattern and the secondary topology pattern. The energy distribution ratio between the projection coefficients of each monitoring node in the dominant topology pattern and the projection coefficients in the secondary topology pattern is calculated, including:

[0031] The global topological feature matrix is ​​constructed based on the degree centrality, betweenness centrality, and edge weight distribution characteristics of all monitoring nodes in the dynamic association graph structure. The rows of the global topological feature matrix correspond to the monitoring nodes, and the columns correspond to the topological feature dimensions. The global topological feature matrix is ​​decomposed by singular value to obtain the left singular vector matrix, the singular value diagonal matrix, and the right singular vector matrix.

[0032] The singular values ​​in the singular value diagonal matrix are sorted according to their energy contribution rates. The left singular vectors corresponding to the first few singular values ​​whose cumulative energy contribution rates reach the preset energy concentration are selected as the dominant topology pattern. The left singular vectors whose energy contribution rates are after the preset energy concentration and whose singular values ​​are non-zero are selected as the secondary topology pattern.

[0033] For each monitoring node, the row vector corresponding to the monitoring node in the global topology feature matrix is ​​extracted, and the row vector is projected into the vector spaces of the dominant topology pattern and the secondary topology pattern respectively. The square of the modulus of the projection vector of the row vector in the vector space of the dominant topology pattern is calculated as the projection coefficient in the dominant topology pattern, and the square of the modulus of the projection vector of the row vector in the vector space of the secondary topology pattern is calculated as the projection coefficient in the secondary topology pattern.

[0034] The energy distribution ratio is obtained by calculating the ratio of the projection coefficient in the dominant topology mode to the projection coefficient in the secondary topology mode.

[0035] Based on the positional relationships of the candidate anomaly nodes in the dynamic association graph structure and the changing direction of the edge weights, the fault propagation path is determined through reverse tracing and path analysis, including:

[0036] Extract all incoming and outgoing edges of the candidate abnormal nodes in the candidate abnormal node set in the dynamic association graph structure, calculate the change magnitude of edge weight of each incoming edge between the current time window and the historical time window, determine the change direction of the edge weight according to the sign of the change magnitude of the edge weight, and select the incoming edge with the largest change magnitude of edge weight and the change direction of enhancement as the main import edge.

[0037] Using the source node of the main import edge as the current trace node, determine whether the current trace node belongs to the candidate abnormal node set. If it does, repeat the operation of extracting the main import edge and use the source node of the main import edge as the new current trace node. If it does not belong, mark the current trace node as the trace termination node.

[0038] Record all node sequences and connection edge sequences traversed during the reverse tracing process from the initial candidate abnormal node to the tracing termination node. Arrange the node sequences in reverse order of the reverse tracing to obtain the forward node sequence, and arrange the connection edge sequences in reverse order of the reverse tracing to obtain the forward connection edge sequence.

[0039] Based on the forward node sequence and the forward connection edge sequence, a directed path is constructed from the trace termination node to the initial candidate abnormal node, which serves as the fault propagation path.

[0040] A second aspect of the present invention provides a fault diagnosis and life prediction system for power transmission network equipment, comprising:

[0041] The first unit is used to acquire time-series electrical characteristic data from multiple monitoring nodes of the power transmission network, wherein the time-series electrical characteristic data contains electrical physical quantities in multiple time-series dimensions;

[0042] The second unit is used to segment the time-series electrical characteristic data into time windows and extract statistical features of the electrical physical quantities within each time window to obtain a segmented statistical feature set.

[0043] The third unit is used to perform cross-port correlation analysis on the segmented statistical feature set based on the spatiotemporal correlation constraint model, establish a dynamic correlation graph structure between different monitoring nodes, and the edge weights of the dynamic correlation graph structure reflect the synchronicity and causal transmission of changes in electrical physical quantities between adjacent monitoring nodes, thereby obtaining spatiotemporal correlation features.

[0044] The fourth unit is used to calculate the abnormal contribution of each monitoring node based on the deviation between the local topological change and the global topological stability of each monitoring node in the spatiotemporal correlation features, and to mark the monitoring nodes whose abnormal contribution exceeds the preset judgment conditions as candidate abnormal nodes.

[0045] The fifth unit is used to determine the fault propagation path by reverse tracing and path analysis based on the positional relationship of the candidate abnormal nodes in the dynamic association graph structure and the changing direction of the edge weights.

[0046] The sixth unit is used to output the fault location result based on the starting node position of the fault propagation path.

[0047] A third aspect of the present invention provides an electronic device, comprising:

[0048] processor;

[0049] Memory used to store processor-executable instructions;

[0050] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0051] 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.

[0052] The beneficial effects of this application are as follows:

[0053] This invention segments the time-series electrical characteristic data of multiple monitoring nodes in the power transmission network into time windows and extracts statistical features. It then establishes a dynamic correlation graph structure by combining a spatiotemporal correlation constraint model. This approach can accurately capture the electrical correlation between devices and the fault propagation pattern, effectively improving the accuracy and reliability of fault diagnosis.

[0054] This invention utilizes a dynamic correlation graph structure to analyze the local topology changes and global stability of monitoring nodes, calculates the anomaly contribution and marks candidate anomaly nodes, and determines the fault propagation path through reverse tracing and path analysis, thereby achieving rapid fault location and root cause analysis, and significantly reducing fault troubleshooting time and maintenance costs.

[0055] The method of this invention combines time-series data analysis with dynamic evolution of graph structures, which can not only monitor equipment status in real time, but also predict potential fault risks and equipment lifespan, providing a scientific basis for power grid operation and maintenance, improving the overall operational safety and reliability of the power transmission network, and having significant economic and social benefits. Attached Figure Description

[0056] Figure 1 This is a flowchart illustrating the method for fault diagnosis and life prediction of power transmission network equipment according to an embodiment of the present invention;

[0057] Figure 2 This is a flowchart illustrating the parallel process for calculating the anomaly contribution of a power transmission network monitoring node in an embodiment of the present invention. Detailed Implementation

[0058] 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.

[0059] 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.

[0060] Figure 1 This is a flowchart illustrating the fault diagnosis and lifespan prediction method for power transmission network equipment according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0061] Acquire time-series electrical characteristic data from multiple monitoring nodes of the power transmission network, wherein the time-series electrical characteristic data contains electrical physical quantities in multiple time-series dimensions;

[0062] The time-series electrical characteristic data is segmented into time windows, and statistical features are extracted from the electrical physical quantities within each time window to obtain a segmented statistical feature set.

[0063] Based on the spatiotemporal correlation constraint model, cross-port correlation analysis is performed on the segmented statistical feature set to establish a dynamic correlation graph structure between different monitoring nodes. The edge weights of the dynamic correlation graph structure reflect the synchronicity and causal transmission of changes in electrical physical quantities between adjacent monitoring nodes, thus obtaining spatiotemporal correlation features.

[0064] Based on the deviation between the local topological change and the global topological stability of each monitoring node in the spatiotemporal correlation features, the abnormal contribution of each monitoring node is calculated, and monitoring nodes whose abnormal contribution exceeds the preset judgment conditions are marked as candidate abnormal nodes.

[0065] Based on the positional relationship of the candidate abnormal nodes in the dynamic association graph structure and the changing direction of the edge weights, the fault propagation path is determined through reverse tracing and path analysis;

[0066] The fault location result is output based on the starting node position of the fault propagation path.

[0067] In one optional implementation, cross-port correlation analysis is performed on the segmented statistical feature set based on a spatiotemporal correlation constraint model to establish a dynamic correlation graph structure between different monitoring nodes, resulting in spatiotemporal correlation features including:

[0068] Calculate the similarity of the electrical physical quantity change trends of any two monitoring nodes in the segmented statistical feature set within the same time window, and the time series-based delayed correlation of electrical physical quantity changes between different monitoring nodes, to obtain the inter-node correlation metric;

[0069] An initial association graph is constructed based on the inter-node association metric. Edges in the initial association graph that do not satisfy physical connectivity are removed according to the topological connection constraints of the power transmission network. At the same time, the retained edges are assigned directional attributes according to the transmission directionality of electrical physical quantities between adjacent monitoring nodes, thus obtaining a topologically constrained directed association graph.

[0070] For each edge in the directed association graph after the topological constraints, the edge weight is calculated. The edge weight is obtained by weighted fusion of the synchronicity measure and the causal transitivity measure of the changes in electrical physical quantities between the two connected monitoring nodes.

[0071] The directed association graph with the topological constraints is combined with the edge weights to form the dynamic association graph structure, and the degree centrality, betweenness centrality and edge weight distribution characteristics of each monitoring node in the dynamic association graph structure are extracted as the spatiotemporal association characteristics.

[0072] When performing cross-port correlation analysis on segmented statistical feature sets based on a spatiotemporal correlation constraint model, the similarity calculation module first receives segmented statistical feature data from each monitoring node. Each node contains a sixteen-dimensional feature vector within a single time window, including the mean voltage amplitude, standard deviation of voltage amplitude, effective current value, current phase angle, active power fluctuation coefficient, reactive power change rate, frequency offset, and harmonic distortion rate. The module calculates the Pearson correlation coefficient for the feature vectors of any two monitoring nodes within the same time window, and analyzes the consistency of the changing trends by comparing each element. During the calculation process, the two feature vectors are first normalized to zero mean, and then the covariance is calculated and divided by the standard deviation product to obtain a similarity value ranging from -1 to 1. When the absolute similarity value is greater than 0.75, the two nodes are considered to have a strong correlation; greater than 0.5 and less than 0.75 is considered a moderate correlation; and less than 0.75 is considered a weak correlation.

[0073] The delay correlation analysis module processes the time-series propagation relationships of electrical physical quantity changes between different monitoring nodes in parallel. The module maintains a sliding window queue with a length of 15 time steps, each time step corresponding to a 1-minute sampling interval. For each pair of monitoring nodes, the module extracts its voltage amplitude change sequence, current amplitude change sequence, and power change sequence as the analysis objects. Delay correlation calculation employs a time-off cross-correlation analysis method, using the electrical physical quantity change sequence of the first node as a reference signal and performing correlation matching with the change sequence of the second node at delays of 0 to 8 time steps. During the calculation, for each delay time point, the correlation coefficient between the two sequences is calculated, and the maximum correlation coefficient and its corresponding delay time are recorded. When the maximum correlation coefficient exceeds 0.6 and the delay time meets the propagation time constraint of electrical signals in the transmission line, a causal correlation relationship is confirmed. The inter-node correlation metric is obtained by weighted fusion of synchronization similarity and delay correlation coefficient, with fusion weights set to 0.6 and 0.4 respectively. The final correlation metric value ranges from 0 to 1.

[0074] The initial association graph construction module creates a graph structure representation based on the calculated association metrics between nodes. The graph structure uses an adjacency list for storage, with each monitored node corresponding to a vertex in the list. When the association metric between two nodes exceeds a preset threshold of 0.4, a connection edge is established between the corresponding vertices. The initial association graph is undirected; the existence of an edge only indicates a statistically significant association between two nodes. Each edge entry in the adjacency list contains detailed information such as the target node identifier, association metric value, similarity component value, delayed association component value, and delay time. The module also maintains a global edge weight matrix, with a size equal to the square of the number of nodes, for quickly querying the association metric between any two nodes.

[0075] The topology connectivity constraint processing module receives static topology information of the transmission network as physical constraints. This topology information includes complete power grid structure data such as the start and end node identifiers of transmission lines, high-voltage and low-voltage side connection points of transformers, the current status of circuit breakers and disconnectors, and bus segmentation information. The constraint processing procedure uses a graph traversal algorithm to check whether there is a valid electrical connection path between the two monitoring nodes corresponding to each edge in the initial association graph. The checking method employs a breadth-first search algorithm, starting from the first node and searching for all paths to the second node in the physical topology graph, with path length limited to five hops. Edges that do not satisfy physical connectivity are marked as non-physical edges and removed from the initial association graph. During removal, the module records the number of removed edges and the corresponding node pair information for subsequent algorithm performance evaluation.

[0076] The direction attribute assignment module directs the retained edges based on the directionality of electrical physical quantities. Direction determination is based on power flow analysis; the module extracts active and reactive power measurement data for the two connected nodes within the observation time window. By calculating the net power transfer between the two nodes, the main power flow direction is determined. When the net active power transfer from node A to node B is positive and lasts for more than 60% of the observation window, the edge direction points from A to B. For cases where the power flow direction changes frequently within the time window, the module calculates the cumulative time of forward and reverse power transfer and selects the direction with the longer cumulative time as the main direction attribute of the edge. For edges with relatively balanced power flows, the bidirectional attribute is maintained, and the weights of the two directions are calculated separately in subsequent processing. After processing, the topologically constrained directed graph ensures both the rationality of physical connections and provides directional information for causal relationship analysis.

[0077] The edge weight calculation module assigns detailed weights to each edge in the directed graph. Synchronization metrics are calculated by analyzing the consistency of changes in electrical physical quantities between two connected monitoring nodes within the same time window. The calculation method first aligns the multidimensional feature vectors of the two nodes in time, then calculates the cosine similarity of the normalized feature vectors. During the cosine similarity calculation, the dot product of the feature vectors is divided by the product of the magnitudes of the two vectors, yielding a synchronization metric value ranging from [-1, 1]. Causality transitivity metrics are calculated based on the Granger causality test time series analysis method.

[0078] The module constructs a vector autoregressive model containing historical electrical physical quantity data for two nodes. By comparing the difference in prediction accuracy of the second node's future value with and without information from the first node, the causal influence of the first node on the second node is quantified. The causal transitivity metric is obtained through F-statistic transformation, with a value range of [0,1]. The edge weights are calculated by weighted linear combination of the synchronicity metric and the causal transitivity metric, with weight coefficients of 0.55 and 0.45, respectively, resulting in a final edge weight value range of [0,1].

[0079] The dynamic association graph structure formation module integrates the processed directed association graph with the calculated edge weights. The graph structure uses a weighted directed graph data structure for storage, including core components such as node sets, edge sets, and edge weight mapping tables. The node set records basic information such as the identifier, geographical location, and device type of all monitored nodes. The edge set stores detailed attributes of each edge, including the starting node, ending node, edge weight, and creation time. The edge weight mapping table provides a fast query interface from node pairs to edge weights and supports real-time weight update operations.

[0080] The spatiotemporal correlation feature extraction module calculates network topology feature indicators for each monitored node from the dynamic correlation graph structure. Degree centrality calculation consists of two components: in-degree centrality and out-degree centrality. In-degree centrality equals the number of edges pointing to that node divided by the total number of nodes in the graph minus 1; out-degree centrality equals the number of edges originating from that node divided by the total number of nodes in the graph minus 1. Betweenness centrality calculation uses the Floyd algorithm to find the shortest path between all pairs of nodes in the graph, counts the number of shortest paths passing through the target node, and then divides this number by the total number of shortest paths in the graph to obtain the normalized betweenness centrality value. Edge weight distribution feature extraction includes statistics such as the arithmetic mean, geometric mean, standard deviation, maximum value, minimum value, and median of all edge weights connected to the node, used to comprehensively describe the strength distribution pattern of the connections around the node.

[0081] In one alternative implementation, the edge weight is obtained by weighted fusion of a synchronicity metric and a causal transitivity metric for changes in electrical physical quantities between the two connected monitoring nodes, including:

[0082] For two connected monitoring nodes, the waveforms of changes in electrical physical quantities of the two monitoring nodes within the same time window are extracted. The morphological similarity and the time difference of the peak arrival time between the waveforms of changes in electrical physical quantities are calculated. The synchronization metric is obtained based on the normalized reciprocal relationship between the morphological similarity and the time difference.

[0083] For two connected monitoring nodes, the temporal sequence of the changes in electrical physical quantities between adjacent monitoring nodes is identified, and the transmission direction of the changes in electrical physical quantities is determined according to the temporal sequence. The response intensity ratio between the amplitudes of the changes in electrical physical quantities along the transmission direction of adjacent monitoring nodes is calculated, and the response intensity ratio is used as the causal transitivity metric.

[0084] The synchronization metric and the causal transitivity metric are respectively assigned adaptive weighting coefficients, which are dynamically adjusted according to the electrical distance and load transmission capacity of the two connected monitoring nodes in the power grid topology.

[0085] The edge weight is obtained by weighting and summing the synchronicity metric and the causal transitivity metric according to the adaptive weight coefficient.

[0086] When the edge weight calculation module extracts waveforms of electrical physical quantity changes for two connected monitoring nodes, the waveform data acquisition interface obtains raw measurement data within a specified time window from the real-time data acquisition server. Each monitoring node provides four types of basic waveform data: voltage amplitude sequence, current amplitude sequence, active power sequence, and reactive power sequence. The sampling frequency is uniformly set to 50Hz, and the length of a single time window is 2 minutes, containing 6000 sampling points. The waveform preprocessing module performs outlier detection and smoothing filtering on the raw data. Outlier detection uses the 3σ criterion, and data points exceeding the normal range are repaired using linear interpolation. Smoothing filtering uses a Hamming window function with a window length of 7 sampling points. The filtered waveform data is stored in a memory buffer for subsequent processing.

[0087] The morphological similarity calculation unit uses a dynamic time warping algorithm to compare and analyze the waveform changes of electrical physical quantities between two monitoring nodes. During the algorithm implementation, a distance matrix is ​​constructed, with a size of 6000×6000. Each matrix element calculates the Euclidean distance between the corresponding time points of the two waveforms. The dynamic programming algorithm starts from the top left corner of the matrix and searches for the optimal alignment path according to the principle of minimum cumulative distance. Path search constraints include monotonicity and continuity constraints to ensure that the causal relationship of the time series is not disrupted. The cumulative distance of the optimal path divided by the path length yields the average alignment distance. Morphological similarity is obtained through a negative exponential function transformation, with the transformation formula being e^(-∞ / ∞). -平均对齐距离 Ensure that the similarity value is within the range of [0,1].

[0088] The peak arrival time detection module employs an adaptive threshold peak recognition algorithm to process the waveform of electrical physical quantity changes. The threshold is set to twice the standard deviation of the waveform within the current time window. Peak detection also requires that the peak point has local maxima within five sampling points to its left and right. For identified candidate peak points, they are sorted in descending order of peak amplitude, and the top five most significant peaks are selected as feature peaks. The time difference is calculated by matching the feature peaks of two nodes using a paired nearest neighbor algorithm. The matching principle is that the relative error of the peak amplitude is less than 20% and the time difference is less than 10 seconds. The arrival time difference of successfully matched peak pairs is calculated, and the arithmetic mean of the time differences of all peak pairs is taken as the overall time difference. The normalization reciprocal relationship is achieved by 1 / (1 + time difference), ensuring that the smaller the time difference, the stronger the synchronization.

[0089] The electrical physical quantity change timing identification module determines the moment of significant change through a rate of change detection algorithm. The rate of change is calculated using the central difference method, with the formula (f(t+Δt)-f(t-Δt)) / (2Δt). The change threshold is set to 5 times the standard deviation of steady-state fluctuations. When the absolute value of the rate of change of three consecutive sampling points exceeds the threshold, it is marked as the moment of change. Steady-state interval identification is achieved through sliding variance analysis, with a sliding window length of 20 sampling points. When the variance within the window is less than 10% of the global variance, the state is considered steady. Precise timing of change is achieved using the second derivative zero-point detection method, which searches for the point of sign change of the second derivative within five sampling points before and after the coarse timing to determine the precise timing of change.

[0090] The time sequence determination module compares the timestamp differences of the changes at two monitoring nodes to determine the transmission direction. The timestamp accuracy is in the millisecond range. When the change at node A is earlier than that at node B and the time difference is between 0.1 s and 5 s, the transmission direction is determined to be from A to B. The propagation delay constraint is determined based on a combination of the geographical and electrical distances between the two nodes. Each km of geographical distance corresponds to approximately 3.3 μs of light propagation delay, and the electrical distance is calculated using the impedance parameters of the transmission lines to determine the electromagnetic wave propagation delay. Cases where changes occur simultaneously or the time difference exceeds a reasonable range are marked as having no clear transmission direction and assigned a neutral weight in subsequent calculations.

[0091] The response intensity calculation module extracts the steady-state values ​​of electrical physical quantities before and after the change for amplitude analysis. The steady-state value before the change is selected as the average of the data within 15 seconds before the change, and the steady-state value after the change is selected as the average of the data within 15 seconds after the change. The validity check of the steady-state interval data is achieved through analysis of variance, requiring that the variance of the data within the steady-state interval be less than 20% of the variance of the change interval. The response intensity is defined as the relative amplitude of the steady-state value change, calculated as |steady-state value after change - steady-state value before change| / steady-state value before change. The response intensity ratio is obtained by dividing the response intensity of the downstream node by the response intensity of the upstream node, with the ratio range limited to [0.2, 5.0]. Values ​​outside this range are truncated. The causal transitivity metric uses the sigmoid function to normalize the response intensity ratio, with the function parameters set so that a ratio of 1 corresponds to a metric value of 0.5.

[0092] The adaptive weighting coefficient calculation module dynamically adjusts the weight allocation based on electrical distance and load transmission capacity. Electrical distance calculation uses the Dijkstra shortest path algorithm to find the path with minimum impedance in the transmission network impedance diagram; the path weight is the product of the impedance per unit length and the length of the line. Electrical distance normalization involves dividing the calculated impedance distance by the network diameter to obtain a normalization coefficient within the range of [0, 1]. Load transmission capacity data is obtained from dispatch automation equipment and includes parameters such as line rated capacity, current load level, and safety margin. The load factor is calculated as the ratio of the current load to the rated capacity, and load factor normalization ensures that the value is within the range of [0, 1].

[0093] The synchronicity weight in the adaptive weighting coefficient is calculated using a linear interpolation method. Interpolation parameters include the normalized electrical distance, load factor, and node type correction coefficient. The basic weight formula is 0.8 - (normalized electrical distance × 0.3) + (load factor × 0.1). The node type correction coefficient is determined based on the functional positioning of the monitoring node in the power grid: 1.1 for generation nodes, 0.9 for load nodes, and 1.0 for substation nodes. The causal transitivity weight is obtained by subtracting the synchronicity weight from 1, ensuring that the sum of the two weights always equals 1. The weight update cycle is set to 5 minutes, and historical weight values ​​are exponentially smoothed with a smoothing coefficient of 0.2 each time an update is performed.

[0094] The edge weight fusion calculation module performs a weighted summation operation on the synchronicity metric and the causal transitivity metric. The weighted summation formula is: synchronicity metric × synchronicity weight + causal transitivity metric × causal transitivity weight. Numerical stability checks are performed during the calculation process, including validating the input metric, verifying the normalization of the weight coefficients, and checking the reasonableness of the output results. Validity verification requires the metric to be within the range [0, 1] and not null, and weight normalization verification requires the error of the weight sum to be less than 0.001. Edge weight post-processing includes smoothing filtering and boundary constraints. The smoothing filter uses a first-order low-pass filter with a cutoff frequency set to 0.1 / min, and the boundary constraints constrain the weights to the range [0.05, 0.95].

[0095] In one optional implementation, the anomaly contribution of each monitoring node is calculated based on the deviation between the local topological change and the global topological stability of each monitoring node in the spatiotemporal correlation features, including:

[0096] For each monitoring node in the spatiotemporal correlation features, a topological feature evolution trajectory of the monitoring node is constructed within multiple consecutive time windows. The topological feature evolution trajectory includes the degree centrality sequence, betweenness centrality sequence, and connection edge weight sequence of the monitoring node. The trajectory deviation distance of the topological feature evolution trajectory relative to the historical steady-state trajectory is calculated to obtain the local topological change amount.

[0097] A global topology feature matrix is ​​constructed based on the dynamic correlation graph structure. Singular value decomposition is performed on the global topology feature matrix to obtain the dominant topology pattern and the secondary topology pattern. The energy distribution ratio between the projection coefficient of each monitoring node in the dominant topology pattern and the projection coefficient in the secondary topology pattern is calculated. When the energy distribution ratio is reversed, the degree of deviation of the monitoring node from the global topology stability is determined.

[0098] The anomaly contribution of each monitoring node is obtained by coupling the trajectory deviation distance in the local topology change with the energy distribution ratio reversal amplitude in the deviation degree of the global topology stability.

[0099] like Figure 2 As shown, the method includes:

[0100] When constructing the topological feature evolution trajectory for each monitoring node in the spatiotemporal correlation features, the anomaly contribution calculation module extracts topological feature data from the time-series data storage within twenty consecutive time windows. Each time window is five minutes long, and the total trajectory time span is one hundred minutes. The degree centrality sequence extraction module obtains the in-degree centrality and out-degree centrality values ​​of each monitoring node within each time window, maintaining a precision of four decimal places and a value range of zero to one. The betweenness centrality sequence extraction module calculates the betweenness centrality of each monitoring node in the dynamic correlation graph and uses the shortest path algorithm to count the proportion of paths passing through that node to the total number of paths. The edge weight sequence extraction module collects the weight values ​​of all edges directly connected to the monitoring node and calculates the arithmetic mean, weighted average, and standard deviation of the edge weights as the edge weight feature values ​​for that time window. The topological feature evolution trajectory is represented by a three-dimensional vector, with vector components representing degree centrality, betweenness centrality, and edge weight feature values, forming a trajectory sequence at twenty time points.

[0101] The historical steady-state trajectory establishment module determines the normal topological feature patterns of each monitoring node by analyzing historical data from the past six months. The steady-state identification algorithm employs a sliding window variance analysis method with a window length of 10 time steps. A steady-state interval is marked when the variance of the topological features within a window is less than 30% of the historical variance for five consecutive windows. The steady-state trajectory calculation obtains the probability density distribution by kernel density estimation of the topological features within all steady-state intervals, selecting the trajectory path with the highest probability density as the historical steady-state trajectory reference benchmark. The trajectory deviation distance calculation uses the Mahalanobis distance algorithm, which considers the correlation and variance differences between various dimensions of the topological features. A covariance matrix is ​​constructed during the distance calculation process, with matrix elements obtained through covariance analysis of historical steady-state data. Matrix inversion uses singular value decomposition to ensure numerical stability. Local topological changes are obtained by time-weighted averaging of the Mahalanobis distances between the current trajectory sequence and the historical steady-state trajectories. The time weights use an exponential decay function with a decay coefficient of 0.9, making deviations from recent time points have a greater impact on the results.

[0102] The global topology feature matrix construction module creates a feature representation matrix based on a dynamic association graph structure. The matrix dimension is the number of monitored nodes multiplied by the number of topology features. Rows represent each monitored node, and columns represent five topology feature indices: degree centrality, betweenness centrality, clustering coefficient, eigenvector centrality, and compact centrality. The clustering coefficient is calculated by analyzing the connection density between a node's neighbors, using the formula: the actual number of connections between neighboring nodes divided by the maximum number of connections between neighboring nodes. Eigenvector centrality uses a power-law iteration algorithm to calculate the principal eigenvectors of the adjacency matrix of the association graph. The convergence condition is that the Euclidean distance between two consecutive iterations is less than 0.001. Compact centrality is obtained by calculating the reciprocal average of the shortest path lengths from the given node to all other nodes. The global topology feature matrix is ​​updated once per time window, and matrix elements are standardized to zero mean to ensure consistent numerical ranges for different feature indices.

[0103] The Singular Value Decomposition (SVD) module performs numerical decomposition on the global topological feature matrix to obtain the dominant and secondary topological modes. SVD employs the Jacobian rotation algorithm, with the decomposition precision set to the square root of the machine precision to ensure the numerical stability of the decomposition results. The decomposition results consist of three components: a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix. The dominant topological mode corresponds to the three largest singular values ​​and their corresponding singular vectors, while the secondary topological modes correspond to the 4th to 8th singular values ​​and their corresponding singular vectors. The energy percentage of each singular value is calculated by dividing the sum of the squares of all singular values ​​by the sum of the squares of all singular values. The dominant mode typically accounts for 70% to 90% of the energy, while the secondary mode accounts for 5% to 20%.

[0104] The projection coefficient calculation module projects the topological feature vectors of each monitoring node into the subspaces of the dominant and secondary topological modes, respectively. Projection calculation uses vector inner product operations; the projection coefficient of the dominant topological mode is the sum of the inner products of the node feature vectors and the dominant mode basis vectors, and the projection coefficient of the secondary topological mode is the sum of the inner products of the node feature vectors and the secondary mode basis vectors. The energy distribution ratio is defined as the sum of the squares of the projection coefficients of the dominant topological mode divided by the sum of the squares of the projection coefficients of the secondary topological mode. Under normal operating conditions, the energy distribution ratio is usually greater than 5; when the ratio is less than 2, it indicates a significant increase in the node's contribution to the secondary mode. Energy distribution ratio reversal is determined by comparing the relative change of the current ratio with the historical average. A reversal is confirmed when the current ratio is less than 40% of the historical average for more than three consecutive time windows. The reversal magnitude is calculated as the difference between the historical average and the current ratio divided by the historical average to obtain the normalized deviation.

[0105] The global topological stability deviation quantification module comprehensively considers the magnitude and duration of energy distribution ratio reversal. The deviation calculation employs a time-decaying cumulative deviation model. Model parameters include reversal magnitude, duration, and time decay coefficient. The cumulative deviation value is obtained by summing the instantaneous deviation magnitudes of each time window multiplied by a time weight. The time weight uses a linear decay function, with the nearest time window having a weight of 1 and the farthest time window having a weight of 0.1. Deviation normalization is achieved through a sigmoid function transformation to ensure the value ranges between 0 and 1. Adjustments to the function parameters ensure that a moderate deviation corresponds to a value of 0.5.

[0106] The anomaly contribution coupling calculation module fuses local topology changes and global topology stability deviations. The coupling calculation employs a nonlinear weighted fusion method, with the fusion weights dynamically adjusted based on the importance of the monitored node in the power grid. Node importance assessment is based on three factors: load capacity, connectivity, and historical fault frequency, with importance scores normalized to the range of 0 to 1. The weight for local topology changes is the node importance score multiplied by 0.6 plus 0.2, while the weight for global topology stability deviation is one minus the local weight. The anomaly contribution calculation formula is the local topology change multiplied by the local weight plus the global deviation multiplied by the global weight. The result is nonlinearly scaled using a cube root transformation to enhance sensitivity to outliers.

[0107] The data storage and update module maintains historical records and real-time status information of abnormal contributions. The storage structure uses a time-series database, and record fields include node identifier, calculation timestamp, local topology change, global deviation, coupling weight, and final abnormal contribution. The data update cycle is synchronized with the topology feature extraction cycle, updating every five minutes. The anomaly detection threshold is set to the historical average abnormal contribution plus three standard deviations. Nodes exceeding this threshold trigger an anomaly alarm and record detailed diagnostic information.

[0108] In one optional implementation, singular value decomposition is performed on the global topology feature matrix to obtain a dominant topology pattern and a secondary topology pattern. The energy distribution ratio between the projection coefficients of each monitoring node in the dominant topology pattern and the projection coefficients in the secondary topology pattern is calculated, including:

[0109] The global topological feature matrix is ​​constructed based on the degree centrality, betweenness centrality, and edge weight distribution characteristics of all monitoring nodes in the dynamic association graph structure. The rows of the global topological feature matrix correspond to the monitoring nodes, and the columns correspond to the topological feature dimensions. The global topological feature matrix is ​​decomposed by singular value to obtain the left singular vector matrix, the singular value diagonal matrix, and the right singular vector matrix.

[0110] The singular values ​​in the singular value diagonal matrix are sorted according to their energy contribution rates. The left singular vectors corresponding to the first few singular values ​​whose cumulative energy contribution rates reach the preset energy concentration are selected as the dominant topology pattern. The left singular vectors whose energy contribution rates are after the preset energy concentration and whose singular values ​​are non-zero are selected as the secondary topology pattern.

[0111] For each monitoring node, the row vector corresponding to the monitoring node in the global topology feature matrix is ​​extracted, and the row vector is projected into the vector spaces of the dominant topology pattern and the secondary topology pattern respectively. The square of the modulus of the projection vector of the row vector in the vector space of the dominant topology pattern is calculated as the projection coefficient in the dominant topology pattern, and the square of the modulus of the projection vector of the row vector in the vector space of the secondary topology pattern is calculated as the projection coefficient in the secondary topology pattern.

[0112] The ratio of the projection coefficient in the dominant topology mode to the projection coefficient in the secondary topology mode is calculated to obtain the energy distribution ratio.

[0113] The global topology feature matrix construction module creates a numerical matrix representation based on the topology feature data of all monitoring nodes in the dynamic association graph structure. The degree centrality feature extraction unit calculates the in-degree centrality and out-degree centrality of each monitoring node from the dynamic association graph. In-degree centrality equals the number of edges pointing to that node divided by the total number of nodes in the graph minus 1, and out-degree centrality equals the number of edges emanating from that node divided by the total number of nodes in the graph minus 1. The betweenness centrality calculation unit uses the Brandes algorithm to calculate the betweenness centrality of each monitoring node. The algorithm achieves efficient computation by accumulating shortest path dependencies, avoiding the overhead of calculating all-pair shortest paths in traditional methods.

[0114] The edge weight distribution feature extraction unit performs statistical analysis on all edge weights connected to each monitoring node, calculating five statistical measures—mean, standard deviation, maximum, minimum, and median—as edge weight distribution features. The global topological feature matrix is ​​constructed by arranging monitoring nodes in numerical order as matrix rows. The topological feature dimensions include eight dimensions: in-degree centrality, out-degree centrality, betweenness centrality, mean edge weight, standard deviation edge weight, maximum edge weight, minimum edge weight, and median edge weight, which are used as matrix columns. The numerical range of matrix elements is normalized using a maximum-minimum normalization method to scale the values ​​of each feature dimension to the [0,1] interval, ensuring the comparability of values ​​for features with different dimensions.

[0115] The singular value decomposition (SVD) module performs numerical decomposition on the global topological feature matrix, resulting in three matrix components. SVD employs a two-stage bidiagonalization algorithm, including Haushold reflection transform to convert the original matrix into a bidiagonal matrix and a QR iterative algorithm to calculate the singular values ​​of the bidiagonal matrix. The algorithm implementation sets the convergence precision to the square root of the machine precision and the maximum number of iterations to 10 times the larger dimension of the matrix. The dimension of the left singular vector matrix is ​​the number of monitoring nodes × the matrix rank, the dimension of the right singular vector matrix is ​​the number of topological feature dimensions × the matrix rank, and the singular value diagonal matrix is ​​a square matrix whose diagonal elements are singular values. Numerical stability is checked during the decomposition process; if the matrix condition number exceeds 10... 12 Regularization is triggered periodically, and the regularization method involves adding a small perturbation term to the diagonal of the matrix. The decomposition results are stored in a compressed sparse row format to reduce storage space requirements and improve the efficiency of subsequent calculations.

[0116] The energy contribution rate calculation module performs energy analysis and sorting operations on the singular values ​​in the singular value diagonal matrix. The energy contribution rate of each singular value is calculated by dividing the square of that singular value by the sum of the squares of all singular values, ensuring that the sum of the energy contribution rates equals 1. The singular values ​​are sorted in descending order of numerical value, and the corresponding left and right singular vectors are rearranged in the same order. The cumulative energy contribution rate is calculated by summing the energy contribution rates of the sorted singular values ​​one by one, and is used to determine the boundary between the dominant and secondary topological patterns. The preset energy concentration parameter is set to 0.85, indicating that the patterns corresponding to the top few singular values ​​with a cumulative energy contribution rate of 85% are selected as the dominant topological patterns. The dominant topological pattern selection algorithm traverses the cumulative energy contribution rate sequence and finds the first index position where the cumulative value exceeds the preset energy concentration. The left singular vectors corresponding to all singular values ​​before this position constitute the vector space basis of the dominant topological pattern.

[0117] The secondary topology pattern selection module processes singular values ​​and their corresponding vectors whose energy contribution rates are below a preset energy concentration level. The selection range for secondary topology patterns includes all left singular vectors whose cumulative energy contribution rate exceeds the preset energy concentration level but whose singular values ​​are not zero. A numerical threshold comparison is used to determine if a singular value is zero; if the singular value is less than the maximum singular value of the matrix × ε × the larger dimension of the matrix, it is considered to have a numerical value of zero. The dimension of the secondary topology pattern vector space is usually significantly smaller than that of the dominant topology pattern, but it contains detailed information about changes in the global topology structure. The orthogonality of the vector space basis is verified by calculating the inner product between the basis vectors; the absolute value of the inner product should be less than the numerical precision threshold. The pattern partitioning results are stored in a separate data structure, including the dominant pattern basis vector matrix, the secondary pattern basis vector matrix, the corresponding singular value sequence, and the energy contribution rate sequence.

[0118] The projection calculation module extracts the row vectors of each monitoring node in the global topological feature matrix for pattern projection analysis. Row vector extraction obtains the 8-dimensional topological feature vector of the specified monitoring node through matrix indexing operations. Vector elements include the node's degree centrality, betweenness centrality, and edge weight distribution values. The dominant topological pattern projection calculation performs matrix multiplication between the node feature vector and the dominant pattern basis vector matrix to obtain the node's coordinate representation in the dominant pattern vector space. The modulus square calculation of the projection vector is obtained by summing the squares of each component of the projection vector; this value represents the energy distribution intensity of the node feature in the dominant topological pattern. The secondary topological pattern projection calculation uses the same matrix multiplication and modulus square calculation methods to project the node feature vector onto the secondary pattern vector space and calculate the corresponding energy distribution intensity.

[0119] The projection coefficient normalization module performs numerical processing and validity verification on the squared magnitude of the calculated projection vector. Projection coefficients in the dominant topological pattern are normalized by dividing the squared magnitude of the projection vector by the dimension of the dominant pattern vector space, eliminating the influence of dimensional differences on the projection coefficients. Projection coefficients in secondary topological patterns are processed using the same normalization method. The validity check of projection coefficients includes numerical range verification and physical meaning verification; projection coefficients should be non-negative real numbers within a reasonable range. Abnormal projection coefficients are handled using a truncation method, cutting off values ​​exceeding the 3σ range to boundary values. The projection calculation results are cached in an in-memory data structure, supporting fast query and update operations.

[0120] The energy distribution ratio calculation module performs the ratio calculation of the projection coefficients of the dominant topology mode and the projection coefficients of the secondary topology mode. Before the ratio calculation, a zero-value check is performed on the denominator. If the projection coefficient of the secondary topology mode is less than the numerical precision threshold, it is replaced with the smallest positive value to avoid division by zero errors. The energy distribution ratio is calculated by dividing the projection coefficient of the dominant topology mode by the projection coefficient of the secondary topology mode. The ratio reflects the relative importance of node characteristics in the dominant and secondary modes. The ratio result is logarithmically transformed by calculating the ratio's ln(), making the ratio distribution closer to a normal distribution and facilitating subsequent statistical analysis. The logarithmically transformed ratio is then standardized by calculating the mean and standard deviation of the logarithmic ratios of all nodes. The standardized energy distribution ratio is obtained by subtracting the mean from the logarithmic ratio of each node and then dividing by the standard deviation.

[0121] The data management and update module maintains the storage and version control of projection coefficients and energy distribution ratios. Data storage uses a relational database structure, with table structures including fields such as node identifier, calculation timestamp, dominant mode projection coefficient, secondary mode projection coefficient, original energy distribution ratio, logarithmic transformation ratio, and normalized ratio. The data update strategy employs incremental updates, recalculating projection coefficients and ratios only for nodes whose topological characteristics have changed. Historical data retention is set to a 3-month retention period; expired data is automatically archived to cold storage. The query interface provides flexible query functions based on node identifiers, time ranges, ratio thresholds, and other criteria.

[0122] In one optional implementation, determining the fault propagation path through reverse tracing and path analysis, based on the positional relationship of the candidate abnormal nodes in the dynamic association graph structure and the changing direction of the edge weights, includes:

[0123] Extract all incoming and outgoing edges of the candidate abnormal nodes in the candidate abnormal node set in the dynamic association graph structure, calculate the change magnitude of edge weight of each incoming edge between the current time window and the historical time window, determine the change direction of the edge weight according to the sign of the change magnitude of the edge weight, and select the incoming edge with the largest change magnitude of edge weight and the change direction of enhancement as the main import edge.

[0124] Using the source node of the main import edge as the current trace node, determine whether the current trace node belongs to the candidate abnormal node set. If it does, repeat the operation of extracting the main import edge and use the source node of the main import edge as the new current trace node. If it does not belong, mark the current trace node as the trace termination node.

[0125] Record all node sequences and connection edge sequences traversed during the reverse tracing process from the initial candidate abnormal node to the tracing termination node. Arrange the node sequences in reverse order of the reverse tracing to obtain the forward node sequence, and arrange the connection edge sequences in reverse order of the reverse tracing to obtain the forward connection edge sequence.

[0126] Based on the forward node sequence and the forward connection edge sequence, a directed path is constructed from the trace termination node to the initial candidate abnormal node, which serves as the fault propagation path.

[0127] In this specific embodiment, for each candidate abnormal node in the dynamic association graph structure, all its incoming and outgoing edge information is extracted. Incoming edges represent the connection of electrical signals flowing from other nodes to the candidate abnormal node, and outgoing edges represent the connection of the candidate abnormal node to other nodes. Taking the A section busbar in a 110kV substation as a candidate abnormal node as an example, this node is connected to 5 incoming edges, which are respectively from the adjacent B1 to B5 protection device monitoring nodes, and 3 outgoing edges are connected to the C1 to C3 load monitoring nodes.

[0128] Calculate the change in edge weight for each incoming edge between the current time window and historical time windows. For example, for the edge weight from B1 to A, extract the edge weight value of 0.83 in the current time window t and the edge weight value of 0.62 in the previous time window t-1, and calculate the change in edge weight as 0.21. Similarly, calculate the changes in edge weight from B2 to B5 as 0.15, -0.07, 0.03, and 0.08, respectively. Determine the direction of edge weight change based on the sign of the change in edge weight; a positive value indicates an increase in edge weight, and a negative value indicates a decrease in edge weight. In this example, the edge weights of incoming edges B1, B2, B4, and B5 change in an increasing direction, while the edge weight of incoming edge B3 changes in a decreasing direction.

[0129] The incoming edge with the largest change in edge weight and the direction of the change being strengthening is selected as the main incoming edge. In this example, the change in edge weight from B1 to A is 0.21, which is the largest change in edge weight and the direction of the change being strengthening among all incoming edges. Therefore, it is selected as the main incoming edge. This indicates that the electrical signal transmission strength from node B1 to node A is significantly increased, which is the main path for fault propagation.

[0130] Taking the source node B1 of the main import edge B1 to A as the current trace node, determine whether node B1 belongs to the candidate abnormal node set. Assuming that B1 is also a candidate abnormal node, continue to perform the main import edge extraction operation on B1. Node B1 has 3 incoming edges, which come from nodes D1, D2 and D3 respectively. The calculated edge weight change ranges of these three incoming edges are 0.32, 0.11 and -0.05 respectively. Among them, the edge weight change range of the incoming edge from D1 to B1 is the largest and the direction is increasing. Therefore, it is selected as the main import edge of B1.

[0131] Using D1 as the new current tracing node, we continue to determine whether D1 belongs to the set of candidate abnormal nodes. If D1 does not belong to the set of candidate abnormal nodes, the system marks D1 as the tracing termination node. This indicates that the fault may have originated from node D1 and propagated to node A through B1.

[0132] Record all node sequences and connection edge sequences traversed during the reverse tracing process from the initial candidate abnormal node A to the tracing termination node D1. The node sequence is [A, B1, D1], and the connection edge sequence is [B1 to A, D1 to B1]. Arrange the node sequence in reverse order of the reverse tracing to obtain the forward node sequence [D1, B1, A]; arrange the connection edge sequence in reverse order of the reverse tracing to obtain the forward connection edge sequence [D1 to B1, B1 to A].

[0133] Based on the forward node sequence and the forward connection edge sequence, a directed path is constructed from the tracing termination node D1 to the initial candidate abnormal node A, i.e., D1→B1→A, as the fault propagation path. This path indicates that the fault initially occurs at node D1, propagates to node B1, and eventually affects node A, leading to system abnormality.

[0134] In practical applications, it is necessary to analyze the specific numerical characteristics of the edge weights. For example, the edge weight from D1 to B1 was 0.45 before the fault occurred, and increased to 0.77 after the fault occurred; the edge weight from B1 to A was 0.62 before the fault occurred, and increased to 0.83 after the fault occurred. The significant increase in the edge weight values ​​indicates that the transmission strength of the electrical signal along the fault propagation path has increased significantly.

[0135] For complex network topologies, multiple candidate propagation paths may be identified. For example, if another candidate anomalous node E is connected to A, and its main import edge terminates at node F, another failure propagation path F→E→A will be obtained. The feasibility of each path is evaluated by comparing the change magnitude of the main import edge weights and the node anomalous contribution of different propagation paths. Paths with larger change magnitudes of edge weights and higher node anomalous contributions are considered more feasible failure propagation paths. For instance, the average change magnitude of edge weights in path D1→B1→A is 0.27, which is higher than the 0.19 in path F→E→A. Therefore, the system determines D1→B1→A as the primary failure propagation path.

[0136] Analyzing fault propagation paths allows for accurate location of fault sources and understanding of how faults propagate within the transmission network, providing valuable decision support information for equipment maintenance and fault handling. In actual transmission network maintenance, workers can directly inspect and repair the corresponding equipment based on the fault source location, significantly improving fault handling efficiency.

[0137] A second aspect of the present invention provides a fault diagnosis and life prediction system for power transmission network equipment, comprising:

[0138] The first unit is used to acquire time-series electrical characteristic data from multiple monitoring nodes of the power transmission network, wherein the time-series electrical characteristic data contains electrical physical quantities in multiple time-series dimensions;

[0139] The second unit is used to segment the time-series electrical characteristic data into time windows and extract statistical features of the electrical physical quantities within each time window to obtain a segmented statistical feature set.

[0140] The third unit is used to perform cross-port correlation analysis on the segmented statistical feature set based on the spatiotemporal correlation constraint model, establish a dynamic correlation graph structure between different monitoring nodes, and the edge weights of the dynamic correlation graph structure reflect the synchronicity and causal transmission of changes in electrical physical quantities between adjacent monitoring nodes, thereby obtaining spatiotemporal correlation features.

[0141] The fourth unit is used to calculate the abnormal contribution of each monitoring node based on the deviation between the local topological change and the global topological stability of each monitoring node in the spatiotemporal correlation features, and to mark the monitoring nodes whose abnormal contribution exceeds the preset judgment conditions as candidate abnormal nodes.

[0142] The fifth unit is used to determine the fault propagation path by reverse tracing and path analysis based on the positional relationship of the candidate abnormal nodes in the dynamic association graph structure and the changing direction of the edge weights.

[0143] The sixth unit is used to output the fault location result based on the starting node position of the fault propagation path.

[0144] A third aspect of the present invention provides an electronic device, comprising:

[0145] processor;

[0146] Memory used to store processor-executable instructions;

[0147] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0148] 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.

[0149] 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.

[0150] 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 therein. Such 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 fault diagnosis and lifespan prediction of power transmission network equipment, characterized in that, include: Acquire time-series electrical characteristic data from multiple monitoring nodes of the power transmission network, wherein the time-series electrical characteristic data contains electrical physical quantities in multiple time-series dimensions; The time-series electrical characteristic data is segmented into time windows, and statistical features are extracted from the electrical physical quantities within each time window to obtain a segmented statistical feature set. Based on the spatiotemporal correlation constraint model, cross-port correlation analysis is performed on the segmented statistical feature set to establish a dynamic correlation graph structure between different monitoring nodes. The edge weights of the dynamic correlation graph structure reflect the synchronicity and causal transmission of changes in electrical physical quantities between adjacent monitoring nodes, thus obtaining spatiotemporal correlation features. Based on the deviation between the local topological change and the global topological stability of each monitoring node in the spatiotemporal correlation features, the abnormal contribution of each monitoring node is calculated, and monitoring nodes whose abnormal contribution exceeds the preset judgment conditions are marked as candidate abnormal nodes. Based on the positional relationship of the candidate abnormal nodes in the dynamic association graph structure and the changing direction of the edge weights, the fault propagation path is determined through reverse tracing and path analysis; The fault location result is output based on the starting node position of the fault propagation path; Based on the spatiotemporal correlation constraint model, cross-port correlation analysis is performed on the segmented statistical feature set to establish a dynamic correlation graph structure between different monitoring nodes, resulting in spatiotemporal correlation features including: Calculate the similarity of the electrical physical quantity change trends of any two monitoring nodes in the segmented statistical feature set within the same time window, and the time series-based delayed correlation of electrical physical quantity changes between different monitoring nodes, to obtain the inter-node correlation metric; An initial association graph is constructed based on the inter-node association metric. Edges in the initial association graph that do not satisfy physical connectivity are removed according to the topological connection constraints of the power transmission network. At the same time, the retained edges are assigned directional attributes according to the transmission directionality of electrical physical quantities between adjacent monitoring nodes, thus obtaining a topologically constrained directed association graph. For each edge in the directed association graph after the topological constraints, the edge weight is calculated. The edge weight is obtained by weighted fusion of the synchronicity measure and the causal transitivity measure of the changes in electrical physical quantities between the two connected monitoring nodes. The directed association graph with the topological constraints is combined with the edge weights to form the dynamic association graph structure, and the degree centrality, betweenness centrality and edge weight distribution characteristics of each monitoring node in the dynamic association graph structure are extracted as the spatiotemporal association characteristics.

2. The method according to claim 1, characterized in that, The edge weight is obtained by weighted fusion of the synchronicity measure and the causal transitivity measure of the changes in electrical physical quantities between the two connected monitoring nodes, including: For two connected monitoring nodes, the waveforms of changes in electrical physical quantities of the two monitoring nodes within the same time window are extracted. The morphological similarity and the time difference of the peak arrival time between the waveforms of changes in electrical physical quantities are calculated. The synchronization metric is obtained based on the normalized reciprocal relationship between the morphological similarity and the time difference. For two connected monitoring nodes, the temporal sequence of the changes in electrical physical quantities between adjacent monitoring nodes is identified, and the transmission direction of the changes in electrical physical quantities is determined according to the temporal sequence. The response intensity ratio between the amplitudes of the changes in electrical physical quantities along the transmission direction of adjacent monitoring nodes is calculated, and the response intensity ratio is used as the causal transitivity metric. The synchronization metric and the causal transitivity metric are respectively assigned adaptive weighting coefficients, which are dynamically adjusted according to the electrical distance and load transmission capacity of the two connected monitoring nodes in the power grid topology. The edge weight is obtained by weighting and summing the synchronicity metric and the causal transitivity metric according to the adaptive weight coefficient.

3. The method according to claim 1, characterized in that, Based on the deviation between the local topological changes and the global topological stability of each monitoring node in the spatiotemporal correlation features, the anomaly contribution of each monitoring node is calculated, including: For each monitoring node in the spatiotemporal correlation features, a topological feature evolution trajectory of the monitoring node is constructed within multiple consecutive time windows. The topological feature evolution trajectory includes the degree centrality sequence, betweenness centrality sequence, and connection edge weight sequence of the monitoring node. The trajectory deviation distance of the topological feature evolution trajectory relative to the historical steady-state trajectory is calculated to obtain the local topological change amount. A global topology feature matrix is ​​constructed based on the dynamic correlation graph structure. Singular value decomposition is performed on the global topology feature matrix to obtain the dominant topology pattern and the secondary topology pattern. The energy distribution ratio between the projection coefficient of each monitoring node in the dominant topology pattern and the projection coefficient in the secondary topology pattern is calculated. When the energy distribution ratio is reversed, the degree of deviation of the monitoring node from the global topology stability is determined. The anomaly contribution of each monitoring node is obtained by coupling the trajectory deviation distance in the local topology change with the energy distribution ratio reversal amplitude in the deviation degree of the global topology stability.

4. The method according to claim 3, characterized in that, Singular value decomposition is performed on the global topology feature matrix to obtain the dominant topology pattern and the secondary topology pattern. The energy distribution ratio between the projection coefficients of each monitoring node in the dominant topology pattern and the projection coefficients in the secondary topology pattern is calculated, including: The global topological feature matrix is ​​constructed based on the degree centrality, betweenness centrality, and edge weight distribution characteristics of all monitoring nodes in the dynamic association graph structure. The rows of the global topological feature matrix correspond to the monitoring nodes, and the columns correspond to the topological feature dimensions. The global topological feature matrix is ​​decomposed by singular value to obtain the left singular vector matrix, the singular value diagonal matrix, and the right singular vector matrix. The singular values ​​in the singular value diagonal matrix are sorted according to their energy contribution rates. The left singular vectors corresponding to the first few singular values ​​whose cumulative energy contribution rates reach the preset energy concentration are selected as the dominant topology pattern. The left singular vectors whose energy contribution rates are after the preset energy concentration and whose singular values ​​are non-zero are selected as the secondary topology pattern. For each monitoring node, the row vector corresponding to the monitoring node in the global topology feature matrix is ​​extracted, and the row vector is projected into the vector spaces of the dominant topology pattern and the secondary topology pattern respectively. The square of the modulus of the projection vector of the row vector in the vector space of the dominant topology pattern is calculated as the projection coefficient in the dominant topology pattern, and the square of the modulus of the projection vector of the row vector in the vector space of the secondary topology pattern is calculated as the projection coefficient in the secondary topology pattern. The ratio of the projection coefficient in the dominant topology mode to the projection coefficient in the secondary topology mode is calculated to obtain the energy distribution ratio.

5. The method according to claim 1, characterized in that, Based on the positional relationships of the candidate anomaly nodes in the dynamic association graph structure and the changing direction of the edge weights, the fault propagation path is determined through reverse tracing and path analysis, including: Extract all incoming and outgoing edges of the candidate abnormal nodes in the candidate abnormal node set in the dynamic association graph structure, calculate the change magnitude of edge weight of each incoming edge between the current time window and the historical time window, determine the change direction of the edge weight according to the sign of the change magnitude of the edge weight, and select the incoming edge with the largest change magnitude of edge weight and the change direction of enhancement as the main import edge. Using the source node of the main import edge as the current trace node, determine whether the current trace node belongs to the candidate abnormal node set. If it does, repeat the operation of extracting the main import edge and use the source node of the main import edge as the new current trace node. If it does not belong, mark the current trace node as the trace termination node. Record all node sequences and connection edge sequences traversed during the reverse tracing process from the initial candidate abnormal node to the tracing termination node. Arrange the node sequences in reverse order of the reverse tracing to obtain the forward node sequence, and arrange the connection edge sequences in reverse order of the reverse tracing to obtain the forward connection edge sequence. Based on the forward node sequence and the forward connection edge sequence, a directed path is constructed from the trace termination node to the initial candidate abnormal node, which serves as the fault propagation path.

6. A fault diagnosis and life prediction system for power transmission network equipment, used to implement the method as described in any one of claims 1-5, characterized in that, include: The first unit is used to acquire time-series electrical characteristic data from multiple monitoring nodes of the power transmission network, wherein the time-series electrical characteristic data contains electrical physical quantities in multiple time-series dimensions; The second unit is used to segment the time-series electrical characteristic data into time windows and extract statistical features of the electrical physical quantities within each time window to obtain a segmented statistical feature set. The third unit is used to perform cross-port correlation analysis on the segmented statistical feature set based on the spatiotemporal correlation constraint model, establish a dynamic correlation graph structure between different monitoring nodes, and the edge weights of the dynamic correlation graph structure reflect the synchronicity and causal transmission of changes in electrical physical quantities between adjacent monitoring nodes, thereby obtaining spatiotemporal correlation features. The fourth unit is used to calculate the abnormal contribution of each monitoring node based on the deviation between the local topological change and the global topological stability of each monitoring node in the spatiotemporal correlation features, and to mark the monitoring nodes whose abnormal contribution exceeds the preset judgment conditions as candidate abnormal nodes. The fifth unit is used to determine the fault propagation path by reverse tracing and path analysis based on the positional relationship of the candidate abnormal nodes in the dynamic association graph structure and the changing direction of the edge weights. The sixth unit is used to output the fault location result based on the starting node position of the fault propagation path.

7. 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 5.

8. 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 5.

Images