Secondary equipment hidden danger mining method and system based on recording wave file and monitoring data

Abstract

The application relates to the technical field of data processing, and discloses a secondary equipment hidden danger mining method and system based on a recording wave file and monitoring data. The method comprises the following steps: constructing a correlation graph by analyzing the recording wave file, combining equipment topology and protection logic to calculate condition probability and intervention probability difference, identifying a cause-effect relationship, performing Fourier transform and wavelet decomposition on a sampling sequence to extract multi-scale features, injecting a small signal test excitation to estimate a transfer function through system identification, fusing monitoring data to form a health feature vector to determine a health state, inputting the correlation graph and the health state into a graph neural network to calculate a weighted attention coefficient, aggregating neighbor information to predict a fault propagation path, establishing a degradation model to correct a failure rate, predicting residual life to generate a graded early warning report. The application realizes the transformation from passive post-event analysis to active predictive maintenance, and improves the timeliness, accuracy and systematicness of secondary equipment hidden danger identification.

Description

Technical Field

[0001] This application relates to the field of data processing technology, and in particular to a method and system for identifying potential hazards in secondary equipment based on waveform recording files and monitoring data. Background Technology

[0002] Secondary equipment is a key component of power systems, responsible for measurement, protection, control, and monitoring. It includes protection devices, measurement and control devices, and fault recorders, and its operational status directly impacts the safety and stability of the power grid. Existing methods for identifying potential problems in secondary equipment primarily rely on post-hoc analysis of fault recording files. This involves parsing COMTRADE format recording files to extract voltage and current waveforms, using Fourier transform to calculate harmonic content, and combining this with equipment action reports to determine abnormal conditions. Some methods introduce wavelet transform for multi-scale decomposition of transient signals or employ neural networks for pattern recognition of historical fault samples to achieve automatic classification of typical fault types. Simultaneously, the system monitors performance parameters such as CPU utilization and memory usage of the collected equipment, triggering alarms by setting thresholds.

[0003] However, existing technologies have the following shortcomings: First, relying solely on post-failure analysis of waveform records is a passive response mode, only identifying problems after a fault occurs. It fails to detect potential issues during the incubation period of equipment degradation, leading to delayed fault detection and missed opportunities for optimal prevention. Second, existing methods analyze equipment as isolated entities, ignoring the complex interrelationships between devices in the secondary system. They only establish association rules based on statistical co-occurrence, failing to distinguish between spurious associations and genuine causal relationships, resulting in an inaccurate understanding of fault propagation mechanisms between devices. Third, traditional methods only extract single frequency or time domain features from waveform records, lacking a comprehensive characterization of the signal's time-frequency joint properties. They are insufficiently sensitive to weak signals indicating early-stage problems, making it difficult to capture subtle signs of equipment performance degradation.

[0004] Further analysis reveals that, due to the lack of an active detection mechanism, existing technologies cannot assess the health status trends of equipment under normal operating conditions. Even when performance monitoring data is collected, only simple threshold judgments are made, failing to integrate waveform characteristics with monitoring data to construct a comprehensive health assessment model. This results in incomplete health status assessments and a lack of predictive capabilities. A deeper problem is that even if an anomaly is identified in a single device, existing technologies cannot predict the propagation path and impact range of the fault in the secondary system network. This is because the lack of accurate modeling of causal relationships between devices and propagation analysis methods based on network topology makes it impossible to assess the possibility of a single-point fault evolving into a systemic risk. Furthermore, existing technologies lack lifetime prediction models based on physical degradation mechanisms, failing to quantify the remaining service life and failure probability of equipment. This leads to a lack of scientific basis for maintenance decisions, relying solely on experience or periodic inspections, resulting in wasted maintenance resources or overlooked failure risks. Summary of the Invention

[0005] This application provides a method and system for secondary equipment hidden danger mining based on waveform recording files and monitoring data. It solves the problems of delayed hidden danger discovery, inaccurate equipment correlation analysis, incomplete health status assessment, unpredictable fault propagation path, and difficulty in quantifying remaining life in the prior art by constructing an equipment association map based on causal reasoning, fusing waveform recording features and active detection health assessment model, and a fault propagation prediction method driven by graph neural network. It realizes the transformation from passive post-event analysis to proactive predictive maintenance, and improves the timeliness, accuracy, and systematicness of secondary equipment hidden danger identification.

[0006] Firstly, this application provides a method for identifying potential hazards in secondary equipment based on waveform recording files and monitoring data. The method includes:

[0007] Step S1: Parse the waveform recording file to obtain the sampling sequence, combine the equipment topology relationship and protection logic relationship, identify the causal relationship between the equipment by calculating the difference between the conditional probability and the intervention probability and quantify the causal strength, and construct a secondary equipment causal association map;

[0008] Step S2: Perform full-cycle Fourier transform on the sampled sequence to calculate the amplitude and phase of each harmonic, and perform wavelet decomposition to extract the energy distribution and statistics of each scale layer to obtain multi-scale waveform feature vectors.

[0009] Step S3: Inject a small signal test excitation into the secondary equipment and record the response output. Estimate the transfer function of the input-output relationship through the system identification method. Extract the frequency domain response characteristics and fuse them with the multi-scale waveform feature vector and equipment performance monitoring data to form a health feature vector. Calculate the distance of the health feature vector from the health benchmark to determine the health status of the equipment.

[0010] Step S4: Input the causal association graph of the secondary equipment and the health status into the graph neural network, calculate the weighted attention coefficient between nodes by combining the causal strength, update the node representation by aggregating neighbor node information, iteratively calculate the fault risk of each node based on the node representation and propagation strength, and identify the fault propagation path;

[0011] Step S5: Establish an equipment degradation model based on the historical evolution trajectory of the health feature vector, predict the remaining service life of the equipment according to the corrected equipment failure rate of the health status, and generate a graded early warning report in combination with the fault propagation path.

[0012] Secondly, this application provides a secondary equipment hidden danger detection system based on waveform recording files and monitoring data, the secondary equipment hidden danger detection system based on waveform recording files and monitoring data includes:

[0013] The parsing module is used to parse the waveform file to obtain the sampling sequence, combine the equipment topology relationship and protection logic relationship, identify the causal relationship between the equipment by calculating the difference between the conditional probability and the intervention probability, quantify the causal strength, and construct a secondary equipment causal association map;

[0014] The extraction module is used to perform full-cycle Fourier transform on the sampled sequence to calculate the amplitude and phase of each harmonic, and to perform wavelet decomposition to extract the energy distribution and statistics of each scale layer, thereby obtaining multi-scale waveform feature vectors.

[0015] The recording module is used to inject small signal test excitation into the secondary equipment and record the response output. It estimates the transfer function of the input-output relationship through the system identification method, extracts the frequency domain response characteristics and fuses them with the multi-scale waveform feature vector and equipment performance monitoring data to form a health feature vector, and calculates the distance of the health feature vector from the health benchmark to determine the health status of the equipment.

[0016] The input module is used to input the causal association graph of the secondary equipment and the health status into the graph neural network, calculate the weighted attention coefficient between nodes by combining the causal strength, update the node representation by aggregating neighbor node information, iteratively calculate the fault risk of each node based on the node representation and propagation strength, and identify the fault propagation path.

[0017] The prediction module is used to establish a device degradation model based on the historical evolution trajectory of the health feature vector, predict the remaining service life of the device according to the corrected device failure rate of the health status, and generate a graded early warning report in combination with the fault propagation path.

[0018] Thirdly, a secondary equipment hidden danger discovery device based on waveform recording files and monitoring data is provided, comprising: a memory and at least one processor, wherein the memory stores instructions; the at least one processor calls the instructions in the memory to cause the secondary equipment hidden danger discovery device based on waveform recording files and monitoring data to execute the above-mentioned secondary equipment hidden danger discovery method based on waveform recording files and monitoring data.

[0019] Fourthly, a computer-readable storage medium is provided, wherein instructions are stored in the computer-readable storage medium, which, when run on a computer, cause the computer to execute the above-described method for identifying potential hazards in secondary equipment based on waveform files and monitoring data.

[0020] The technical solution provided in this application obtains sampling sequences by parsing waveform recording files and combines them with equipment topology and protection logic relationships. It identifies causal relationships between devices and quantifies causal strength by calculating the difference between conditional probability and intervention probability, constructing a secondary device causal correlation map. This technical feature introduces a causal reasoning mechanism, distinct from existing statistical correlation analysis, accurately distinguishing between false correlations and true causal relationships, eliminating interference from confounding factors, and making the identification of device correlations more reliable. By performing a full-cycle Fourier transform on the sampling sequences to calculate the amplitude and phase of each harmonic, and performing wavelet decomposition to extract the energy distribution and statistics at each scale level, a multi-scale waveform feature vector is obtained. This technical feature combines frequency domain analysis with time-frequency domain analysis. Compared to existing single feature extraction methods, it can simultaneously capture the steady-state harmonic components and transient abrupt changes of the signal, exhibiting higher sensitivity to weak signals indicating early-stage potential problems. By injecting small-signal test stimuli into secondary equipment and recording the response output, the transfer function is estimated using a system identification method. Frequency domain response characteristics are extracted and fused with multi-scale waveform feature vectors and equipment performance monitoring data to form a health feature vector. The distance of the health feature vector from the health benchmark is calculated to determine the health status of the equipment. This technology realizes a mode shift from passively waiting for faults to actively detecting hidden dangers. By using small-signal excitation to obtain dynamic response characteristics without affecting the normal operation of the equipment, and by integrating multi-source heterogeneous data to construct a comprehensive health assessment model, the health status assessment becomes more comprehensive and predictive, enabling the early detection of hidden dangers during the incubation period of equipment deterioration.

[0021] This technology inputs the causal relationship map and health status of secondary equipment into a graph neural network, calculates weighted attention coefficients between nodes based on causal strength, updates node representations by aggregating neighbor node information, and iteratively calculates the failure risk of each node based on node representations and propagation strength to identify failure propagation paths. This technique uses graph neural networks to capture complex nonlinear interactions between devices. By integrating data-driven attention mechanisms with knowledge-driven causal strength, it makes failure propagation path prediction more accurate, identifies indirect propagation paths ignored by traditional methods, and assesses the possibility of single-point failures evolving into systemic risks. Furthermore, it establishes an equipment degradation model based on the historical evolution trajectory of health feature vectors, predicts the remaining service life of equipment based on the equipment failure rate adjusted by health status, and generates tiered early warning reports by combining these with failure propagation paths. This technique introduces a physical mechanism-driven life prediction method, using health status as a correction factor to adjust the failure rate, ensuring that the remaining life prediction considers the actual degree of equipment degradation. The prediction results are more consistent with the actual equipment condition. Simultaneously, by associating with failure propagation paths to generate hierarchical early warning information, it provides a scientific basis for maintenance decisions, realizing a shift from experience-based judgment to quantitative assessment, improving the utilization efficiency of maintenance resources and reducing failure risks. Attached Figure Description

[0022] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 This is a schematic diagram of an embodiment of the secondary equipment hidden danger discovery method based on waveform recording files and monitoring data in this application.

[0024] Figure 2 This is a schematic diagram illustrating the device degradation trend analysis based on the historical evolution trajectory of health feature vectors in an embodiment of this application.

[0025] Figure 3 This is a schematic diagram of an embodiment of the secondary equipment hidden danger detection system based on waveform files and monitoring data in this application.

[0026] Figure 4 This is a schematic block diagram of the secondary equipment hidden danger discovery device based on waveform recording files and monitoring data in an embodiment of the present invention. Detailed Implementation

[0027] This application provides a method and system for identifying potential hazards in secondary equipment based on waveform recordings and monitoring data. The terms "first," "second," "third," "fourth," etc. (if present)," in the specification, claims, and accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in a sequence other than that illustrated or described herein. Furthermore, the terms "comprising" or "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or devices.

[0028] For ease of understanding, the specific process of the embodiments of this application is described below. Please refer to [link / reference]. Figure 1 One embodiment of the secondary equipment hidden danger discovery method based on waveform recording files and monitoring data in this application includes:

[0029] Step S1: Parse the waveform recording file to obtain the sampling sequence, combine the equipment topology relationship and protection logic relationship, identify the causal relationship between the equipment by calculating the difference between the conditional probability and the intervention probability and quantify the causal strength, and construct a secondary equipment causal association map;

[0030] Step S2: Perform full-cycle Fourier transform on the sampled sequence to calculate the amplitude and phase of each harmonic, and perform wavelet decomposition to extract the energy distribution and statistics of each scale layer to obtain multi-scale waveform feature vectors.

[0031] Step S3: Inject a small signal test excitation into the secondary equipment and record the response output. Estimate the transfer function of the input-output relationship through the system identification method. Extract the frequency domain response characteristics and fuse them with the multi-scale waveform feature vector and equipment performance monitoring data to form a health feature vector. Calculate the distance of the health feature vector from the health benchmark to determine the health status of the equipment.

[0032] Step S4: Input the causal relationship graph and health status of the secondary equipment into the graph neural network, calculate the weighted attention coefficient between nodes by combining the causal strength, update the node representation by aggregating neighbor node information, iteratively calculate the fault risk of each node based on the node representation and propagation strength, and identify the fault propagation path;

[0033] Step S5: Establish an equipment degradation model based on the historical evolution trajectory of health feature vectors, predict the remaining service life of the equipment based on the corrected equipment failure rate according to the health status, and generate a graded early warning report in combination with the fault propagation path.

[0034] It is understood that the executing entity of this application can be a secondary equipment hidden danger detection system based on waveform files and monitoring data, or it can be a terminal or a server; the specific implementation is not limited here. This application's embodiment uses a server as an example for illustration.

[0035] Specifically, a framework for analyzing equipment degradation trends is constructed based on historical time-series data of health feature vectors. A physical mechanism-driven life prediction method is used to achieve quantitative assessment and graded early warning of potential hazards. First, an equipment degradation trend function is established based on the historical evolution trajectory of the health feature vectors. Specifically, the health feature vectors obtained in step three are arranged chronologically to form a time series. The time derivative of each dimension feature in this series is calculated, i.e., the difference in feature values ​​between adjacent moments divided by the time interval, to obtain the degradation rate vector representing the rate of decline in equipment performance. The Arrhenius-Weibull function is used to describe the change in equipment failure rate over time. This function includes two key parameters: a shape parameter β and a scale parameter η. The shape parameter determines the failure mode type, and the scale parameter represents the characteristic life. The maximum likelihood estimation method is used to statistically analyze the failure times and censored data in the historical operating data. After constructing the likelihood function, the Newton iteration method is used to solve for the parameter estimates that maximize the likelihood function, thus obtaining the functional expression describing the equipment degradation trend.

[0036] Next, the equipment failure rate is adjusted based on the determined equipment health status. This process introduces health as a correction factor. The Mahalanobis distance between the current health feature vector and the baseline health mean vector is denoted as D². The difference between this distance and the critical failure threshold is normalized to obtain the current health value. The closer the health value is to 1, the healthier the equipment; the closer it is to 0, the closer the equipment is to the edge of failure. A proportional hazards function is used to describe the impact mechanism of health on the failure rate. The baseline failure rate is multiplied by a natural exponential function with the health correction term as the exponent. The correction term is calculated by subtracting the current health value from 1 and then multiplying by the health influence coefficient. The health influence coefficient is estimated from historical failure cases using a Cox regression model. This coefficient reflects the accelerating effect of deteriorating equipment health on failure risk. The adjusted actual failure rate can truly reflect the current failure tendency of the equipment.

[0037] Then, the remaining useful life of the equipment is predicted based on the corrected equipment failure rate. The prediction process involves substituting the corrected actual failure rate into the exponential integral expression of the reliability function, which describes the probability that the equipment will still be able to operate normally after a given time. The remaining useful life in the mathematical expectation sense is obtained by performing time-domain integration of the reliability function from the current time to infinity. Considering the uncertainty of parameter estimation, the Monte Carlo simulation method is used to evaluate the probability distribution of the remaining useful life. Specifically, random sampling is performed within the confidence intervals of the shape and scale parameters. After obtaining a set of parameter values ​​for each sampling, the corresponding remaining useful life is calculated. After repeating the sampling 10,000 times, all calculation results are statistically analyzed. The 50th percentile is extracted as the point estimate of the remaining useful life, and the 90th percentile is extracted as the upper confidence bound. The difference between the two reflects the degree of uncertainty in the prediction.

[0038] Finally, the predicted remaining service life of the equipment is correlated with the fault propagation path identified in step four to generate a graded early warning report. The correlation process involves traversing each equipment node in the fault propagation path, querying the predicted remaining service life of the corresponding equipment, and classifying the equipment into three risk levels according to the remaining service life value: equipment with a remaining service life of less than 30 days is marked as red (severe hazard level); equipment with a remaining service life between 30 and 60 days is marked as orange (moderately abnormal level); and equipment with a remaining service life greater than 60 days but a health level below 0.8 is marked as yellow (slightly deteriorated level). For cases where multiple high-risk equipment exist in the fault propagation path, a comprehensive risk index for the entire substation is calculated by weighted summation of network centrality weight and fault consequence severity. This index comprehensively reflects the overall safety level of the substation. The graded early warning report includes equipment-level early warning information listing the current health level, deterioration trend, predicted remaining service life, and hazard category of each piece of equipment; path-level propagation information plotting the fault propagation path and marking the trigger probability and expected propagation delay of each node; and a system-level comprehensive risk index providing the overall risk score and risk level determination results for the entire substation.

[0039] In one specific embodiment, step S1 includes:

[0040] According to the configuration file of the standard analytical waveform recording file, obtain the number of analog channels, the number of digital channels, and the sampling frequency. Read the time series sampling point data of each channel from the data file according to the file format type to form a sampling sequence.

[0041] For any two device nodes, the conditional probability is calculated by statistically analyzing the co-occurrence frequency of device anomalies from historical operating data. The intervention probability is calculated by observing changes in device status while controlling other variables. When the difference between the conditional probability and the intervention probability exceeds a set threshold, an instrumental variable is introduced. The causal effect coefficient is estimated as the causal strength using the two-stage least squares method.

[0042] A directed graph is constructed as a causal relationship graph of secondary devices, where nodes represent secondary devices, directed edges represent causal relationships, and edge weights are obtained by weighted summation of statistical co-occurrence probability, causal strength, and the satisfaction of protection logic relationship.

[0043] Specifically, the sampling sequence is obtained by parsing the waveform recording file according to the IEEE COMTRADE standard. The waveform recording file consists of two parts: a configuration file and a data file. The configuration file uses the .cfg extension to store channel configuration information, and the data file uses the .dat extension to store the actual sampled data. When parsing the configuration file, the text content is read line by line. The first line obtains the substation name and equipment number; the second line obtains the total number of channels and channel type identifiers, where the channel type identifiers include the number of analog channels and the number of digital channels; the third line obtains the sampling frequency value; and subsequent lines read the name, unit, transformation ratio, and offset parameters of each channel. Based on the sampling frequency and channel number information obtained from the configuration file, the sampling point data is read from the data file in either ASCII or binary format. In ASCII format, each line of data is split using commas to obtain the timestamp and the value of each channel. In binary format, the values ​​at the corresponding positions are read according to the byte offset and data type definition. The read raw sampling point data is then linearly transformed into engineering quantity values ​​according to the transformation ratio and offset parameters in the configuration file and arranged in chronological order to form a sampling sequence.

[0044] Next, we identify causal relationships and quantify the strength of causality between any two device nodes. This process differs from traditional statistical association analysis by introducing a causal inference mechanism to determine the true causal effect between devices. First, we calculate the conditional probability by statistically analyzing the co-occurrence frequency of device anomalies from historical operational data. We iterate through all historical records to find the sets of times when device A and device B are abnormal, and count the number of times both device A and device B are abnormal. Dividing the co-occurrence count by the total number of times device A is abnormal yields the conditional probability of device B being abnormal under the condition that device A is abnormal. Then, we calculate the intervention probability by observing changes in device status while controlling for other variables. The intervention probability represents the probability that device B will become abnormal after forcibly setting device A to an abnormal state. This probability needs to exclude the interference of confounding factors, which are third-party factors that simultaneously affect both device A and device B. We calculate the difference between the conditional probability and the intervention probability. When the absolute value of the difference exceeds a set threshold, it indicates the presence of a significant confounding factor leading to a spurious association. At this point, we introduce an instrumental variable to eliminate the influence of the confounding factor. The instrumental variable needs to meet three conditions: it is related to device A, it only affects device B through device A, and it is not related to the confounding factor.

[0045] A two-stage least squares method is used to estimate the causal effect coefficient as the causal strength. This method involves parameter estimation in two stages. In the first stage, a regression equation is established with the instrumental variable as the independent variable and the state of device A as the dependent variable. The regression coefficients are estimated using the least squares method to obtain the predicted value of the state of device A. This predicted value includes the influence of the instrumental variable on device A but excludes the influence of confounding factors. In the second stage, the predicted value of the state of device A is used as the independent variable and the state of device B as the dependent variable to establish a regression equation. The regression coefficients are estimated using the least squares method to obtain the causal effect coefficient, which reflects the strength of the pure causal effect of device A on device B. False causal relationships are filtered out by incorporating the physical constraints of the power system. These physical constraints include power flow equation constraints, protection logic criterion constraints, and secondary loop connection constraints. The method checks whether the identified causal relationships violate the voltage balance law or the current conservation law, and checks whether there are physical connection paths between devices. Causal relationships that violate physical laws are eliminated.

[0046] Finally, a directed graph is constructed as a causal relationship graph for secondary equipment. Each node in the graph represents a secondary device, and the node attributes include static information such as device type, device number, and commissioning time. Directed edges represent the causal relationship between devices, pointing from the cause node to the effect node, and the direction of the edge indicates the propagation direction of the causal effect. The edge weight is obtained by weighted summation of three parts. The first part is the statistical co-occurrence probability, which represents the historical frequency of the simultaneous occurrence of anomalies in two devices. The normalized co-occurrence probability value is obtained by dividing the number of co-occurrences by the total number of anomaly events. The second part is the causal strength, which is the causal effect coefficient estimated by the two-stage least squares method mentioned above. This coefficient quantifies the magnitude of the direct causal effect of the cause device on the effect device. The third part is the satisfaction of the protection logic relationship. According to the logic configuration scheme of the protection device, it is checked whether there is linkage logic or interlocking logic between the two devices. The satisfaction value is 1 when linkage logic exists and 0 when it does not exist. The final edge weight value is obtained by multiplying the three parts by their respective weighting coefficients and then summing them. The sum of the weighting coefficients equals 1 to ensure weight normalization. The specific value of the weighting coefficients is adjusted according to the actual application scenario. In general, the weighting coefficient of causal strength is set to the maximum value to highlight the role of causal reasoning.

[0047] In one specific embodiment, step S2 includes:

[0048] The full-cycle Fourier algorithm is applied to the sampled sequence to calculate the sine and cosine coefficients of the sampled points within one power frequency cycle. The amplitude and phase of each harmonic are then calculated based on the sine and cosine coefficients.

[0049] The sampling sequence is decomposed into multiple wavelet functions by selecting the mother wavelet function, and the scale parameter is discretized into different decomposition layers. Each decomposition layer corresponds to a specific frequency band range, and the wavelet coefficients of each decomposition layer are calculated.

[0050] For the wavelet coefficients of each decomposition layer, the energy value is calculated as the energy distribution, and the mean, standard deviation, skewness and kurtosis are calculated as statistics.

[0051] The amplitude, phase, energy distribution, and statistics of each harmonic are spliced ​​together to form a multi-scale waveform feature vector.

[0052] Specifically, the full-cycle Fourier transform algorithm is applied to the sampled sequence to extract frequency domain features. This algorithm, specifically designed for power system frequency signals, differs from the traditional discrete Fourier transform by utilizing sampled data from a complete power frequency cycle to eliminate spectral leakage caused by non-integer period truncation. For a 50Hz power frequency system with a sampling frequency of 1200Hz, one power frequency cycle contains 24 sampling points. The sampled sequence is segmented according to the power frequency cycle, and 24 consecutive sampling points are extracted from each segment to form an analysis window. For each analysis window, the sine and cosine coefficients of each harmonic are calculated. The sine coefficients are calculated by iterating through all sampling points within the window, multiplying the value of each sampling point by the corresponding sine function value at that time. The angular frequency of the sine function is the harmonic order multiplied by the power frequency angular frequency. All products are summed and then multiplied by a normalization coefficient to obtain the sine coefficients. The cosine coefficient is calculated by iterating through all sampling points within the window, multiplying the value of each sampling point by the cosine function value at the corresponding time. The angular frequency of the cosine function is also calculated by multiplying the harmonic order by the power frequency angular frequency. All products are summed, and after special processing of the first and last sampling points, the result is multiplied by a normalization coefficient to obtain the cosine coefficient. The amplitude and phase of each harmonic are then calculated based on the sine and cosine coefficients. The amplitude is calculated by adding the squares of the sine and cosine coefficients, taking the square root, and then dividing by the square root of 2 for effective value conversion. The phase is calculated by performing an arctangent operation on the ratio of the cosine and sine coefficients to obtain the initial phase angle.

[0053] The sampled sequence is decomposed into time-frequency domain features using a multi-level wavelet decomposition process, employing a wavelet mother function. The Daubechies4 function is chosen as the mother function due to its compact support characteristics and good frequency selectivity, making it suitable for analyzing transient signals in power systems. The wavelet decomposition process divides the signal into detail and approximation components in different frequency bands. Each decomposition level corresponds to a specific frequency range, determined by both the sampling frequency and the number of decomposition levels. The scale parameter is discretized into different decomposition levels, increasing in powers of 2: 2 for the first level, 4 for the second, 8 for the third, and so on. The maximum number of decomposition levels is set to the binary logarithm of the total number of sampling points; for a sequence with 1024 sampling points, this equates to 10 binary logarithms. The frequency band range corresponding to each decomposition layer is calculated by dividing the sampling frequency by a power of the scale parameter. The frequency band range corresponding to the first layer is from 1200Hz (2^1) to 600Hz (2^2), which is 300Hz. The frequency band range corresponding to the second layer is 300Hz to 150Hz, and the frequency band range corresponding to the third layer is 150Hz to 75Hz. The process of calculating the wavelet coefficients of each decomposition layer involves performing a convolution operation on the sampled sequence using a wavelet mother function. In the convolution operation, the wavelet mother function is translated on the time axis and scaled at different scales. The wavelet coefficients at that position and scale are obtained by calculating the inner product of the signal and the scaled / translated wavelet function.

[0054] Then, energy values ​​are calculated for the wavelet coefficients of each decomposition layer as the energy distribution. The energy value reflects the intensity of the signal component in that frequency band. The calculation method is to iterate through all wavelet coefficients of the decomposition layer, square each wavelet coefficient, and sum them to obtain the total energy value of the layer. The mean, standard deviation, skewness, and kurtosis of the wavelet coefficients of each decomposition layer are calculated as statistical quantities to describe the statistical characteristics of the signal. The mean is calculated by summing all wavelet coefficients of the layer and dividing by the number of coefficients. The standard deviation is calculated by squarening the difference between each wavelet coefficient and the mean, summing the results, dividing by the number of coefficients, and then taking the square root. The skewness is calculated by dividing the difference between each wavelet coefficient and the mean by the standard deviation, cubed the result, summing all the cubed values, and then dividing by the number of coefficients. The skewness value reflects the symmetry of the wavelet coefficient distribution; a positive skewness indicates a right-skewed distribution, and a negative skewness indicates a left-skewed distribution. The kurtosis calculation method is to divide the difference between each wavelet coefficient and the mean by the standard deviation and then take the fourth power. Sum all the fourth power values ​​and then divide by the number of coefficients. The kurtosis value reflects the sharpness of the wavelet coefficient distribution. A kurtosis greater than 3 indicates that the distribution is sharper than the normal distribution, while a kurtosis less than 3 indicates that the distribution is flatter than the normal distribution.

[0055] Finally, the amplitude and phase, energy distribution and statistics of each harmonic are concatenated to form a multi-scale waveform feature vector. The concatenation process is arranged in order of feature type. First, the amplitude values ​​of the fundamental wave to the 10th harmonic are arranged in a total of 10 values. Then, the phase values ​​of the fundamental wave to the 10th harmonic are arranged in a total of 10 values. Next, the energy distribution values ​​of the 1st to the 10th layers are arranged in a total of 10 values. Then, the mean, standard deviation, skewness and kurtosis values ​​of the 1st to the 10th layers are arranged in order. All values ​​are concatenated column by column to form a one-dimensional feature vector. The total dimension of this vector is 10 + 10 + 10 + 10 + 10 + 10 + 10 = 70 dimensions. Each dimension of the feature vector corresponds to a specific frequency domain or time-frequency domain characteristic of the sampling sequence. The fundamental amplitude reflects the effective value level of voltage or current, the amplitude of each harmonic reflects the degree of high-frequency pollution, the low-frequency layer energy reflects the intensity of the power frequency component, the high-frequency layer energy reflects the intensity of the transient impulse component, the mean of the wavelet coefficients reflects the DC bias of the signal in that frequency band, the standard deviation reflects the fluctuation amplitude of the signal in that frequency band, the skewness reflects the waveform symmetry of the signal in that frequency band, and the kurtosis reflects the pulse characteristics of the signal in that frequency band.

[0056] In one specific embodiment, a wavelet mother function is selected to perform multi-level wavelet decomposition on the sampled sequence, discretizing the scale parameter into different decomposition layers. Each decomposition layer corresponds to a specific frequency band range, and the wavelet coefficients of each decomposition layer are calculated, including:

[0057] Select the wavelet mother function to perform continuous wavelet transform on the time series signal in the sampled sequence, where the scale parameter controls the frequency resolution and the translation parameter controls the time positioning.

[0058] The scale parameter is discretized according to a power of 2. The maximum number of decomposition layers is set to the binary logarithm of the total number of sampling points. The frequency range corresponding to each layer is calculated by dividing the sampling frequency by a power of the scale parameter.

[0059] Wavelet mother function is applied to each decomposition layer for convolution operation, and wavelet transform coefficients of the sampling sequence at each time position and each scale layer are calculated as wavelet coefficients.

[0060] The amplitude sequences of wavelet coefficients from each decomposition layer are extracted for subsequent energy distribution and statistical calculations.

[0061] Specifically, the Daubechies4 wavelet function is chosen as the mother wavelet function to perform continuous wavelet transform on the time series signal in the sampled sequence. The mother wavelet function is a basis function with oscillatory and rapidly decaying characteristics. The Daubechies4 function has four vanishing moments and a tight support property. The vanishing moment represents the order of the wavelet function's orthogonality to the polynomial; a higher vanishing moment indicates a stronger ability to extract signal details. The tight support property indicates that the wavelet function takes a value of zero outside a finite interval, which reduces computational complexity. The continuous wavelet transform includes two variables: a scaling parameter and a translation parameter. The scaling parameter controls the degree of scaling of the wavelet function, thereby controlling the frequency resolution. A larger scaling parameter results in a wider wavelet function with lower frequency resolution but higher time resolution, while a smaller scaling parameter results in a narrower wavelet function with higher frequency resolution but lower time resolution. The translation parameter controls the position of the wavelet function on the time axis, thereby controlling the time positioning. The translation parameter slides from the start time to the end time of the sampling sequence. At each translation position, the inner product of the signal and the wavelet function is calculated to obtain the wavelet transform value at that time and scale.

[0062] Next, the scale parameters are discretized according to powers of 2 to convert the continuous wavelet transform into a discrete wavelet transform. During discretization, the scale parameter values ​​are limited to integer powers of 2. The scale parameter for the first layer is 2 to the power of 1 (2), the scale parameter for the second layer is 2 to the power of 2 (4), the scale parameter for the third layer is 2 to the power of 3 (8), and so on. The maximum number of decomposition layers is set to the binary logarithm of the total number of sampling points. The binary logarithm is calculated by taking the base-2 logarithm of the total number of sampling points. When the total number of sampling points is 512, the binary logarithm is 9 layers; when the total number of sampling points is 1024, the binary logarithm is 10 layers; and when the total number of sampling points is 2048, the binary logarithm is 11 layers. The frequency range corresponding to each layer is calculated by dividing the sampling frequency by a power of the scale parameter. The upper limit frequency of the frequency range is the sampling frequency divided by the scale parameter of that layer, and the lower limit frequency is the sampling frequency divided by the scale parameter of the next layer. The upper limit of the frequency range of the first layer is the sampling frequency 1200Hz divided by 2, which equals 600Hz, and the lower limit is 1200Hz divided by 4, which equals 300Hz. The first layer corresponds to the frequency range of 300Hz to 600Hz. The upper limit of the frequency range of the second layer is 300Hz, and the lower limit is 1200Hz divided by 8, which equals 150Hz. The second layer corresponds to the frequency range of 150Hz to 300Hz.

[0063] Then, wavelet coefficients are calculated by applying a wavelet mother function to each decomposition layer through convolution. The convolution process involves scaling the wavelet mother function according to the scale parameter of that layer. The width of the scaled wavelet function is equal to the original width multiplied by the scale parameter. The scaled wavelet function is then slid along the time axis, and at each time point, the value of the sampled sequence is multiplied by the value of the wavelet function at the corresponding position, and the sum is obtained to obtain the convolution result at that position. The wavelet transform coefficients of the sampled sequence at each time point and at each scale layer are calculated as wavelet coefficients. The wavelet coefficients are a two-dimensional array, with the first dimension corresponding to the time position index and the second dimension corresponding to the scale layer index. Each wavelet coefficient value reflects the signal strength of the sampled sequence within a specific frequency band at a specific time point. The wavelet coefficients of the first layer reflect the time-varying signal components in the 300Hz to 600Hz frequency band, the wavelet coefficients of the second layer reflect the time-varying signal components in the 150Hz to 300Hz frequency band, the wavelet coefficients of the low-frequency layer reflect the power frequency and its low-order harmonic components, and the wavelet coefficients of the high-frequency layer reflect transient impulses and high-frequency noise components.

[0064] Finally, the amplitude sequences of wavelet coefficients from each decomposition layer are extracted for subsequent calculations. The method for extracting the amplitude sequences involves iterating through the wavelet coefficients at all time positions within the decomposition layer, taking the absolute value of each wavelet coefficient to obtain its amplitude, and arranging the amplitudes at all time positions in chronological order to form the amplitude sequence for that layer. The amplitude sequence retains the intensity information of the wavelet coefficients but discards the sign information. The length of the amplitude sequence is equal to the number of time sampling points of the wavelet coefficients at that layer. The length of the amplitude sequence varies between different decomposition layers; the higher the decomposition layer, the shorter the amplitude sequence. The length of the amplitude sequence for the first layer is approximately half the length of the original sampling sequence, and the length of the amplitude sequence for the second layer is approximately one-quarter of the length of the original sampling sequence. The amplitude sequences of each decomposition layer serve as input data for subsequent energy distribution and statistical calculations. Energy distribution calculations require summing the squares of each element in the amplitude sequence, while statistical calculations require calculating the mean, standard deviation, skewness, and kurtosis of the amplitude sequence.

[0065] In one specific embodiment, step S3 includes:

[0066] A sine wave with an amplitude of 0.1%-1% of the rated operating level is generated as a small signal test excitation. The small signal test excitation is injected into the secondary equipment, and the equipment output signal is recorded by the data acquisition device as the response output.

[0067] An autoregressive moving average model is established to describe the input-output relationship. The recursive least squares algorithm is used to identify parameters of small signal test excitation and response output, estimate the polynomial coefficients of the model, and calculate the transfer function based on the polynomial coefficients.

[0068] The transfer function is transformed in the frequency domain, and the amplitude frequency response and phase frequency response are calculated at the preset characteristic frequency points. The amplitude frequency response value and phase frequency response value at each characteristic frequency point are combined to form a frequency domain eigenvector as the frequency domain response characteristic.

[0069] Frequency domain response characteristics, multi-scale waveform feature vectors, and equipment performance monitoring data are concatenated column by column to form a health feature vector. Multiple sets of health feature vectors from healthy equipment are collected to calculate the mean vector and covariance matrix as a health benchmark model. The Mahalanobis distance between the health feature vector and the mean vector of the equipment to be evaluated is calculated. When the Mahalanobis distance exceeds the detection threshold set based on the chi-square distribution, it is determined to be an abnormal state as the equipment health status.

[0070] Specifically, a sine wave with an amplitude of 0.1% to 1% of the rated operating level is generated as a small-signal test excitation. The rated operating level refers to the nominal input voltage or current level of the secondary equipment during normal operation. For voltage-type protection devices, the rated operating level is 100V, and for current-type protection devices, it is 5A. The amplitude is set within the range of 0.1% to 1% of the rated operating level. If the value is too small, the signal-to-noise ratio will be insufficient, leading to increased measurement errors; if the value is too large, it will interfere with normal operation or even trigger protection actions. The angular frequency of the sine wave is selected to be the power frequency or an integer multiple thereof. The power frequency is twice pi multiplied by 50Hz, which equals 314 radians per second. Integer multiples include twice the power frequency, three times the power frequency, up to ten times the power frequency, corresponding to the 100Hz to 500Hz frequency band. When generating the sine wave, the initial phase angle is randomly selected within the range of 0 to twice pi. Randomizing the initial phase angle avoids synchronization lock between the test excitation and the internal periodic signals of the equipment. Small-signal test stimuli are injected into the secondary equipment via a dedicated test interface or auxiliary input channel. The injection process is performed under normal equipment operation, and the injection duration is set to at least 10 power frequency cycles to ensure the statistical stability of the signal. The equipment output signal is recorded as the response output by a data acquisition device. The sampling frequency of the data acquisition device is set to at least 10kHz to meet the requirements for capturing high-frequency components. The number of acquisition channels corresponds to the number of equipment output ports, and all output signals are acquired synchronously to avoid phase differences.

[0071] Next, an autoregressive moving average (ARMAX) model is established to describe the input-output relationship. This model consists of three parts: an autoregressive term, a moving average term, and an external input term. The autoregressive term describes the dependence of the output signal on its historical values, the moving average term describes the dependence of the output signal on historical noise values, and the external input term describes the driving effect of the input signal on the output signal. A recursive least squares (RLS) algorithm is used to identify parameters for the small-signal test excitation and response output. RLS is an online parameter estimation method that updates the parameter estimates each time a new input-output data pair is received, unlike batch least squares which requires accumulating all data before calculation. The recursive process includes three steps: prediction error calculation, gain matrix update, and parameter vector update. The prediction error is the difference between the actual output and the model's predicted output. The gain matrix reflects the contribution weight of new data to the parameter update and is obtained by matrix multiplication and scalar division of the covariance matrix and the regression vector. The parameter vector update formula is to add the product of the gain matrix and the prediction error to the old parameter vector to obtain the new parameter vector. The covariance matrix update formula introduces a forgetting factor to adjust the weights of the old and new data. The forgetting factor ranges from 0.95 to 1; the smaller the value, the faster the forgetting of historical data. The polynomial coefficients of the estimation model include the coefficients of the autoregressive polynomial and the coefficients of the external input polynomial. The order of the autoregressive polynomial indicates how many historical time points the output signal depends on, and the order of the external input polynomial indicates how many time points the influence of the input signal is delayed in reaching the output. The transfer function is calculated based on the polynomial coefficients. The transfer function is the ratio of the external input polynomial to the autoregressive polynomial. The numerator polynomial coefficients come from the estimation results of the external input term, and the denominator polynomial coefficients come from the estimation results of the autoregressive term with a 1 added to the constant term.

[0072] Then, a frequency domain transformation is performed on the transfer function to obtain its frequency domain response characteristics. This transformation converts the time-domain transfer function into a frequency-domain complex function by replacing the delay operator in the transfer function with a complex exponential function. The amplitude-frequency response and phase-frequency response are calculated at preset characteristic frequency points. These characteristic frequency points are selected from 50Hz (power frequency) and its 2nd to 10th harmonics, corresponding to frequencies from 100Hz to 500Hz (10 frequency points in total). The amplitude-frequency response is calculated by taking the modulus of the complex value of the frequency-domain complex function at that frequency point. The modulus reflects the amplification or attenuation factor of the input signal at that frequency; a modulus greater than 1 indicates amplification, and a modulus less than 1 indicates attenuation. The phase-frequency response is calculated by taking the argument of the complex value of the frequency-domain complex function at that frequency point. The argument reflects the phase lead or lag of the output signal relative to the input signal; a positive argument indicates lead, and a negative argument indicates lag. The amplitude and phase response values ​​at each characteristic frequency point are used to form a frequency domain feature vector as the frequency domain response characteristic. The frequency domain feature vector contains 20 elements, including 10 amplitude response values ​​and 10 phase response values, arranged from low to high frequency.

[0073] Finally, the frequency domain response characteristics, multi-scale waveform feature vectors, and equipment performance monitoring data are concatenated column-wise to form a health feature vector. The multi-scale waveform feature vectors are derived from the Fourier transform and wavelet decomposition results in step two. The equipment performance monitoring data includes four indicators: CPU utilization, memory usage, network traffic, and response time. Column-wise concatenation refers to joining the three vectors end-to-end into a single long vector. The total dimension of the concatenated health feature vector equals the frequency domain response characteristic dimension plus the multi-scale waveform feature vector dimension plus the equipment performance monitoring data dimension. Multiple sets of health feature vectors from healthy devices are collected to calculate the mean vector and covariance matrix as a health baseline model. Healthy devices refer to newly commissioned equipment that has passed factory testing. The number of collection groups is set to more than 50 to ensure statistical reliability. Each set of health feature vectors corresponds to one test stimulus injection and response collection process. The mean vector is calculated by summing all health feature vectors by dimension and then dividing by the number of collection groups. The covariance matrix is ​​calculated by subtracting the mean vector from each set of health feature vectors, multiplying it by its transpose, and then summing all the outer product matrices and dividing by the number of collection groups minus 1. The Mahalanobis distance between the health feature vector and the mean vector of the device to be evaluated is calculated. The Mahalanobis distance calculation process involves subtracting the mean vector from the health feature vector to obtain the deviation vector. This deviation vector is then multiplied on the left by the inverse of the covariance matrix, and then on the right by the transpose of the deviation vector, performing a quadratic operation. The result of this quadratic operation is the squared value of the Mahalanobis distance. Mahalanobis distance considers the correlation and variance differences between features of different dimensions, making it more suitable for anomaly detection in multidimensional correlated data compared to Euclidean distance. When the Mahalanobis distance exceeds a detection threshold set based on a chi-square distribution, the device is considered to be in an abnormal state and its health status is determined. The degrees of freedom of the chi-square distribution are equal to the dimensions of the health feature vector. The detection threshold is the quantile of the chi-square distribution at a significance level of 0.01. A significance level of 0.01 means that the probability of a normal device being misjudged as abnormal is controlled within 1%.

[0074] In one specific embodiment, step S4 includes:

[0075] The initial features of the nodes in the causal relationship graph of the secondary device are set as the health feature vectors of the corresponding devices. The initial features of each node are linearly transformed by the shared weight matrix. The attention function is applied to calculate the original attention score of the node to its neighboring nodes. The original attention score is normalized to obtain the attention weight. The attention weight is multiplied by the causal strength of the corresponding edge and normalized to obtain the weighted attention coefficient.

[0076] The transformation features of neighboring nodes are weighted and summed based on the weighted attention coefficients. After processing by the activation function, the current node features are updated. The neighbor information aggregation and node feature update process is repeated by stacking multiple graph convolutional layers to obtain the final node representation.

[0077] The propagation strength of the source node is initialized to 1 and the propagation strength of the other nodes is 0. The propagation strength vector is propagated iteratively using a normalized adjacency matrix. A restart probability is introduced to prevent over-propagation. After a preset number of iterations, the converged propagation strength is obtained.

[0078] The node representation, propagation intensity, and prediction time window are concatenated and then the fault risk of each node is calculated through a fully connected layer. The node sequence with fault risk exceeding a set probability threshold is extracted as the fault propagation path.

[0079] Specifically, the initial features of the nodes in the causal correlation graph of secondary devices are set as the health feature vectors of the corresponding devices. The causal correlation graph contains nodes and directed edges based on the construction results. Each node corresponds to a secondary device, and the initial features of the node are the health feature vectors calculated in step three. The health feature vectors include frequency domain response characteristics, multi-scale waveform characteristics, and device performance monitoring data. A linear transformation is performed on the initial features of each node through a shared weight matrix. The shared weight matrix is ​​a two-dimensional matrix, and all nodes use the same weight matrix for feature transformation. The linear transformation process involves multiplying the initial feature vector of the node with the weight matrix. Matrix multiplication is calculated according to the rule of summing corresponding element-wise multiplication of rows and columns. The dimension of the transformed features is determined by the number of columns in the weight matrix. An attention function is applied to calculate the original attention score of a node to its neighboring nodes. A neighboring node is defined as a node in the causal relationship graph that has a directed edge connection to the current node. The attention function is implemented using a single-layer feedforward neural network, concatenating the transformed features of the current node with those of its neighboring nodes. The concatenated vector is then used for an inner product operation with a learnable parameter vector. The inner product result is processed by the LeakyReLU activation function to obtain the original attention score. The LeakyReLU function remains unchanged for positive numbers and multiplies negative numbers by a slope coefficient less than 1. The original attention scores are then normalized to obtain the attention weights. Normalization is implemented using the softmax function, which takes the natural exponent of the original attention score of the current node to a neighboring node and divides it by the sum of the natural exponents of the original attention scores of the current node to all its neighboring nodes. The sum of the normalized attention weights equals 1. The attention weights are multiplied by the causal strength of the corresponding edges and normalized to obtain the weighted attention coefficients. The causal strength comes from the causal effect coefficients estimated by the two-stage least squares method in step one. The multiplication operation integrates the data-driven attention weights with the knowledge-driven causal strength. After multiplication, all neighbors of the current node are normalized again to ensure that the sum of the weighted attention coefficients equals 1.

[0080] Next, the transformed features of neighboring nodes are weighted and summed based on the weighted attention coefficients. This weighted summation process iterates through all neighboring nodes of the current node, multiplying the transformed feature vector of each neighboring node by its corresponding weighted attention coefficient. All products are then summed dimensionally to obtain the aggregated feature vector. The current node's features are updated after processing with an activation function. The chosen activation function is the ELU function, which remains unchanged for positive numbers and applies a natural exponentiation minus 1 to negative numbers. This non-linear transformation enhances the model's expressive power, and the activation function output serves as the feature representation of the current node in the next layer. Multiple graph convolutional layers are stacked to repeat the neighbor information aggregation and node feature update process. Three graph convolutional layers are used: the first transforms the initial node features and aggregates neighbor information; the second transforms the output features of the first layer again and aggregates information from more distant neighbors; and the third performs the final transformation on the output features of the second layer. Each graph convolutional layer repeats the entire process of linear transformation, attention calculation, weighted summation, and activation function. The final node representation is obtained as the node feature vector output by the third layer graph convolution. The final node representation integrates the node's own features and the topological structure information of its multi-hop neighbors.

[0081] Then, the propagation intensity of the source node is initialized to 1, while that of the other nodes is initialized to 0. The source node refers to the device node in the causal relationship graph that has been identified as being in an abnormal state. The propagation intensity vector is a one-dimensional vector with the same length as the number of nodes. The source node is assigned a value of 1 to indicate that the fault signal started propagating from that node, while the other nodes are assigned a value of 0 to indicate that they have not yet been affected by the fault. A normalized adjacency matrix is ​​used to iteratively propagate the propagation intensity vector. The adjacency matrix is ​​a matrix representation of the causal relationship graph. A matrix element of 1 indicates that there is a directed edge between the corresponding nodes, and a matrix element of 0 indicates that there is no directed edge. The normalized adjacency matrix is ​​symmetrically normalized by the degree matrix. The degree matrix is ​​a diagonal matrix, and the diagonal elements are the degree of the corresponding node, i.e., the number of connecting edges. The iterative propagation process multiplies the normalized adjacency matrix with the propagation intensity vector to obtain the propagation result. The multiplication operation spreads the propagation intensity value of each node to its neighboring nodes. The matrix-vector multiplication is repeated to achieve multi-hop propagation. To prevent over-propagation, a restart probability is introduced, set to 0.15. The iterative propagation formula multiplies the propagation result by 1, subtracts the restart probability, and adds the initial propagation intensity vector multiplied by the restart probability. This mechanism ensures that a certain proportion of the fault signal always returns to the source node during propagation, avoiding excessive accumulation of propagation intensity at distant nodes. After a preset number of iterations (20), a converged propagation intensity is obtained. After 20 iterations, the change in the propagation intensity vector decreases to a negligible level. Each element in the converged propagation intensity vector represents the intensity of the corresponding node's impact from the source node's fault.

[0082] Finally, the node representation, propagation strength, and prediction time window are concatenated and used to calculate the failure risk of each node through a fully connected layer. The node representation comes from the output feature vector of the aforementioned graph convolutional network, the propagation strength comes from the convergence result of the aforementioned iterative propagation algorithm, and the prediction time window is a scalar value representing how long in the future the failure risk is predicted. The concatenation process combines the node representation vector, the corresponding propagation strength scalar, and the prediction time window scalar in sequence into a new vector. The length of the new vector is equal to the node representation dimension plus 2. The fully connected layer contains a weight matrix and a bias vector. The concatenated vector is multiplied by the weight matrix and then the bias vector is added. The result is then mapped to the interval 0 to 1 using the sigmoid activation function. The mapped value is used as the failure risk probability value of that node. The sigmoid function is calculated by dividing 1 by 1 and adding the negative power of the natural exponent. Extract the node sequence with fault risk exceeding a set probability threshold as the fault propagation path. The probability threshold is set to 0.3. Traverse the fault risk probability value of all nodes and sort the nodes with probability values ​​greater than 0.3 from high to low fault risk. The sorted node sequence forms the fault propagation path according to the directed edge connection relationship in the causal relationship graph. The starting point of the path is the source node, i.e. the initial abnormal device, the intermediate nodes of the path are the devices that may be affected, and the ending point of the path is the end device of the fault propagation chain.

[0083] In one specific embodiment, step S5 includes:

[0084] The time derivative of the health feature vector is calculated from the historical time series data. The shape and scale parameters of the Arrhenius-Weibull function are estimated by the maximum likelihood estimation method to obtain the equipment degradation trend function.

[0085] The current health value is calculated based on the health status. The current health value is then substituted into the proportional risk function to exponentially correct the baseline failure rate, thereby obtaining the corrected equipment failure rate.

[0086] The reliability function is calculated by time-domain integration of the corrected equipment failure rate, and the median and upper confidence bound are extracted by Monte Carlo simulation to obtain the remaining service life of the equipment.

[0087] The remaining service life of the equipment is correlated with the fault propagation path, and risk levels are divided and labeled according to the remaining service life value to generate a graded early warning report.

[0088] Specifically, the degradation rate information is obtained by calculating the time derivative of the historical time series data of the health feature vector. The historical time series data of the health feature vector is arranged in chronological order to form a vector sequence of multiple time points. The method for calculating the time derivative is to perform a difference operation on the health feature vectors of two adjacent time points. The difference operation calculates the feature value of the next time point minus the feature value of the previous time point in each dimension. The difference result is divided by the time interval between the two time points to obtain the time derivative vector. The time derivative vector reflects the rate of change of the health feature over time, i.e., the degradation rate. The shape parameter and scale parameter of the Arrhenius-Weibull function are estimated by the maximum likelihood estimation method. The Arrhenius-Weibull function is a physical model describing the change of equipment failure rate over time. The shape parameter determines the failure mode type: a shape parameter less than 1 indicates early failure, a shape parameter equal to 1 indicates random failure, and a shape parameter greater than 1 indicates wear failure. The scale parameter represents the characteristic lifetime, i.e., the time when the failure rate reaches a specific level. The maximum likelihood estimation method constructs a likelihood function and solves for the parameter combination that maximizes the likelihood function. The likelihood function is constructed based on the equipment failure time and censored data from historical operating data. The failure time refers to the actual point in time when the equipment fails, and the censored data refers to the point in time when the observation ends before the equipment fails. Taking the logarithm of the likelihood function converts it to a log-likelihood function to simplify the calculation. The partial derivatives of the shape and scale parameters are calculated separately for the log-likelihood function and set to zero. The estimated values ​​of the shape and scale parameters are obtained by solving the nonlinear equations using Newton's iteration method. The estimated values ​​are then substituted into the Arrhenius-Weibull function to obtain the equipment degradation trend function.

[0089] Next, the current health score is calculated based on the health status, which is derived from the Mahalanobis distance-based determination in step three. The current health score is calculated using the ratio of the Mahalanobis distance to the critical failure threshold. The calculation method is 1 minus the ratio of the squared Mahalanobis distance to the critical failure threshold. When the Mahalanobis distance is close to zero, a health score close to 1 indicates that the equipment is fully healthy; when the Mahalanobis distance is close to the critical failure threshold, a health score close to 0 indicates that the equipment is nearing the edge of failure. The current health score is then substituted into the proportional hazards function to exponentially correct the baseline failure rate. The proportional hazards function describes the impact mechanism of health status on the failure rate. The baseline failure rate is the failure rate without considering health status correction, determined by the value of the aforementioned Arrhenius-Weibull function at the current moment. The index correction process involves multiplying the baseline failure rate by a natural exponential function. The exponential term of this function is composed of the health impact coefficient and 1 minus the current health value. The health impact coefficient reflects the accelerating effect of deteriorating health on failure risk. This coefficient is estimated from historical failure cases using a Cox regression model. The Cox regression model analyzes the relationship between changes in health and the time of failure, and the health impact coefficient is calculated using partial likelihood estimation. The corrected equipment failure rate is the actual failure rate considering the current health status. The corrected failure rate increases with decreasing health, reflecting the accelerating effect of equipment degradation on failure risk.

[0090] Then, the reliability function is calculated by time-domain integration of the corrected equipment failure rate. The reliability function describes the probability that the equipment will still function normally after a given time. The calculation method involves definite integration of the corrected equipment failure rate from the current time to a future time, taking the negative result, and then performing natural exponentiation to obtain the reliability function value. The expected value of the remaining useful life is calculated by time-domain integration of the reliability function from the current time to infinity; this integral value is the expected value of the remaining useful life. The median and upper confidence bound are extracted using Monte Carlo simulation to account for the uncertainty in parameter estimation. The Monte Carlo simulation method performs random sampling within the confidence intervals of the shape and scale parameters. The confidence intervals are calculated based on the asymptotic normality property of maximum likelihood estimation and the Fisher information matrix. Each sampling yields a set of parameter values, which are then input into the reliability function and lifetime calculation formula. This process is repeated 10,000 times to obtain 10,000 predicted remaining lifetime values. These 10,000 predicted values ​​are then sorted from smallest to largest. The 5,000th value after sorting is the median, representing the point estimate of the remaining lifetime. The 9,000th value after sorting is the 90% confidence upper bound, representing the uncertainty boundary of the remaining lifetime. The remaining useful life of the equipment is thus obtained by combining the point estimate and the confidence upper bound.

[0091] Finally, the remaining service life of the equipment is correlated with the fault propagation path. The fault propagation path comes from the node sequence predicted by the graph neural network in step four. The correlation process traverses each equipment node in the fault propagation path and queries the predicted remaining service life value of the corresponding equipment. Risk levels are categorized according to the remaining service life value: red (severe risk), orange (moderately abnormal), and yellow (slightly deteriorated). The classification rule is as follows: equipment with a remaining service life point estimate of less than 30 days is marked as red (severe risk); equipment with a remaining service life point estimate between 30 and 60 days is marked as orange (moderately abnormal); and equipment with a remaining service life point estimate greater than 60 days but a health level below 0.8 is marked as yellow (slightly deteriorated). The hierarchical labeling process labels each equipment node in the fault propagation path with its risk level color and remaining service life value. The labeling content also includes the current health level, the mean of the time derivative of the degradation trend, the predicted remaining service life point estimate, and the confidence upper bound. The generated tiered early warning report includes three levels: device-level early warning information, path-level propagation information, and system-level comprehensive risk indicators. The device-level early warning information lists all devices with a health level below 0.8 and their detailed parameters. The path-level propagation information draws a topology map of the fault propagation path and marks the trigger probability and expected propagation delay of each node. The system-level comprehensive risk indicators calculate the overall risk score of the entire site through weighted summation.

[0092] Figure 2 This is a schematic diagram illustrating the device degradation trend analysis based on the historical evolution trajectory of health feature vectors in an embodiment of this application. Figure 2 This paper illustrates the evolution of health status values ​​for three typical secondary devices over a 30-day operating cycle. Device A exhibits a normal degradation trend, with its health status steadily decreasing from an initial 1.0 to approximately 0.85, consistent with the wear failure mode described by the Arrhenius-Weibull function. Device B shows accelerated degradation after day 15, with its health status plummeting below 0.4, entering the red critical vulnerability zone. This inflection point can be accurately captured by calculating the time derivative of the health feature vector. Device C maintains stable operation, with its health status fluctuating within the healthy range above 0.9. This figure verifies that the health feature vector constructed by fusing waveform characteristics, frequency domain response characteristics, and device performance monitoring data can effectively quantify the degree of device degradation, providing a reliable basis for health status assessment and remaining life prediction based on Mahalanobis distance, and realizing a shift from passive post-event analysis to proactive predictive maintenance.

[0093] The above describes the method for identifying potential hazards in secondary equipment based on waveform recordings and monitoring data in the embodiments of this application. The following describes the system for identifying potential hazards in secondary equipment based on waveform recordings and monitoring data in the embodiments of this application. Please refer to [link / reference]. Figure 3 One embodiment of the secondary equipment hidden danger detection system based on waveform recording files and monitoring data in this application includes:

[0094] The parsing module is used to parse the waveform file to obtain the sampling sequence, combine the equipment topology relationship and protection logic relationship, identify the causal relationship between the equipment by calculating the difference between the conditional probability and the intervention probability, quantify the causal strength, and construct a secondary equipment causal association map;

[0095] The extraction module is used to perform full-cycle Fourier transform on the sampled sequence to calculate the amplitude and phase of each harmonic, and to perform wavelet decomposition to extract the energy distribution and statistics of each scale layer, thereby obtaining multi-scale waveform feature vectors.

[0096] The recording module is used to inject small signal test excitation into the secondary equipment and record the response output. It estimates the transfer function of the input-output relationship through the system identification method, extracts the frequency domain response characteristics and fuses them with the multi-scale waveform feature vector and equipment performance monitoring data to form a health feature vector, and calculates the distance of the health feature vector from the health benchmark to determine the health status of the equipment.

[0097] The input module is used to input the causal association graph of the secondary equipment and the health status into the graph neural network, calculate the weighted attention coefficient between nodes by combining the causal strength, update the node representation by aggregating neighbor node information, iteratively calculate the fault risk of each node based on the node representation and propagation strength, and identify the fault propagation path.

[0098] The prediction module is used to establish a device degradation model based on the historical evolution trajectory of the health feature vector, predict the remaining service life of the device according to the corrected device failure rate of the health status, and generate a graded early warning report in combination with the fault propagation path.

[0099] above Figure 3 The secondary equipment hidden danger discovery system based on waveform recording files and monitoring data in this embodiment of the invention will be described in detail from the perspective of modular functional entities. The secondary equipment hidden danger discovery device based on waveform recording files and monitoring data in this embodiment of the invention will be described in detail from the perspective of hardware processing.

[0100] Reference Figure 4 This invention also provides a secondary equipment hazard discovery device based on waveform recording files and monitoring data. This device can be a server, and its internal structure can be as follows: Figure 4As shown, the secondary equipment hazard detection device based on waveform recording files and monitoring data includes a processor, memory, display screen, input device, network interface, and database connected via a system bus. The processor, designed as a computer, provides computing and control capabilities. The memory of the device includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system, computer programs, and database. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The database stores the data corresponding to this embodiment. The network interface of the device is used for communication with external terminals via a network connection. When the computer program is executed by the processor, it implements the above-described method.

[0101] Those skilled in the art will understand that Figure 4 The structure shown is merely a block diagram of a portion of the structure related to the present invention and does not constitute a limitation on the secondary equipment hidden danger discovery equipment based on waveform files and monitoring data applied to the present invention.

[0102] The present invention also provides a computer-readable storage medium, which can be a non-volatile computer-readable storage medium or a volatile computer-readable storage medium, wherein the computer-readable storage medium stores instructions that, when the instructions are executed on a computer, cause the computer to perform the steps of the secondary equipment hidden danger discovery method based on waveform files and monitoring data.

[0103] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0104] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a secondary equipment hazard detection device (which can be a personal computer, server, or network device, etc.) based on waveform files and monitoring data to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0105] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. 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 of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A secondary equipment hidden danger mining method based on recording file and monitoring data, characterized in that, The method comprises: S1 step: parse the recording file to obtain a sampling sequence, combine the device topology relationship and the protection logic relationship, identify the causal relationship between devices by calculating the difference between the conditional probability and the intervention probability, and quantify the causal strength, construct a secondary device causal correlation graph, including: according to the configuration file of the standard recording file, the number of analog channels, the number of switch channels and the sampling frequency are obtained, according to the file format type, the time sequence sampling point data of each channel is read from the data file to form the sampling sequence; for any two device nodes, the conditional probability is calculated by statistically analyzing the abnormal co-occurrence frequency of the device from the historical operation data, the intervention probability is calculated by observing the change of the device state under the condition of controlling other variables, when the difference between the conditional probability and the intervention probability exceeds the set threshold, the instrumental variable is introduced, and the causal effect coefficient is estimated as the causal strength by using the two-stage least squares method; a directed graph is constructed as the secondary device causal correlation graph, wherein the node represents the secondary device, the directed edge represents the causal relationship, and the edge weight is calculated by weighted summation of the statistical co-occurrence probability, the causal strength and the satisfaction degree of the protection logic relationship; S2 step: full cycle Fourier transform is performed on the sampling sequence to calculate the amplitude and phase of each harmonic, and wavelet decomposition is performed to extract the energy distribution and statistical quantity of each scale layer to obtain a multi-scale waveform feature vector; S3 step: small signal test excitation is introduced into the secondary device, and the response output is recorded, the transfer function of the input-output relationship is estimated by system identification method, the frequency domain response characteristics are extracted, and the multi-scale waveform feature vector and the device performance monitoring data are fused to form a health feature vector, and the distance of the health feature vector deviating from the health benchmark is calculated to determine the health state of the device; S4 step: input the secondary device causal correlation graph and the health state into the graph neural network, calculate the weighted attention coefficient between nodes combined with the causal strength, update the node representation by aggregating the neighbor node information, and iteratively calculate the fault risk of each node based on the node representation and the propagation strength to identify the fault propagation path; S5 step: establish a device degradation model based on the historical evolution trajectory of the health feature vector, calculate the modified device failure rate according to the health state, predict the remaining useful life of the device, and generate a graded warning report combined with the fault propagation path, including: calculate the time derivative of the historical time series data of the health feature vector, estimate the shape parameter and scale parameter of the Arrhenius-Weibull function by maximum likelihood estimation method, and obtain the device degradation trend function; calculate the current health degree value according to the health state, substitute the current health degree value into the proportional risk function to exponentially modify the benchmark failure rate, and obtain the modified device failure rate; calculate the reliability function by time domain integration of the modified device failure rate, extract the median and confidence upper bound by Monte Carlo simulation, and obtain the remaining useful life of the device; correlate the remaining useful life of the device with the fault propagation path, divide the risk level according to the remaining life value and label it in layers, and generate the graded warning report.

2. The method according to claim 1, wherein, The S2 step comprises: applying a full Fourier algorithm to the sampling sequence to calculate the sine term coefficient and the cosine term coefficient of the sampling points in a power frequency cycle, and calculating the amplitude and phase of each harmonic according to the sine term coefficient and the cosine term coefficient; selecting a wavelet mother function to perform multi-layer wavelet decomposition on the sampling sequence, discretizing the scale parameter into different decomposition layers, each decomposition layer corresponding to a specific frequency range, and calculating the wavelet coefficient of each decomposition layer; calculating the energy value as the energy distribution according to the wavelet coefficient of each decomposition layer, and calculating the mean, standard deviation, skewness and kurtosis as the statistics; splicing the amplitude and phase of each harmonic, the energy distribution and the statistics to form the multi-scale waveform feature vector.

3. The method according to claim 2, wherein, The selection of the wavelet mother function to perform multi-layer wavelet decomposition on the sampling sequence, discretize the scale parameter into different decomposition layers, each decomposition layer corresponding to a specific frequency range, and calculate the wavelet coefficient of each decomposition layer comprises: selecting a wavelet mother function to perform continuous wavelet transform on the time series signal in the sampling sequence, wherein the scale parameter controls the frequency resolution and the translation parameter controls the time positioning; discretize the scale parameter according to the power of 2, set the maximum number of decomposition layers to the binary logarithm value of the total number of sampling points, and calculate the frequency range corresponding to each layer by dividing the sampling frequency by the power of the scale parameter; convolve the wavelet mother function on each decomposition layer to calculate the wavelet transform coefficient of the sampling sequence at each time position and scale layer as the wavelet coefficient; extract the amplitude sequence of the wavelet coefficient of each decomposition layer for subsequent energy distribution and statistics calculation.

4. The method according to claim 1, wherein, The S3 step comprises: generating a sine waveform with an amplitude of 0.1%-1% of the rated working level as the small signal test excitation, inputting the small signal test excitation into the secondary setting, and recording the device output signal as the response output through the data acquisition device; establishing an autoregressive moving average model to describe the input-output relationship, identifying the parameters of the small signal test excitation and the response output using the recursive least squares algorithm, estimating the polynomial coefficients of the model, and calculating the transfer function according to the polynomial coefficients; performing frequency domain transformation on the transfer function, calculating the amplitude-frequency response and phase-frequency response at the preset characteristic frequency points, and combining the amplitude-frequency response values and phase-frequency response values at each characteristic frequency point to form a frequency domain feature vector as the frequency domain response characteristic; splicing the frequency domain response characteristic, the multi-scale waveform feature vector and the device performance monitoring data to form the health feature vector, collecting multiple sets of health feature vectors of healthy devices to calculate the mean vector and covariance matrix as a health baseline model, calculating the Mahalanobis distance between the health feature vector of the device to be evaluated and the mean vector, and determining the abnormal state when the Mahalanobis distance exceeds the detection threshold set based on the chi-square distribution as the device health state.

5. The method according to claim 1, wherein, The S4 step comprises: The node initial feature of the secondary equipment causal correlation graph is set as the health feature vector of the corresponding equipment, linear transformation is performed on each node initial feature through a shared weight matrix, an attention function is applied to calculate the original attention score of a node to its neighbor nodes, the original attention score is normalized to obtain an attention weight, the attention weight is multiplied by the causal strength of the corresponding edge and normalized to obtain the weighted attention coefficient; The transformed features of the neighbor nodes are weighted and summed according to the weighted attention coefficient, and the current node feature is updated after being processed by an activation function, the neighbor information aggregation and node feature updating process are repeatedly performed by stacking multiple graph convolution layers to obtain the final node representation; The propagation strength of the source node is initialized as 1 and the propagation strength of the rest nodes is 0, the propagation strength vector is iteratively propagated by using a normalized adjacency matrix, a restart probability is introduced to prevent overdiffusion, and the converged propagation strength is obtained after a preset number of iterations; The node representation, the propagation strength and a prediction time window are spliced to calculate the fault risk of each node through a fully connected layer, and a node sequence with a fault risk exceeding a set probability threshold is extracted as the fault propagation path.

6. A secondary equipment hidden danger mining system based on recording file and monitoring data, characterized in that, The secondary equipment hidden danger mining system based on the recording file and monitoring data comprises: An analysis module is configured to analyze the recording file to obtain a sampling sequence, identify the causal relationship between devices and quantify the causal strength by calculating the difference between the conditional probability and the intervention probability in combination with the device topology relationship and the protection logic relationship, and construct a secondary equipment causal correlation graph; An extraction module is configured to perform full-cycle Fourier transform on the sampling sequence to calculate the amplitude and phase of each harmonic, and perform wavelet decomposition to extract the energy distribution and statistical quantity of each scale layer to obtain a multi-scale waveform feature vector; A recording module is configured to input a small signal test excitation into a secondary device and record the response output, estimate the transfer function of the input-output relationship by a system identification method, extract the frequency domain response characteristic, and fuse the multi-scale waveform feature vector and the device performance monitoring data to form a health feature vector, and calculate the distance of the health feature vector from a health benchmark to determine the health state of the device; An input module is configured to input the secondary equipment causal correlation graph and the health state into a graph neural network, calculate the weighted attention coefficient between nodes in combination with the causal strength, update the node representation by aggregating the neighbor node information, and iteratively calculate the fault risk of each node based on the node representation and the propagation strength to identify the fault propagation path; A prediction module is configured to establish a device degradation model based on the historical evolution trajectory of the health feature vector, predict the remaining useful life of the device according to the modified device failure rate of the health state, and generate a graded early warning report in combination with the fault propagation path.

7. A device for mining secondary equipment hidden dangers based on recording file and monitoring data, characterized in that, The computer program is run on the processor, and the processor executes the computer program to realize the secondary equipment hidden danger mining method based on the recording file and monitoring data according to any one of claims 1 to 5.

8. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is run on the processor, and the processor executes the computer program to realize the secondary equipment hidden danger mining method based on the recording file and monitoring data according to any one of claims 1 to 5.

Images