Method for detecting winding faults of power transformer based on vibration spectrum analysis

Abstract

The application discloses a power transformer winding fault detection method based on vibration spectrum analysis and relates to the technical field of power equipment detection, which comprises the following steps: collecting original vibration acceleration signals of a transformer and preprocessing the original vibration acceleration signals to form a signal frame sequence; applying an improved Hilbert-Huang transform to each signal frame to obtain a set of intrinsic mode functions, calculating a component instantaneous frequency sequence and extracting a frequency band with an energy proportion exceeding a preset threshold as a candidate fault characteristic frequency band; extracting a subband signal of each frequency band and calculating a high-order statistical quantity feature matrix; loading a health state reference feature matrix matched with a transformer model, capacity and voltage grade; comparing the two in an element-by-element manner to identify an abnormal feature submatrix; spatially clustering the abnormal feature submatrix to form a suspected fault mode cluster; and outputting a specific winding fault mode diagnosis result. The method can accurately extract fault characteristics and improve the accuracy and reliability of fault detection.

Description

Technical Field

[0001] This invention belongs to the field of power equipment testing technology, specifically a method for detecting faults in power transformer windings based on vibration spectrum analysis. Background Technology

[0002] Power transformers are core energy conversion devices in power systems. Their windings, as critical components, are susceptible to loosening, deformation, and displacement due to factors such as electromagnetic forces, temperature changes, and mechanical vibrations. Failure to detect and diagnose these faults in a timely manner can lead to equipment outages or even power grid accidents. Currently, transformer winding fault detection based on vibration spectrum analysis is one of the mainstream technologies. This type of technology typically collects transformer vibration signals, processes the signals using traditional Hilbert-Huang transforms or Fourier transforms to extract frequency band features, and then combines this with preset health reference data for fault identification.

[0003] In traditional techniques, the Hilbert-Huang transform is prone to mode aliasing when processing complex transformer vibration signals, resulting in weakly targeted frequency band features. This makes it difficult to effectively filter out frequency bands directly related to faults, and a large number of irrelevant frequency bands can interfere with the accuracy of fault identification. Furthermore, existing technologies often use a generic health status reference feature matrix without matching it to the specific transformer model, capacity, and voltage level. This leads to significant errors in feature comparison, making it difficult to accurately identify subtle abnormal features. Moreover, the identified abnormal features lack effective classification processing, failing to clearly distinguish different types of winding fault modes, and easily resulting in misdiagnosis or missed fault detection.

[0004] Transformer vibration signals contain a large amount of equipment operating noise and irrelevant interference signals. How to accurately filter out the fault-related frequency bands, and at the same time achieve feature comparison with the equipment's own parameters and effective classification of abnormal features, has become a key challenge to improve the accuracy of winding fault detection. Summary of the Invention

[0005] This invention aims to solve at least one of the technical problems existing in the prior art;

[0006] Therefore, this invention proposes a power transformer winding fault detection method based on vibration spectrum analysis, including:

[0007] The original vibration acceleration signal of the transformer is acquired and preprocessed to obtain the preprocessed vibration time-domain signal, forming multiple signal frame sequences;

[0008] For each of the signal frames, the improved Hilbert-Huang transform is applied to calculate the set of eigenmode functions of the signal frame;

[0009] For each eigenmode function component in the set of eigenmode functions, calculate its instantaneous frequency sequence;

[0010] From the instantaneous frequency sequence of all intrinsic mode function components, frequency bands with energy proportions exceeding a preset threshold are extracted as a set of candidate fault feature frequency bands;

[0011] For each frequency band in the candidate fault feature frequency band set, the corresponding sub-band signal is extracted from the preprocessed vibration time domain signal, and its higher-order statistical feature matrix is ​​calculated.

[0012] Load a health status reference feature matrix that matches the current transformer model, capacity, and voltage level;

[0013] The high-order statistical feature matrix of all candidate fault feature bands is calculated and compared element by element with the health status reference feature matrix to identify abnormal feature sub-matrices with differences exceeding a preset threshold.

[0014] Spatial clustering is performed on the abnormal feature sub-matrices to group the strongly correlated abnormal feature sub-matrices into one class, forming several independent suspected fault mode clusters. Based on the central features and frequency band distribution of each cluster, specific winding fault mode diagnosis results are output.

[0015] Furthermore, the original vibration acceleration signal of the acquired transformer is preprocessed to obtain a preprocessed vibration time-domain signal, forming multiple signal frame sequences, including:

[0016] Vibration acceleration sensors are installed at designated locations on the outer wall of the oil tank of the power transformer to collect the original vibration acceleration signals of the transformer under steady-state operating conditions.

[0017] The original vibration acceleration signal is preprocessed, including removing the DC component, bandpass filtering, and trend term elimination, to obtain the preprocessed vibration time-domain signal.

[0018] The preprocessed vibration time-domain signal is divided into equal-length frames to form a continuous sequence of multiple signal frames, specifically including:

[0019] Set a fixed time window length as the frame length;

[0020] Set a window movement step size smaller than the frame length as the frame shift;

[0021] Starting from the beginning of the preprocessed vibration time-domain signal, the first segment of the signal is extracted with a frame length as its width, and this segment is taken as the first signal frame.

[0022] Move the capture window backward by one frame shift length, capture the second segment of the signal, and use it as the second signal frame;

[0023] Repeat the window movement and capture operation until the end of the signal, thus forming a sequence of multiple signal frames that are continuous in time and partially overlap.

[0024] For the last segment of a signal frame whose length is less than the frame length, zero values ​​are used to pad the end so that its length reaches the set frame length.

[0025] Furthermore, the improved Hilbert-Huang transform includes:

[0026] Local extremum detection is performed on the input signal frame to identify all maximum and minimum points;

[0027] Construct the upper envelope using all maxima and the lower envelope using all minima.

[0028] Calculate the mean of the upper and lower envelopes to obtain the original mean envelope;

[0029] Subtracting the original mean envelope from the input signal frame yields a preliminary decomposed component;

[0030] For the initially decomposed components, calculate their instantaneous amplitude at each sampling point, and perform extreme point detection and envelope construction on the instantaneous amplitude sequence again to obtain the amplitude modulation envelope of the components.

[0031] The amplitude modulation envelope is separated from the initially decomposed components, thereby decomposing the original signal into the product of a pure frequency modulation component and an amplitude modulation component.

[0032] The pure frequency-modulated component is output as an eigenmode function component that satisfies the condition.

[0033] For the separated amplitude modulation envelope, the steps of detecting from local extrema and component decomposition are repeated until it satisfies monotonicity or a preset stopping criterion. The resulting series of pure frequency modulation components constitute the intrinsic mode function set.

[0034] Further, for each eigenmode function component in the set of eigenmode functions, its instantaneous frequency sequence is calculated, including:

[0035] For any of the pure frequency modulation intrinsic mode function components, perform a Hilbert transform to obtain its orthogonal analytic signal components;

[0036] A complex analytic signal is constructed from the eigenmode function components of the pure frequency modulation and its analytic signal components;

[0037] Calculate the phase angle of the complex analytic signal at each sampling point;

[0038] The phase angle between consecutive sampling points is differentially calculated to obtain the change in phase angle;

[0039] Divide the change in phase angle by the product of twice pi and the sampling time interval to obtain the instantaneous frequency value of the intrinsic mode function component at each sampling point.

[0040] Following the steps of Hilbert transform, constructing complex analytic signals, calculating phase angles, calculating phase angle changes using differential methods, and calculating instantaneous frequency values, each component in the intrinsic mode function set is processed sequentially to obtain the instantaneous frequency sequence of all components.

[0041] Furthermore, from the instantaneous frequency sequences of all intrinsic mode function components, frequency bands with an energy proportion exceeding a preset threshold are extracted as a set of candidate fault feature frequency bands, including:

[0042] For each component in the set of intrinsic mode functions, calculate its energy over the entire signal frame time;

[0043] The energy of all intrinsic mode function components is normalized to obtain the energy percentage of each component.

[0044] A threshold for energy percentage is set, and intrinsic mode function components whose energy percentage exceeds the threshold are selected as dominant energy components.

[0045] For each dominant energy component, the histogram distribution of its instantaneous frequency sequence across the entire frequency range is statistically analyzed.

[0046] In the histogram distribution, identify the frequency range where the number of frequency points reaches a peak;

[0047] The frequency ranges corresponding to each dominant energy component are merged, and the overlapping ranges are combined to form the candidate fault feature frequency band set.

[0048] Furthermore, for each frequency band in the candidate fault feature frequency band set, the corresponding sub-band signal is extracted from the preprocessed vibration time-domain signal, and its higher-order statistical feature matrix is ​​calculated, including:

[0049] The preprocessed vibration time-domain signal is filtered using a bandpass filter with a center frequency equal to the center frequency of the frequency band and a bandwidth equal to the width of the frequency band to obtain the sub-band signal of the frequency band.

[0050] For the sub-band signal, calculate its third moment, i.e., skewness;

[0051] For the sub-band signal, calculate its fourth moment, i.e., kurtosis;

[0052] For the sub-band signal, calculate its fifth moment, i.e., the hyperkurtosis;

[0053] For the sub-band signal, calculate its sixth moment;

[0054] From the sub-band signal, extract multiple segments of equal length, calculate the amplitude square of each segment, and then calculate the kurtosis of the amplitude square value to obtain the amplitude square kurtosis.

[0055] From the subband signal, extract multiple segments of equal length, calculate the kurtosis of the envelope spectrum of each segment, and obtain the envelope kurtosis;

[0056] The calculated skewness, kurtosis, hyperkurtosis, sixth moment, amplitude squared kurtosis, and envelope kurtosis are arranged in rows to form a six-row, one-column feature matrix, which serves as the higher-order statistical feature matrix of the frequency band.

[0057] Furthermore, a health status reference feature matrix matching the current transformer model, capacity, and voltage level is loaded, including:

[0058] Load a priori knowledge base containing vibration characteristic data collected and calculated from transformers of various models, capacities, and voltage levels under confirmed healthy conditions;

[0059] Based on the nameplate information of the transformer being tested, including model, capacity, and voltage level, a matching query is performed in the prior knowledge base;

[0060] Extract the set of high-order statistical feature matrices of all healthy samples that are completely consistent with the current transformer parameters from the knowledge base;

[0061] Calculate the average value of all feature matrices in the set of higher-order statistics feature matrices for each feature dimension;

[0062] The average values ​​of each feature dimension are combined into a six-row, one-column average feature matrix, which serves as the health status reference feature matrix for the current detection task.

[0063] Furthermore, the high-order statistical feature matrix of all candidate fault feature bands calculated is compared element-by-element with the health status reference feature matrix, including:

[0064] Let the high-order statistical feature matrix of a candidate frequency band to be detected be denoted as matrix D, and let the health status reference feature matrix be denoted as matrix H. Matrix D and matrix H have the same dimension.

[0065] Calculate the absolute difference between the elements at corresponding positions of matrix D and matrix H;

[0066] Divide the absolute difference by the absolute value of the element at the corresponding position in matrix H to obtain the relative difference of each feature dimension;

[0067] Compare the calculated relative difference with the preset difference threshold;

[0068] If the relative difference of a certain feature dimension exceeds the difference threshold, then the feature dimension is marked as an abnormal feature dimension.

[0069] Extract the row vectors containing all the feature dimensions marked as abnormal from matrix D, and combine the row vectors into a new submatrix, which is the abnormal feature submatrix.

[0070] Further, spatial clustering is performed on the abnormal feature submatrix, including:

[0071] Each abnormal feature submatrix corresponds to one or more abnormal feature dimensions in a candidate fault feature frequency band;

[0072] Each anomalous feature submatrix is ​​considered as a feature point in a multidimensional space, and its feature dimension coordinates are composed of the row vector values ​​of the anomalous feature submatrix.

[0073] Calculate the cosine similarity between any two anomalous feature submatrices as a measure of the distance between them;

[0074] Set a similarity threshold as a condition for merging clusters;

[0075] Initialize all anomalous feature submatrices as their own independent clusters;

[0076] Calculate the average cosine similarity between all pairs of clusters, and find the pair of clusters with the highest average similarity;

[0077] If the average similarity of the pair of clusters with the highest average similarity is greater than the similarity threshold, then the two clusters are merged into a new cluster, and the center feature point of the new cluster is updated to the weighted average of the center feature points of the original two clusters.

[0078] Repeat the steps of calculating the average cosine similarity between all pairs of clusters, finding the pair of clusters with the highest average similarity, and merging clusters when the conditions are met, until the average similarity between any two clusters is higher than the aforementioned similarity threshold.

[0079] Each remaining cluster is then considered an independent cluster representing a suspected fault mode.

[0080] Furthermore, based on the central characteristics and frequency band distribution of each cluster, specific winding fault mode diagnostic results are output, including:

[0081] For each cluster of suspected fault modes that is ultimately formed, calculate the central eigenvector of all its member anomalous feature submatrices;

[0082] Extract the center frequency of the candidate fault feature band corresponding to all members in the suspected fault mode cluster;

[0083] The distribution range of the center frequency is statistically analyzed, and its minimum and maximum center frequencies are identified to define the fault characteristic frequency coverage range of the suspected fault mode cluster.

[0084] The central feature vector and the fault feature frequency coverage interval are matched with a preset fault feature rule base.

[0085] The fault feature rule base stores different winding fault modes, such as inter-turn short circuit, winding deformation, insulation deterioration, and loosening of clamping force, as well as their corresponding typical feature vector patterns and feature frequency ranges.

[0086] Search the rule base for the fault mode entry that is closest to the central feature vector and whose fault feature frequency coverage range overlaps with the fault mode entry.

[0087] The searched fault mode entries are output as the diagnostic results of the suspected fault mode cluster.

[0088] Compared with the prior art, the beneficial effects of the present invention are:

[0089] An improved Hilbert-Huang transform is applied to the acquired and preprocessed vibration signal frame sequence to calculate the set of intrinsic mode functions (EMFs) of the signal frame. Then, by analyzing the instantaneous frequency sequences of each EMF component, frequency bands with energy proportions exceeding a preset threshold are extracted as a set of candidate fault feature frequency bands. This technique effectively improves the mode aliasing problem inherent in the traditional Hilbert-Huang transform, accurately filters out frequency bands related to winding faults, excludes frequency bands corresponding to equipment operating noise and irrelevant interference signals, making fault features more targeted, avoiding interference from irrelevant frequency bands in the fault identification process, and improving the accuracy of fault feature extraction.

[0090] This approach extracts corresponding sub-band signals from candidate fault feature frequency bands and calculates high-order statistical feature matrices. These matrices are then compared element-by-element with a health status reference feature matrix matching the current transformer model, capacity, and voltage level. Abnormal feature sub-matrices with differences exceeding a preset threshold are identified. Spatial clustering of these abnormal feature sub-matrices is then performed, grouping strongly correlated abnormal feature sub-matrices into one category, forming several independent clusters of suspected fault modes, and outputting specific diagnostic results. This technical approach solves the comparison error problem caused by the mismatch between the general health reference matrix and the equipment's own parameters. Element-by-element comparison can capture subtle feature differences, accurately identifying abnormal features inconsistent with the health status. Spatial clustering can distinguish abnormal features corresponding to different types of faults, clearly presenting various suspected fault modes and avoiding misjudgments and omissions caused by confusion of different fault features, thus achieving accurate differentiation and diagnosis of winding fault modes. Attached Figure Description

[0091] Figure 1 This is a flowchart illustrating the steps of the power transformer winding fault detection method based on vibration spectrum analysis described in this invention.

[0092] Figure 2 A flowchart for calculating the instantaneous frequency sequence of intrinsic mode function components;

[0093] Figure 3 Instantaneous frequency distribution curve of the dominant IMF component;

[0094] Figure 4 A bar chart comparing the higher-order statistical characteristics of transformer vibration signals under healthy and fault conditions;

[0095] Figure 5 Heatmap for matching fault clusters with fault modes in the rule base. Detailed Implementation

[0096] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0097] See Figure 1 The original vibration acceleration signal of the transformer is acquired and preprocessed to obtain the preprocessed vibration time-domain signal, forming multiple signal frame sequences. For each signal frame, an improved Hilbert-Huang transform is applied to calculate the intrinsic mode function (EMF) set. For each EMF component in the EMF set, its instantaneous frequency sequence is calculated. From the instantaneous frequency sequences of all EMF components, frequency bands with energy proportions exceeding a preset threshold are extracted as a candidate fault feature frequency band set. For each frequency band in the candidate fault feature frequency band set, the corresponding sub-band signal is extracted from the preprocessed vibration time-domain signal, and its higher-order statistical feature matrix is ​​calculated. A health status reference feature matrix matching the current transformer model, capacity, and voltage level is loaded. The higher-order statistical feature matrices of all candidate fault feature frequency bands are compared element-wise with the health status reference feature matrix to identify abnormal feature sub-matrices with differences exceeding a preset threshold. Spatial clustering is performed on the abnormal feature submatrices to group the strongly correlated abnormal feature submatrices into one class, forming several independent suspected fault mode clusters. Based on the central features and frequency band distribution of each cluster, specific winding fault mode diagnosis results are output.

[0098] In one embodiment of the present invention, a vibration acceleration sensor is installed at a predetermined location on the outer wall of the power transformer tank to collect the raw vibration acceleration signal of the transformer under steady-state operating conditions. The raw vibration acceleration signal is preprocessed, including DC component removal, bandpass filtering, and trend term elimination, to obtain a preprocessed vibration time-domain signal. The preprocessed vibration time-domain signal is then divided into equal-length frames to form a continuous sequence of multiple signal frames. Specifically, a fixed time window length is set as the frame length, and a window movement step size smaller than the frame length is set as the frame shift. Starting from the beginning of the preprocessed vibration time-domain signal, a first segment of the signal is extracted with a width equal to the frame length, serving as the first signal frame. The extraction window is moved backward by one frame shift length to extract a second segment of the signal, serving as the second signal frame. The window movement and extraction operations are repeated until the end of the signal, thus forming a sequence of multiple signal frames that are continuous in time and partially overlap. For the last signal frame segment whose length is less than the frame length, zero values ​​are added to the end to make its length reach the predetermined frame length.

[0099] In specific implementation, the process of acquiring the original vibration acceleration signal of the transformer and preprocessing it to obtain the preprocessed vibration time-domain signal, forming multiple signal frame sequences, is as follows: A high-sensitivity vibration acceleration sensor is installed at a flat and sturdy location on the outer wall of the power transformer tank, with the sensitive axis of the sensor perpendicular to the tank wall surface. At a sampling frequency of 100kHz, the original vibration acceleration signal sequence of the power transformer under rated voltage and rated load steady-state operating conditions is continuously acquired for at least 60 seconds. The original vibration acceleration signal sequence contains complex vibration information resulting from the coupling of electromagnetic force, mechanical vibration, and environmental interference. In some embodiments, the amplitude range of the original vibration acceleration signal is between -5g and +5g, containing significant power frequency and higher harmonic components.

[0100] In practice, the original vibration acceleration signal sequence undergoes preprocessing to remove the DC component. The arithmetic mean of the entire signal sequence is calculated, and this arithmetic mean is subtracted from the value at each sampling point. The signal after removing the DC component is filtered using a Butterworth bandpass filter with a passband frequency of 100Hz to 2000Hz to retain the characteristic frequency bands related to the winding's mechanical vibration while suppressing low-frequency mechanical interference and high-frequency electrical noise. Optionally, the trend term curve of the bandpass-filtered signal is fitted using the least squares method, and this trend term curve is subtracted from the signal, ultimately obtaining a preprocessed vibration time-domain signal with zero mean and no trend term. It can be understood that preprocessing effectively eliminates the interference of sensor baseline drift and low-frequency background vibration on subsequent analysis.

[0101] In some embodiments, the preprocessed vibration time-domain signal is framed, with a fixed time window length (frame length) of 0.1 seconds, corresponding to 10,000 sampling points. A window movement step size smaller than the frame length (frame shift) is set, which can be calculated according to specific needs; for example, the frame shift equals the frame length multiplied by an overlap ratio coefficient of 0.5. Starting from the beginning of the preprocessed vibration time-domain signal, a first segment of the signal is extracted with a width equal to the frame length; this segment is the first signal frame. The extraction window is moved backward by one frame shift length, and a second segment of the signal is extracted; this segment is the second signal frame, which overlaps with the first signal frame by 50%. The window movement and extraction operations are repeated until the end of the signal, thus forming a sequence of multiple signal frames that are temporally continuous and partially overlapping. The generation process of the signal frame sequence can be formally described, and the value of the frame shift is determined by the following formula:

[0102]

[0103] Where: FrameShift represents the number of sampling points for frame shift, FrameLength represents the number of sampling points for frame length, and OverlapRatio represents the preset signal frame overlap ratio. In the last capture window, if the remaining signal length is less than the set frame length, zero values ​​are used to pad the end of the signal to make the last signal frame the full frame length. Optionally, the zero-value padding length is 300 sampling points. The zero-value padding operation ensures that all signal frames have a uniform length, facilitating subsequent uniform transformation and feature calculation. It can be understood that frame segmentation transforms continuous vibration signals into a series of short-time stationary signal segments, which is beneficial for capturing time-varying local spectral features in the signal.

[0104] In one embodiment of the present invention, local extrema detection is performed on the input signal frame to identify all maxima and minima. An upper envelope is constructed using all maxima, and a lower envelope is constructed using all minima. The mean of the upper and lower envelopes is calculated to obtain the original mean envelope. The original mean envelope is subtracted from the input signal frame to obtain a preliminary decomposed component. For the preliminary decomposed component, its instantaneous amplitude at each sampling point is calculated, and the instantaneous amplitude sequence is again subjected to extrema detection and envelope construction to obtain the amplitude modulation envelope of the component. The amplitude modulation envelope is separated from the preliminary decomposed component, thereby decomposing the original signal into the product of a pure frequency-modulated component and an amplitude-modulated component. The pure frequency-modulated component is output as an eigenmode function component that satisfies certain conditions. For the separated amplitude modulation envelope, the steps from local extremum detection to component decomposition are repeated until it satisfies monotonicity or a preset stopping criterion. The resulting series of pure frequency-modulated components constitute the set of eigenmode functions of the signal frame. (See also...) Figure 2For any purely frequency-modulated eigenmode function component, a Hilbert transform is performed to obtain its orthogonal analytic signal component. A complex analytic signal is constructed from this purely frequency-modulated eigenmode function component and its analytic signal component. The phase angle of this complex analytic signal at each sampling point is calculated. The phase angle between consecutive sampling points is differentially calculated to obtain the change in phase angle. The change in phase angle is divided by the product of twice pi and the sampling time interval to obtain the instantaneous frequency value of the eigenmode function component at each sampling point. Following the above steps, each component in the eigenmode function set is processed sequentially to obtain the instantaneous frequency sequence of all components.

[0105] In practical implementation, the process of applying the improved Hilbert-Huang transform to calculate the set of intrinsic mode functions (EMFs) of a signal frame and calculating the instantaneous frequency sequence for each EMF component is as follows: For a single signal frame extracted from a preprocessed vibration time-domain signal sequence, the signal frame is a digital sequence containing N discrete sampling points. Local extremum detection is performed on the input signal frame. Each sampling point in the signal frame is traversed, and points with values ​​greater than the values ​​of the preceding and following sampling points are identified as local maxima, while points with values ​​less than the values ​​of the preceding and following sampling points are identified as local minima. Using all identified local maxima, a smooth upper envelope is constructed using cubic spline interpolation; using all identified local minima, a smooth lower envelope is constructed using cubic spline interpolation. The arithmetic mean of the upper and lower envelopes at each sampling point is calculated to obtain an original mean envelope that runs through the entire signal frame. Subtracting the value of the original mean envelope at the corresponding sampling point from the value of each sampling point in the input signal frame yields a preliminary decomposed component. In some embodiments, the preliminary decomposed component may still contain undesirable amplitude modulation.

[0106] In practical implementation, the instantaneous amplitude of each sample point of the initially decomposed component is calculated. This instantaneous amplitude is calculated by performing a Hilbert transform on the initially decomposed component to obtain its orthogonal component and then calculating its magnitude. The calculated instantaneous amplitude sequence is then subjected to local extremum detection and envelope construction to obtain the amplitude modulation envelope corresponding to the initially decomposed component. The operation of separating the amplitude modulation envelope from the initially decomposed component involves dividing the value of each sample point of the initially decomposed component by the value of the amplitude modulation envelope at the corresponding sample point, thereby decomposing the original signal into the product of a pure frequency-modulated component and an amplitude-modulated component. The pure frequency-modulated component is defined as a component that satisfies the narrowband condition and has a unit amplitude; it is output as an eigenmode function component that satisfies the condition. Optionally, the process of separating the amplitude modulation envelope can be expressed as decomposing the signal into a form where amplitude and frequency are decoupled. The product relationship between the pure frequency-modulated component and the amplitude-modulated component is described by the following formula:

[0107]

[0108] in: Indicates the components of the initial decomposition. Represents the amplitude modulation envelope. This represents the pure frequency-modulated components. For the separated amplitude modulation envelope, the steps from local extremum detection to component decomposition are repeated until the amplitude modulation envelope satisfies the monotonicity condition or the energy of the amplitude modulation envelope is lower than a preset percentage of the energy of the input signal frame, thus stopping the process. The resulting series of pure frequency-modulated components constitutes the set of eigenmode functions of the input signal frame. It can be understood that this iterative decomposition process aims to adaptively decompose complex non-stationary vibration signals into a series of physically meaningful narrowband components.

[0109] In some embodiments, an instantaneous frequency sequence is calculated for any pure frequency-modulated (PMMC) eigenmode function (EMF) component. A Hilbert transform is then performed on the PMMC component, which is achieved by convolving the PMMC component with an ideal 90-degree phase shifter to obtain orthogonal analytic signal components of the PMMC component. A corresponding complex analytic signal is constructed from the PMMC component and its analytic signal components. The real part of the complex analytic signal is the original PMMC component, and the imaginary part is its Hilbert transform result. The phase angle of the complex analytic signal at each sampling point is calculated. The phase angle is obtained by calculating the arctangent function of the four-quadrant ratio of the imaginary to real parts of the complex analytic signal. The phase angle between consecutive sampling points is differentially calculated to obtain the phase angle change between adjacent sampling points. Optionally, the phase angle change is divided by the product of twice pi and the sampling time interval to obtain the instantaneous frequency value of the EMF component at each sampling point in Hertz. Following the steps described above—Hilbert transform, constructing a complex analytic signal, calculating the phase angle, calculating the phase angle change using differential calculus, and calculating the instantaneous frequency value—each component in the intrinsic mode function set is processed sequentially to obtain the instantaneous frequency sequence of all components. It can be understood that the instantaneous frequency sequence describes the evolution of the frequency components of each intrinsic mode function over time.

[0110] In one embodiment of the present invention, the energy of each component in the intrinsic mode function set is calculated over the entire signal frame time. The energy of all intrinsic mode function components is normalized to obtain the energy proportion of each component. An energy proportion threshold is set, and intrinsic mode function components with energy proportions exceeding this threshold are selected as dominant energy components. For each dominant energy component, its instantaneous frequency sequence is statistically analyzed using a histogram distribution over the entire frequency range. The frequency intervals where the number of frequency points peaks are identified in the histogram distribution. The frequency intervals corresponding to each dominant energy component are merged, and overlapping intervals are combined to form a candidate fault feature frequency band set.

[0111] In specific implementation, the process of extracting frequency bands with energy proportions exceeding a preset threshold from the instantaneous frequency sequences of all intrinsic mode function (IMF) components as a candidate fault feature frequency band set is as follows: For an IMF set obtained after a signal frame undergoes an improved Hilbert-Huang transform, the IMF set contains M pure frequency-modulated (FM) IMF components. For each FM IMF component in the IMF set, its energy is calculated over the entire signal frame time. The energy calculation is completed by squaring the value of each sampling point of the FM IMF component and summing it over the entire time series. The total energy value of the signal frame in the decomposition domain is obtained by summing the energies of all M FM IMF components. Optionally, the energy of each FM IMF component is normalized by dividing the energy of each component by the total energy value to obtain the energy proportion of each component. The process of calculating the energy proportion of all components can be expressed in the following form:

[0112]

[0113] in: This represents the energy percentage of the k-th pure frequency-modulated intrinsic mode function component. This represents the energy of the k-th pure frequency-modulated eigenmode function component. This represents the sum of the energies of all pure frequency-modulated intrinsic mode function (PMEM) components. A specific energy percentage threshold is set, which is a preset constant between 0 and 1. In practice, all PMEM components with an energy percentage greater than the preset threshold are selected and marked as dominant energy components. It can be understood that the dominant energy components carry the main vibrational energy in the signal frame and are more closely associated with potential fault characteristics.

[0114] In some embodiments, for each pure frequency-modulated intrinsic mode function (EMF) component labeled as the dominant energy, the histogram distribution of its corresponding instantaneous frequency sequence across the entire frequency range is statistically analyzed. Statistical analysis of the instantaneous frequency histogram distribution first requires determining the frequency statistics range and bin width. The frequency statistics range covers the Nyquist frequency from 0 Hz to half the sampling frequency, and the bin width is set to a fixed frequency value. Each instantaneous frequency value in the instantaneous frequency sequence of the dominant energy EMF component is traversed, the frequency bin to which each instantaneous frequency value belongs is determined, and the count of that frequency bin is incremented by one, thus forming a histogram of the instantaneous frequency values ​​on the frequency axis. In the formed histogram distribution, the frequency intervals where local peaks occur in the frequency bin counts are identified. A local peak is defined as a frequency bin count that is simultaneously greater than the counts of both the preceding and following frequency bins. Optionally, for each dominant energy EMF component, the center frequency and bandwidth of the frequency intervals corresponding to all count peaks in its instantaneous frequency histogram can be recorded. The frequency intervals corresponding to each dominant energy pure frequency-modulated intrinsic mode function component are merged to process overlapping frequency intervals. This merging operation involves taking the union of intersecting frequency intervals to form a non-overlapping set of frequency bands. This set of frequency bands is the candidate fault feature frequency band set. It can be understood that extracting feature frequency bands from the instantaneous frequency distribution of the dominant energy components effectively focuses on the frequency bands with prominent energy in the signal, providing a target range for subsequent feature extraction.

[0115] See Figure 3 This is a graph showing the instantaneous frequency distribution of the dominant IMF components. It's a visualization result from the "candidate fault feature frequency band extraction" stage in power transformer winding fault detection, used to analyze the frequency energy concentration range of each dominant intrinsic mode function (IMF) after vibration signal decomposition. The graph shows the frequency distribution of five dominant energy components: IMF components 1, 2, 3, 4, and 8. Each component exhibits a single-peak or multi-peak distribution, with the peak positions corresponding to frequency ranges where fault features may be concentrated. These peak frequency ranges are the core source of subsequent candidate fault feature frequency bands; by merging overlapping ranges, characteristic frequency bands related to winding faults can be located. The count of non-peak regions is close to 0, indicating that the energy proportion of these frequency bands is extremely low and they do not belong to the dominant energy components, thus they can be excluded from fault feature analysis.

[0116] In one embodiment of the present invention, for each frequency band in the candidate fault feature frequency band set, the corresponding sub-band signal is extracted from the preprocessed vibration time-domain signal, and its higher-order statistical feature matrix is ​​calculated. The extraction process involves filtering the preprocessed vibration time-domain signal using a bandpass filter with a center frequency equal to the center frequency of the frequency band and a bandwidth equal to the width of the frequency band, to obtain the sub-band signal of that frequency band. The third moment, i.e., skewness, of the sub-band signal is calculated. The fourth moment, i.e., kurtosis, of the sub-band signal is calculated. The fifth moment, i.e., hyperkurtosis, of the sub-band signal is calculated. The sixth moment of the sub-band signal is calculated. From the sub-band signal, multiple segments of equal length are extracted, the squared amplitude of each segment is calculated, and the kurtosis of these squared amplitude values ​​is then calculated to obtain the amplitude squared kurtosis. From the sub-band signal, multiple segments of equal length are extracted, and the kurtosis of the envelope spectrum of each segment is calculated to obtain the envelope kurtosis. The calculated skewness, kurtosis, hyperkurtosis, sixth moment, amplitude squared kurtosis, and envelope kurtosis are arranged in rows to form a 6x1 feature matrix, which serves as the higher-order statistical feature matrix for this frequency band. A priori knowledge base containing vibration characteristic data collected and calculated from transformers of various models, capacities, and voltage levels under confirmed health conditions is loaded and stored. Based on the nameplate information of the currently tested transformer, including model, capacity, and voltage level, a matching query is performed in the priori knowledge base. A set of higher-order statistical feature matrices of all healthy samples that perfectly match the parameters of the current transformer is extracted from the knowledge base. The average value of all feature matrices in this set is calculated for each feature dimension. The average values ​​for each feature dimension are combined into a 6x1 average feature matrix, which serves as the health status reference feature matrix for the current detection task.

[0117] In practical implementation, the process of extracting the corresponding sub-band signal from the preprocessed vibration time-domain signal and calculating its higher-order statistical feature matrix for each frequency band in the candidate fault feature frequency band set, as well as loading the health state reference feature matrix, is as follows: The candidate fault feature frequency band set contains several specific frequency ranges defined by the center frequency and bandwidth, for example, a frequency band with a center frequency of 300 Hz and a bandwidth of 100 Hz. In practical implementation, a bandpass filter centered at the specific frequency band's center frequency and with a specific bandwidth is applied to the preprocessed vibration time-domain signal. The bandpass filter can be designed using a finite-length unit impulse response filter or an infinite-length unit impulse response filter. After filtering by the bandpass filter, a sub-band signal is obtained that retains only the frequency components within that frequency band. Optionally, the length of the sub-band signal is consistent with the length of the original preprocessed vibration time-domain signal, but its spectral energy is concentrated within the target frequency band. The third moment, i.e., the skewness, is calculated for the sub-band signal; the skewness characterizes the asymmetry of the sub-band signal's amplitude distribution. The fourth moment, or kurtosis, of the subband signal is calculated. Kurtosis characterizes the sharpness of the subband signal distribution curve. The fifth moment, or hyperkurtosis, of the subband signal is calculated. The sixth moment of the subband signal is calculated. Multiple equal-length, non-overlapping segments are extracted from the subband signal. The squared amplitude of all sampling points within each segment is calculated, and the kurtosis of these squared value sequences is then calculated to obtain the amplitude squared kurtosis. Multiple equal-length, non-overlapping segments are extracted from the subband signal. The envelope spectrum of each segment is calculated. The envelope spectrum is obtained by performing a Hilbert transform on the segment, taking the modulus, and then performing a Fourier transform. The kurtosis of the envelope spectrum is then calculated to obtain the envelope kurtosis. The calculation of these six eigenvalues ​​can be summarized in a feature list, as shown in Table 1.

[0118] Table 1: List of characteristics of higher-order statistics and their calculation descriptions

[0119]

[0120] In some embodiments, the calculation of the amplitude squared kurtosis involves segmenting the signal, the number of segments being denoted as... , No. The kurtosis of the squared amplitude sequence of a segment is calculated as follows:

[0121]

[0122] in: Indicates the first kurtosis of the squared-amplitude sequence Indicates the first The amplitude square sequence of the segment signal, This represents the mean of the squared magnitude sequence. This represents the mathematical expectation operator, and the final amplitude squared kurtosis eigenvalue is the average of all segmented kurtosis values. The calculated skewness, kurtosis, hyperkurtosis, sixth moment, amplitude squared kurtosis, and envelope kurtosis are arranged in rows to form a 6x1 feature matrix. This feature matrix is ​​the higher-order statistical feature matrix corresponding to the candidate fault characteristic frequency band. It can be understood that the higher-order statistical feature matrix characterizes the distribution characteristics of the sub-band signal from multiple statistical dimensions.

[0123] In practice, a health status reference feature matrix matching the current transformer model, capacity, and voltage level is loaded. A priori knowledge base, stored in a structured database format, is invoked, containing vibration characteristic data collected and calculated from transformers of various models, capacities, and voltage levels under confirmed health conditions. A precise matching query is performed in the priori knowledge base based on the model, rated capacity, and voltage level parameters recorded on the nameplate of the transformer under test. A set of high-order statistical feature matrices of all historical healthy samples that perfectly match the current transformer's model, capacity, and voltage level parameters is extracted from the priori knowledge base. The number of historical healthy samples is multiple. The average value of all feature matrices in the high-order statistical feature matrix set is calculated for each feature dimension. For example, for the skewness row, the arithmetic mean of the values ​​in the first row of all healthy sample feature matrices is calculated. Optionally, the average value calculation process involves aligning the corresponding rows of all healthy sample feature matrices and calculating the arithmetic mean row by row. The average values ​​calculated for each feature dimension are combined into a new six-row, one-column average feature matrix, which serves as the health status reference feature matrix for the current testing task. It is understandable that the health status reference feature matrix represents the vibration characteristic benchmark of the same type of transformer under fault-free conditions.

[0124] See Figure 4 This is a bar chart comparing the higher-order statistical features of transformer vibration signals under healthy and faulty states. It's used to identify anomalous feature dimensions in winding faults and represents a visualization result from the "abnormal feature submatrix identification" stage of fault detection. The sixth moment shows the largest difference, with the fault state feature value being 2.69 times that of the healthy state, making it the most sensitive fault feature dimension. Except for skewness, the other five higher-order statistical features show a significant increase under fault conditions, indicating a significant enhancement in the non-Gaussianity and impulse characteristics of the vibration signal, consistent with typical winding fault behavior. Skewness shows almost no change, indicating that this dimension is insensitive to the current fault mode and can be excluded from the abnormal feature submatrix. These significantly different feature dimensions will be extracted into the abnormal feature submatrix for subsequent spatial clustering and fault mode matching.

[0125] In one embodiment of the present invention, the high-order statistical feature matrix of a candidate frequency band to be detected is denoted as matrix D, and the health status reference feature matrix is ​​denoted as matrix H. Matrix D and matrix H have the same dimension. The absolute difference between the elements at corresponding positions of matrix D and matrix H is calculated. The absolute difference is divided by the absolute value of the element at the corresponding position of matrix H to obtain the relative difference of each feature dimension. The calculated relative difference is compared with a preset difference threshold. If the relative difference of a feature dimension exceeds the difference threshold, the feature dimension is marked as an abnormal feature dimension. All rows containing the marked abnormal feature dimensions are extracted from matrix D, and these row vectors are combined into a new submatrix, which is the abnormal feature submatrix. Each abnormal feature submatrix corresponds to one or more abnormal feature dimensions in a candidate fault feature frequency band. Each abnormal feature submatrix is ​​regarded as a feature point in a multi-dimensional space, and its feature dimension coordinates are composed of the row vector values ​​of the abnormal feature submatrix. The cosine similarity between any two abnormal feature submatrixes is calculated as a measure of their distance. A similarity threshold is set as a merging condition for clustering. Initialize all anomalous feature submatrices as independent clusters. Calculate the average cosine similarity between all pairs of clusters and find the pair with the highest average similarity. If the average similarity of this pair of clusters is greater than a similarity threshold, merge these two clusters into a new cluster and update the center feature point of the new cluster to the weighted average of the center feature points of the original two clusters. Repeat the steps of calculating the average cosine similarity between all pairs of clusters, finding the pair with the highest average similarity, and merging clusters when the condition is met, until no two clusters have an average similarity higher than the similarity threshold. Each remaining cluster is an independent suspected fault pattern cluster. For each final suspected fault pattern cluster, calculate the center feature vector of all its member anomalous feature submatrices. Extract the center frequency of the candidate fault feature frequency bands corresponding to all members of the suspected fault pattern cluster. Statistically analyze the distribution range of these center frequencies, find their minimum and maximum center frequencies, and use these to define the fault feature frequency coverage interval of the suspected fault pattern cluster. Match the center feature vector and the fault feature frequency coverage interval with a pre-set fault feature rule base. The fault feature rule base stores typical feature vector patterns and characteristic frequency ranges corresponding to different winding fault modes. The rule base is searched for fault mode entries that are closest to the central feature vector and whose fault feature frequency ranges overlap. The searched fault mode entries are then output as the diagnostic results for that suspected fault mode cluster.

[0126] In specific implementation, the process of comparing and analyzing the high-order statistical feature matrices of all candidate fault feature bands with the health status reference feature matrix element by element and identifying abnormal feature sub-matrices is as follows: The high-order statistical feature matrix of a candidate fault feature band to be detected is denoted as matrix D, and the health status reference feature matrix is ​​denoted as matrix H. Matrix D and matrix H have the same six-row, one-column dimension. Each row of matrix D and matrix H corresponds to the same high-order statistical feature; for example, the first row consists of skewness feature values. The absolute difference between the elements at corresponding positions in matrix D and matrix H is calculated. The absolute difference is calculated by subtracting the element value at the i-th row and 1-th column of matrix D from the element value at the i-th row and 1-th column of matrix H, and then taking the absolute value. The calculated absolute difference is divided by the absolute value of the element at the corresponding position in matrix H to obtain the relative difference for each feature dimension. The relative difference is a unitless ratio. Optionally, the relative difference can be expressed as:

[0127]

[0128] in: This represents the relative difference of the i-th feature dimension. This represents the value of the element in the i-th row and 1-th column of the matrix D to be detected. This represents the value of the element in the i-th row and 1-th column of the health reference matrix H. A preset difference threshold is set, for example, 0.3. The relative difference of each feature dimension is compared with this threshold. If the relative difference of a feature dimension is greater than the preset threshold, this feature dimension is marked as an abnormal feature dimension. The row vectors containing all the marked abnormal feature dimensions are extracted from matrix D. These row vectors are then combined in their original row order to form a new feature submatrix with fewer rows and the same number of columns. This feature submatrix is ​​the abnormal feature submatrix corresponding to the current candidate fault feature band.

[0129] In some embodiments, spatial clustering is performed on all identified anomalous feature submatrices. Each anomalous feature submatric corresponds to one or more anomalous feature dimensions in a candidate fault feature frequency band. For example, an anomalous feature submatric containing features from rows 2 and 5 corresponds to anomalies in kurtosis and amplitude squared kurtosis features in a specific frequency band. Each anomalous feature submatric is considered as a feature point in a multidimensional space, where the coordinate dimension of the feature point is formed by the sequential arrangement of the row vector values ​​contained in the anomalous feature submatric. The cosine similarity between any two anomalous feature submatrices is calculated as a measure of their distance. The cosine similarity is obtained by calculating the cosine of the angle between the vectors of the two anomalous feature submatrices; the closer the value is to 1, the closer the directions of the two vectors are. A similarity threshold is set as a merging condition for clustering; the similarity threshold is a preset constant between 0.5 and 1. All anomalous feature submatrices are initialized into independent clusters, each initially containing only one anomalous feature submatric member. The average cosine similarity between all pairs of members in all clusters is calculated. The average cosine similarity is obtained by calculating the arithmetic mean of the cosine similarities between all members of the anomalous feature submatrices in two clusters. Among all cluster pairs, find the pair with the highest mean cosine similarity. This implies that the pair with the highest mean cosine similarity has the most similar distribution patterns in the feature space.

[0130] In practice, if the average similarity of the pair of clusters with the highest average cosine similarity is greater than a preset similarity threshold, these two clusters are merged into a new cluster. After merging, the center feature point of the new cluster is updated to the weighted average of the center feature points of the original two clusters, and the weight can be determined by the number of anomalous feature sub-matrices contained in the original cluster. The steps of calculating the average cosine similarity between all pairs of clusters, finding the pair of clusters with the highest average similarity, and merging clusters when the condition is met are repeated until the average similarity between any two clusters is no higher than the preset similarity threshold. Each remaining independent cluster is a suspected fault mode cluster. For each final suspected fault mode cluster, the center eigenvector of all its member anomalous feature sub-matrices is calculated. The center eigenvector can be obtained by taking the arithmetic mean of the coordinates of all feature points within the cluster. The center frequencies of the candidate fault feature frequency bands corresponding to all member anomalous feature sub-matrices in the suspected fault mode cluster are extracted. The distribution range of these center frequencies is statistically analyzed, and the minimum and maximum center frequency values ​​are identified to define the fault feature frequency coverage interval of the suspected fault mode cluster. Optionally, the central feature vector and the fault feature frequency coverage interval are matched with a pre-set fault feature rule base. This rule base stores typical feature vector patterns and feature frequency intervals corresponding to different winding fault modes in a data structure. The rule base is searched for fault mode entries that have the closest Euclidean distance to the central feature vector and whose fault feature frequency coverage intervals overlap. The matched fault mode entries are output as the specific winding fault mode diagnosis result for that suspected fault mode cluster. The diagnosis result may be, for example, "winding axial clamping force loose" or "winding inter-turn short circuit." It can be understood that clustering and rule matching summarize scattered abnormal features into fault modes with clear physical meaning.

[0131] See Figure 5 This is a heatmap showing the matching of fault clusters with fault modes in a rule base, used in the final stage of transformer winding fault diagnosis. The similarity between the detected clusters and known fault modes is quantified by the matching distance. Cluster 1 (inter-turn short circuit) and Cluster 2 (winding deformation) both achieve perfect matches with their labeled fault modes, indicating a high degree of consistency between the clustering results and the rule base, making the diagnostic results reliable. The diagnostic results of Cluster 1 and Cluster 2 can be directly accepted, corresponding to inter-turn short circuit and winding deformation faults, respectively. Cluster 3 (insulation degradation) shows a significant deviation in its matching results, with a matching distance of 0 to the labeled "insulation degradation" mode, but the largest distance to the "loosening clamping force" mode, followed by the "inter-turn short circuit" mode. Excluding self-matching, its distance to the "winding deformation" mode is relatively small, suggesting that the cluster features may partially overlap with winding deformation, requiring further verification.

[0132] The above embodiments are only used to illustrate the technical methods of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical methods of the present invention without departing from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for detecting winding faults in power transformers based on vibration spectrum analysis, characterized in that, include: The original vibration acceleration signal of the transformer is acquired and preprocessed to obtain the preprocessed vibration time-domain signal, forming multiple signal frame sequences; For each of the signal frames, the improved Hilbert-Huang transform is applied to calculate the set of eigenmode functions of the signal frame; For each eigenmode function component in the set of eigenmode functions, calculate its instantaneous frequency sequence; From the instantaneous frequency sequence of all intrinsic mode function components, frequency bands with energy proportions exceeding a preset threshold are extracted as a set of candidate fault feature frequency bands; For each frequency band in the candidate fault feature frequency band set, the corresponding sub-band signal is extracted from the preprocessed vibration time domain signal, and its higher-order statistical feature matrix is ​​calculated. Load a health status reference feature matrix that matches the current transformer model, capacity, and voltage level; The high-order statistical feature matrix of all candidate fault feature bands is calculated and compared element by element with the health status reference feature matrix to identify abnormal feature sub-matrices with differences exceeding a preset threshold. Spatial clustering is performed on the abnormal feature sub-matrices to group the strongly correlated abnormal feature sub-matrices into one class, forming several independent suspected fault mode clusters. Based on the central features and frequency band distribution of each cluster, specific winding fault mode diagnosis results are output. The improved Hilbert-Huang transform includes: Local extremum detection is performed on the input signal frame to identify all maximum and minimum points; Construct the upper envelope using all maxima and the lower envelope using all minima. Calculate the mean of the upper and lower envelopes to obtain the original mean envelope; Subtracting the original mean envelope from the input signal frame yields a preliminary decomposed component; For the initially decomposed components, calculate their instantaneous amplitude at each sampling point, and perform extreme point detection and envelope construction on the instantaneous amplitude sequence again to obtain the amplitude modulation envelope of the components. The amplitude modulation envelope is separated from the initially decomposed components, thereby decomposing the original signal into the product of a pure frequency modulation component and an amplitude modulation component. The pure frequency-modulated component is output as an eigenmode function component that satisfies the condition. For the separated amplitude modulation envelope, the steps of detecting from local extrema and component decomposition are repeated until it satisfies monotonicity or a preset stopping criterion. The resulting series of pure frequency modulation components constitute the intrinsic mode function set. For each frequency band in the candidate fault feature frequency band set, the corresponding sub-band signal is extracted from the preprocessed vibration time-domain signal, and its higher-order statistical feature matrix is ​​calculated, including: The preprocessed vibration time-domain signal is filtered using a bandpass filter with a center frequency equal to the center frequency of the frequency band and a bandwidth equal to the width of the frequency band to obtain the sub-band signal of the frequency band. For the sub-band signal, calculate its third moment, i.e., skewness; For the sub-band signal, calculate its fourth moment, i.e., kurtosis; For the sub-band signal, calculate its fifth moment, i.e., the hyperkurtosis; For the sub-band signal, calculate its sixth moment; From the sub-band signal, extract multiple segments of equal length, calculate the amplitude square of each segment, and then calculate the kurtosis of the amplitude square to obtain the amplitude square kurtosis; From the subband signal, extract multiple segments of equal length, calculate the kurtosis of the envelope spectrum of each segment, and obtain the envelope kurtosis; The calculated skewness, kurtosis, hyperkurtosis, sixth moment, amplitude squared kurtosis, and envelope kurtosis are arranged in rows to form a six-row, one-column feature matrix, which serves as the higher-order statistical feature matrix of the frequency band.

2. The method for detecting power transformer winding faults based on vibration spectrum analysis according to claim 1, characterized in that, The original vibration acceleration signal of the acquired transformer is preprocessed to obtain a preprocessed vibration time-domain signal, forming multiple signal frame sequences, including: Vibration acceleration sensors are installed at designated locations on the outer wall of the oil tank of the power transformer to collect the original vibration acceleration signals of the transformer under steady-state operating conditions. The original vibration acceleration signal is preprocessed, including removing the DC component, bandpass filtering, and trend term elimination, to obtain the preprocessed vibration time-domain signal. The preprocessed vibration time-domain signal is divided into equal-length frames to form a continuous sequence of multiple signal frames, specifically including: Set a fixed time window length as the frame length; Set a window movement step size smaller than the frame length as the frame shift; Starting from the beginning of the preprocessed vibration time-domain signal, the first segment of the signal is extracted with a frame length as its width, and this segment is taken as the first signal frame. Move the capture window backward by one frame shift length, capture the second segment of the signal, and use it as the second signal frame; Repeat the window movement and capture operation until the end of the signal, thus forming a sequence of multiple signal frames that are continuous in time and partially overlap. For the last segment of a signal frame whose length is less than the frame length, zero values ​​are used to pad the end so that its length reaches the set frame length.

3. The method for detecting power transformer winding faults based on vibration spectrum analysis according to claim 2, characterized in that, For each eigenmode function component in the set of eigenmode functions, calculate its instantaneous frequency sequence, including: For any of the pure frequency modulation intrinsic mode function components, perform a Hilbert transform to obtain its orthogonal analytic signal components; A complex analytic signal is constructed from the eigenmode function components of the pure frequency modulation and its analytic signal components; Calculate the phase angle of the complex analytic signal at each sampling point; The phase angle between consecutive sampling points is differentially calculated to obtain the change in phase angle; Divide the change in phase angle by the product of twice pi and the sampling time interval to obtain the instantaneous frequency value of the intrinsic mode function component at each sampling point. Following the steps of Hilbert transform, constructing complex analytic signals, calculating phase angles, calculating phase angle changes using differential methods, and calculating instantaneous frequency values, each component in the intrinsic mode function set is processed sequentially to obtain the instantaneous frequency sequence of all components.

4. The power transformer winding fault detection method based on vibration spectrum analysis according to claim 3, characterized in that, From the instantaneous frequency sequence of all intrinsic mode function components, frequency bands with an energy proportion exceeding a preset threshold are extracted as a set of candidate fault feature frequency bands, including: For each component in the set of intrinsic mode functions, calculate its energy over the entire signal frame time; The energy of all intrinsic mode function components is normalized to obtain the energy percentage of each component. A threshold for energy percentage is set, and intrinsic mode function components whose energy percentage exceeds the threshold are selected as dominant energy components. For each dominant energy component, the histogram distribution of its instantaneous frequency sequence across the entire frequency range is statistically analyzed. In the histogram distribution, identify the frequency range where the number of frequency points reaches a peak; The frequency ranges corresponding to each dominant energy component are merged, and the overlapping ranges are combined to form the candidate fault feature frequency band set.

5. The power transformer winding fault detection method based on vibration spectrum analysis according to claim 4, characterized in that, Load a health status reference feature matrix that matches the current transformer model, capacity, and voltage level, including: Load a priori knowledge base containing vibration characteristic data collected and calculated from transformers of various models, capacities, and voltage levels under confirmed healthy conditions; Based on the nameplate information of the transformer being tested, including model, capacity, and voltage level, a matching query is performed in the prior knowledge base; Extract the set of high-order statistical feature matrices of all healthy samples that are completely consistent with the current transformer parameters from the knowledge base; Calculate the average value of all feature matrices in the set of higher-order statistics feature matrices for each feature dimension; The average values ​​of each feature dimension are combined into a six-row, one-column average feature matrix, which serves as the health status reference feature matrix for the current detection task.

6. The method for detecting power transformer winding faults based on vibration spectrum analysis according to claim 5, characterized in that, The high-order statistical feature matrix of all candidate fault feature bands calculated is compared element-by-element with the health status reference feature matrix, including: Let the high-order statistical feature matrix of a candidate frequency band to be detected be denoted as matrix D, and let the health status reference feature matrix be denoted as matrix H. Matrix D and matrix H have the same dimension. Calculate the absolute difference between the elements at corresponding positions of matrix D and matrix H; Divide the absolute difference by the absolute value of the element at the corresponding position in matrix H to obtain the relative difference of each feature dimension; Compare the calculated relative difference with the preset difference threshold; If the relative difference of a certain feature dimension exceeds the difference threshold, then the feature dimension is marked as an abnormal feature dimension. Extract the row vectors containing all the feature dimensions marked as abnormal from matrix D, and combine the row vectors into a new submatrix, which is the abnormal feature submatrix.

7. The method for detecting power transformer winding faults based on vibration spectrum analysis according to claim 6, characterized in that, Spatial clustering of the anomalous feature submatrix includes: Each abnormal feature submatrix corresponds to one or more abnormal feature dimensions in a candidate fault feature frequency band; Each anomalous feature submatrix is ​​considered as a feature point in a multidimensional space, and its feature dimension coordinates are composed of the row vector values ​​of the anomalous feature submatrix. Calculate the cosine similarity between any two anomalous feature submatrices as a measure of the distance between them; Set a similarity threshold as a condition for merging clusters; Initialize all anomalous feature submatrices as their own independent clusters; Calculate the average cosine similarity between all pairs of clusters, and find the pair of clusters with the highest average similarity; If the average similarity of the pair of clusters with the highest average similarity is greater than the similarity threshold, then the two clusters are merged into a new cluster, and the center feature point of the new cluster is updated to the weighted average of the center feature points of the original two clusters. Repeat the steps of calculating the average cosine similarity between all pairs of clusters, finding the pair of clusters with the highest average similarity, and merging clusters when the conditions are met, until the average similarity between any two clusters is higher than the aforementioned similarity threshold. Each remaining cluster is then considered an independent cluster representing a suspected fault mode.

8. The method for detecting power transformer winding faults based on vibration spectrum analysis according to claim 7, characterized in that, Based on the central characteristics and frequency band distribution of each cluster, specific winding fault mode diagnostic results are output, including: For each cluster of suspected fault modes that is ultimately formed, calculate the central eigenvector of all its member anomalous feature submatrices; Extract the center frequency of the candidate fault feature band corresponding to all members in the suspected fault mode cluster; The distribution range of the center frequency is statistically analyzed, and its minimum and maximum center frequencies are identified to define the fault characteristic frequency coverage range of the suspected fault mode cluster. The central feature vector and the fault feature frequency coverage interval are matched with a preset fault feature rule base. The fault feature rule base stores different winding fault modes, including inter-turn short circuit, winding deformation, insulation deterioration and loosening of clamping force, and their corresponding typical feature vector patterns and feature frequency ranges. Search the rule base for the fault mode entry that is closest to the central feature vector and whose fault feature frequency coverage range overlaps with the fault mode entry. The searched fault mode entries are output as the diagnostic results of the suspected fault mode cluster.

Images