Transformer fault prediction method based on time-space-spectrum multi-source information fusion

Abstract

The application relates to the technical field of transformer fault prediction, and discloses a transformer fault prediction method based on time-space-spectrum multi-source information fusion, which has the technical scheme points of data collection, data preprocessing, feature extraction, model construction and training, and fault prediction; through the fusion of time-space-spectrum multi-source information, the time-varying autoregressive moving average model (TV-ARMA) is adopted to dynamically capture the electrical parameter change, the adaptive wavelet packet transform (AWPT) is combined to optimize the frequency spectrum feature extraction, and the whale optimization algorithm (WOA) is used to train the deep belief network (DBN), so that the high-precision prediction of the transformer fault is realized, and the recognition accuracy and reliability of the winding short circuit, the core grounding and the insulation aging and the like are obviously improved.

Description

Technical Field

[0001] This invention relates to the field of transformer fault prediction technology, and more specifically, to a transformer fault prediction method based on the fusion of spatiotemporal and spectral multi-source information. Background Technology

[0002] Although multi-source information fusion technology has made some progress in the application of transformer fault prediction, existing multi-source information fusion technology has significant defects in data preprocessing, feature extraction, model optimization and judgment rules, making it difficult to meet the requirements of transformer fault prediction for high accuracy, strong generalization and dynamic adaptability.

[0003] Therefore, this invention provides a transformer fault prediction method based on the fusion of time-space-spectrum multi-source information, which improves the above-mentioned technical problems. Summary of the Invention

[0004] This invention aims to address the shortcomings of existing technologies by providing a transformer fault prediction method based on the fusion of spatiotemporal and spectral multi-source information. By fusing spatiotemporal and spectral multi-source information, this invention uses a time-varying autoregressive moving average model (TV-ARMA) to dynamically capture changes in electrical parameters, combines adaptive wavelet packet transform (AWPT) to optimize spectral feature extraction, and utilizes the whale optimization algorithm (WOA) to train a deep belief network (DBN). This achieves high-precision prediction of transformer faults and significantly improves the accuracy and reliability of identifying faults such as winding short circuits, core grounding, and insulation aging.

[0005] The above-mentioned technical objective of this invention is achieved through the following technical solution: a transformer fault prediction method based on time-space-spectral multi-source information fusion, comprising the following steps: S1, data acquisition, which includes: time information acquisition, spatial information acquisition, and spectrum information acquisition; S2, data preprocessing, which includes: normalization processing and noise reduction operation; S3, feature extraction, which includes: time feature extraction, spatial feature extraction, and spectrum feature extraction; S4, model construction and training, which includes: constructing a deep belief network (DBN) and using the whale optimization algorithm (WOA) to optimize the DBN weight parameters; S5, fault prediction, which includes: inputting test data into the trained DBN model, outputting the fault probability, and combining the threshold judgment rules of time, space, and spectrum features to comprehensively judge winding short circuit, core grounding, and insulation aging faults.

[0006] As a preferred technical solution of the present invention, the process of acquiring time information is as follows: Voltage signal V(t), current signal I(t), and power factor PF(t) during transformer operation are acquired in real time using high-precision sensors, and timestamps are recorded; a time-varying autoregressive moving average (TV-ARMA) model is used to dynamically model the voltage and current signals, and the power factor is calculated using an improved power factor formula; the process of acquiring spatial information is as follows: sensors are arranged on the transformer windings, core, and tank walls to acquire spatial distribution data of temperature, vibration, and partial discharge parameters; spatial distribution characteristics are analyzed using the spatial autocorrelation function I(T) and variogram γ(d); the probability density function is used to understand the spatial distribution probability of partial discharge and determine high-incidence areas of partial discharge; the process of acquiring spectrum information is as follows: spectrum analysis is performed on the neutral point current, winding current, and partial discharge signal; variational mode decomposition (VMD) combined with adaptive wavelet packet transform (AWPT) algorithm is used to extract spectrum features, including: frequency band energy proportion P. n,m and center frequency f n,m .

[0007] As a preferred embodiment of the present invention, the formula for calculating the spatial autocorrelation function I(T) is as follows:

[0008]

[0009] in, It is the average of all temperature data; spatial weight matrix d ij The sensor spacing is σ, and the weight attenuation parameter is σ. The TV-ARMA model solves for the time-varying coefficient α using the alternating direction multiplier method (ADMM). i (t) and c i (t), and updates it in real time to reflect dynamic changes caused by load mutations.

[0010] As a preferred technical solution of the present invention, the normalization process is as follows: the time, space and spectrum data are respectively normalized by Z-score and mapped to the [-1,1] interval; the denoising operation is as follows: a multi-scale morphological filtering algorithm is used to perform opening and closing operations on the data using structuring elements of different scales to remove noise.

[0011] As a preferred embodiment of the present invention, the multi-scale morphological filtering uses two-dimensional structuring elements to process spatial data, and the erosion and dilation calculation formulas are as follows:

[0012]

[0013] As a preferred embodiment of the present invention, the process of extracting the time features includes: extracting statistical features, autocorrelation coefficients, cross-correlation coefficients, and power spectral density of voltage, current, and power factor; the statistical features include: mean, variance, kurtosis, and skewness; the process of extracting the spatial features includes: calculating the temperature gradient. Spatial difference in vibration intensity σ A The partial discharge spatial location coordinates (x, y); the process of extracting the spectral features is as follows: extract the characteristic frequency f char Frequency band energy E i Center frequency Bandwidth B i And the spectral entropy value S.

[0014] As a preferred technical solution of the present invention, the characteristic frequency f char The search criteria are:

[0015] and And X(f)>τ

[0016] in, and These are the first and second derivatives of the spectral data X(f) with respect to frequency f, respectively, and τ is a pre-set threshold.

[0017] Frequency band energy E i The calculation formula is:

[0018]

[0019] Where, k i and k i+1 These are the frequency bands [f] i ,f i+1 The corresponding discrete frequency point index, where Δf is the frequency resolution;

[0020] Center frequency The calculation formula is:

[0021]

[0022] Bandwidth B i The calculation formula is: B i =f high -f low

[0023] The formula for calculating the spectral entropy S is:

[0024]

[0025] Where M is the number of discrete points in the spectrum data. This represents the proportion of energy at the i-th frequency point to the total energy.

[0026] As a preferred technical solution of the present invention, the process of constructing a deep belief network (DBN) is as follows: it consists of multiple stacked restricted Boltzmann machine (RBM) layers, the number of input layer nodes matches the number of features, and the output layer corresponds to the fault type or probability; the process of optimizing the DBN weight parameters using the whale optimization algorithm (WOA) includes: initializing the population, fitness evaluation, hunting behavior simulation and iterative optimization until convergence.

[0027] As a preferred embodiment of the present invention, the output layer of the DBN uses a softmax function for multi-class classification, as shown in the formula:

[0028]

[0029] Among them, o k It is the output value of the kth node in the output layer, representing the probability that the sample belongs to the kth type of fault; These are the weights connecting the last hidden layer unit j of the RBM to the output layer unit k. It is the bias of the output layer unit k, and C is the total number of fault types;

[0030] The parameter update formula for the Whale Optimization Algorithm (WOA) includes: Enclosure behavior: X i (t+1)=X * (t)-A·|C·X * (t)-X i (t)|; Spiral behavior: X i (t+1)=X * (t)·e bl ·cos(2πl)+X * (t); random search behavior: X i (t+1)=X rand -A·|C·X rand -X i (t)|.

[0031] As a preferred embodiment of the present invention, the fault determination rule in the fault prediction includes: winding short circuit fault: DBN output probability P w >0.7, and temperature gradient Frequency band energy E w >E w-thresh Core grounding fault: DBN output probability P c >0.6, and the power factor skewness |S PF |>S PF-thresh Spectral characteristic frequency amplitude Ach ar-c >A c-th resh Insulation aging fault: DBN output probability P a >0.5, and the partial discharge spectral entropy value S p >S p-th resh Power spectral density change rate

[0032] In summary, the present invention has the following beneficial effects:

[0033] Firstly, multi-dimensional information collaboration enhances diagnostic capabilities: A time-varying autoregressive moving average (TV-ARMA) model is employed to dynamically capture the transient characteristics of voltage and current, addressing the poor adaptability of traditional fixed-parameter models to load surges and improving the sensitivity of abnormal state detection. Spatial autocorrelation and variogram functions are used to quantitatively analyze temperature and vibration distribution patterns, accurately locating fault areas and avoiding misjudgments caused by sensor sparsity in traditional methods. Variational mode decomposition (VMD) and adaptive wavelet packet transform (AWPT) are used to extract features such as frequency band energy and spectral entropy, improving the accuracy of spectral analysis for non-stationary signals and effectively suppressing noise interference.

[0034] Secondly, efficient data preprocessing and feature fusion: Z-score normalization and multi-scale morphological filtering are used to eliminate heterogeneity of multi-source data, improve feature fusion efficiency, and accelerate model training convergence. By combining spatiotemporal and spectral statistical features (such as kurtosis, temperature gradient, and frequency band energy) with dynamic features (such as time-varying coefficients and autocorrelation coefficients), complex fault modes are comprehensively characterized, enhancing feature representation capabilities.

[0035] Third, intelligent model optimization and high-precision prediction: The Whale Optimization (WOA) algorithm is used to train a Deep Belief Network (DBN), improving parameter tuning efficiency compared to genetic algorithms, enhancing model generalization ability, and reducing false alarm rates under complex operating conditions such as load fluctuations and ambient temperature changes. Multi-dimensional threshold criteria (such as DBN output probability, temperature gradient, and spectral entropy) are combined to comprehensively determine the fault type. Attached Figure Description

[0036] Figure 1 The flowchart illustrates a transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion, as provided in an embodiment of the present invention. Detailed Implementation

[0037] This disclosure addresses the significant shortcomings of existing multi-source information fusion technologies in data preprocessing, feature extraction, model optimization, and decision rules, which make it difficult to meet the high accuracy, strong generalization, and dynamic adaptability requirements of transformer fault prediction. Therefore, this disclosure proposes a transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion. Through deep fusion of multi-source information, dynamic modeling optimization, and intelligent criterion design, it solves the technical bottlenecks of traditional methods, such as handling data heterogeneity, feature partiality, and insufficient model generalization ability, providing a highly reliable and low-cost transformer fault prediction solution for power systems.

[0038] Please refer to Figure 1 , Figure 1 A flowchart of the transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion according to an embodiment of this disclosure is shown. The overall process mainly includes the following five steps:

[0039] S1, Data Acquisition.

[0040] S1.1 Time Information Acquisition: High-precision sensors are used to collect voltage, current, and power factor during transformer operation in real time, and timestamps are recorded accurately.

[0041] In actual power systems, the changes of these electrical parameters over time contain a wealth of information and can reflect the operating status of transformers.

[0042] For the acquired voltage signal V(t) and current signal I(t), considering the harmonic interference and signal fluctuation characteristics present in the actual power grid, a time-varying autoregressive moving average (TV-ARMA) model is used to describe their variation over time. Traditional ARMA models assume that the model parameters are constant, but during transformer operation, due to load changes, environmental factors, and other influences, the electrical parameters exhibit time-varying characteristics. The expression for the TV-ARMA model is:

[0043]

[0044] Among them, a i (t), b j (t), c i (t), d j (t) are the time-varying autoregressive coefficients and moving average coefficients, respectively, which reflect the dynamic characteristics of the voltage and current signal changes over time t; p and q are the autoregressive order and moving average order, respectively; ∈(t) is a white noise sequence with a mean of 0, representing unpredictable random interference.

[0045] The power factor PF(t) can be calculated using the phase difference between voltage and current. Considering the errors in actual measurement and the time-varying characteristics of the signal, the following improved formula is adopted:

[0046]

[0047] in, It is the phase difference between voltage V(t) and current I(t); ΔV(t) and ΔI(t) are the differences between the voltage and current at the current moment and the previous moment, respectively, reflecting the rate of change of voltage and current; and These are the average values ​​of voltage and current over a period of time, respectively; k and m are weighting coefficients adjusted according to the actual situation, used to correct measurement errors and consider the impact of signal changes on the power factor.

[0048] These formulas allow for a more precise description of the changes in voltage, current, and power factor over time. Accurate timestamp recording is crucial when acquiring this time-series data. Employing a high-precision clock synchronization system to accurately capture minute changes in electrical parameters over short periods provides strong support for analyzing transient changes in transformer operating conditions.

[0049] When a transformer experiences a sudden load change, the time-varying coefficients α of voltage and current can be observed using the TV-ARMA model. i (t), c i (t) will change rapidly, reflecting the adjustment of the signal change pattern; at the same time, by using the improved power factor formula, the change of power factor during load changes can be calculated more accurately, providing a more accurate basis for judging the operating status of the transformer.

[0050] S1.2 Spatial Information Acquisition: Sensors are placed on the transformer windings, core, and tank walls to acquire spatial distribution information of temperature, vibration, and partial discharge parameters, and to determine the state differences at different locations.

[0051] To accurately describe the spatial distribution characteristics of these parameters, the concepts of spatial autocorrelation function (SAF) and variogram are introduced.

[0052] For temperature parameters, let's take data collected by temperature sensors located at different positions on the transformer as an example. Let the coordinates of the temperature sensors be (x...). i ,y i ), i = 1, 2, ..., n, the temperature value at the corresponding location is T(x i ,y i The spatial autocorrelation function is used to measure the spatial correlation of temperature, and its formula is:

[0053]

[0054] in, It is the average of all temperature data; spatial weight matrix It defines the position (x) i ,y i ) and (x j ,y j The spatial relationship between ) is given by σ, which is a parameter controlling the rate of spatial weight decay. The distance between two points is represented by I(T). The value of I(T) ranges from -1 to 1. A value closer to 1 indicates a stronger positive correlation in temperature distribution, meaning that temperature values ​​at adjacent locations are more similar; a value closer to -1 indicates a stronger negative correlation; and a value close to 0 indicates no spatial correlation. By calculating the spatial autocorrelation function, it is possible to determine whether there is a clustering phenomenon in the temperature distribution of a transformer. For example, if the I(T) value is high in a certain area, it indicates that the temperature change in that area is relatively uniform, which may indicate a risk of localized overheating.

[0055] The variogram is used to describe the degree of temperature variation with spatial distance, and the formula is:

[0056]

[0057] Where d is the spatial lag distance; N(d) is the number of sample point pairs with a distance of d; (x i+d ,y i+d ) is with (x i ,y i The location is defined by a distance *d*. The variogram γ(d) changes with *d*. When *d* is small, if *γ(d)* increases slowly, it indicates that the temperature changes little over a short distance, indicating good spatial continuity. Conversely, if *γ(d)* increases rapidly, it indicates that the temperature changes significantly over a short distance, potentially indicating a localized temperature anomaly. By analyzing the variogram curve, the spatial scale of temperature changes can be determined, i.e., the spatial range within which temperatures are similar or different. This is crucial for locating fault areas.

[0058] A similar method can be used for vibration parameters. Let the vibration sensor be located at different positions (x... m ,y m The collected vibration amplitude is A(x) m ,y m To account for the propagation characteristics of vibrations in different directions, an anisotropic variability function is introduced:

[0059]

[0060] Where θ represents direction; (x i+d,θ ,y i+d,θ ) is with (x i ,y i The position is at a distance d and in the direction θ; N θ (d) represents the number of sample point pairs at a distance d in direction θ. By calculating the anisotropic variability function in different directions, the propagation law and attenuation of vibration in different directions can be analyzed. For example, if γ in a certain direction... θ (d) Rapid growth indicates that the vibration changes significantly in that direction, which may indicate the presence of structural defects or faults related to that direction.

[0061] For partial discharge parameters, due to their strong randomness and locality, a spatial distribution based on probability density is used to describe them. Let the partial discharge sensor be located at position (x... k ,y k The detected partial discharge quantity is Q(x) k ,y k Define the probability density function p(Q) for partial discharge:

[0062]

[0063] Where δ is the Dirac function, when Q = Q(x) k ,y k When ), δ(QQ(x) k ,y k If the probability density function is 1, then the probability density function is 0; otherwise, it is 0. n is the total number of detected partial discharge events. The probability density function provides an intuitive understanding of the spatial distribution probability of partial discharge, allowing for the identification of high-incidence areas. Furthermore, by combining spatial autocorrelation and variogram analysis methods, the spatial correlation and variation patterns of partial discharge can be further studied, providing a basis for determining the type and severity of partial discharge faults.

[0064] These formulas and methods allow for in-depth analysis of the spatial distribution of transformer temperature, vibration, and partial discharge parameters, providing a more comprehensive understanding of the transformer's operating status from a spatial perspective and offering richer evidence for fault prediction.

[0065] S1.3 Spectrum Information Acquisition: Acquire transformer neutral point current, winding current and partial discharge signals and perform spectrum analysis. Use variational mode decomposition combined with adaptive wavelet packet transform algorithm to obtain the spectrum characteristics of the signal in different frequency bands.

[0066] First, variational mode decomposition (VMD) is performed on the acquired signal x(t). VMD decomposes the signal into K modal components u. k (t), each modal component has its corresponding center frequency ω.k Its core is solving a variational problem, the Lagrange expression of which is:

[0067]

[0068] Where α is a quadratic penalty factor used to balance signal reconstruction accuracy and modal bandwidth; λ(t) is a Lagrange multiplier used to ensure that the constraints are met; δ(t) is the Dirac function; * denotes convolution operation; and ‖·‖2 denotes the L2 norm.

[0069] Solving the above variational problem using the Alternating Direction Multiplier Method (ADMM) (optimizing the time-varying coefficients through the ADMM algorithm) yields the various modal components u. k (t). The original signal x(t) is decomposed into K modal components with different center frequencies, thus achieving preliminary separation of different frequency components of the signal.

[0070] For each modal component u k (t) Perform Adaptive Wavelet Packet Transform (AWPT). In the wavelet packet transform, the signal is recursively decomposed and reconstructed. Let... Let W be the scaling function, ψ(t) be the wavelet function, and W be the wavelet packet function for the nth layer and the mth node. n,m (t), we have:

[0071]

[0072] Where h(l) and g(l) are the coefficients of the low-pass filter and the high-pass filter, respectively, and Z represents the set of integers.

[0073] In adaptive wavelet packet transform, the optimal wavelet packet decomposition tree structure is selected based on the characteristics of the signal. Here, the information entropy criterion is introduced to select the optimal node. For a signal segment y(t) corresponding to a node, its information entropy E is defined as:

[0074]

[0075] During wavelet packet decomposition, the information entropy of each node is calculated, and the node with the lowest information entropy is selected for further decomposition until a preset decomposition depth is reached or other stopping conditions are met. This allows the wavelet packet decomposition results to better match the frequency characteristics of the signal and more accurately capture the signal's features in different frequency bands.

[0076] After adaptive wavelet packet transform, the wavelet packet coefficients c for each frequency band are obtained. n,m To extract more representative spectral features, the energy percentage P of each frequency band was calculated. n,m The formula is:

[0077]

[0078] Energy percentage P n,m This reflects the proportion of each frequency band in the total signal energy, highlighting the main frequency components of the signal. Simultaneously, the center frequency f of each frequency band is calculated. n,m :

[0079]

[0080] Where f(i) is the frequency value corresponding to the frequency sampling point. The center frequency f n,m It can more accurately reflect the frequency position of each frequency band, which helps to analyze changes in the frequency components of the signal.

[0081] The spectral features, such as energy proportion and center frequency, obtained by using the variational mode decomposition combined with the adaptive wavelet packet transform algorithm and related formulas can more comprehensively and accurately describe the characteristics of transformer neutral point current, winding current, and partial discharge signals in different frequency bands. When these spectral features are fused with temporal and spatial information, they can provide richer and more accurate data for transformer fault prediction, improving the accuracy and reliability of fault prediction.

[0082] S2, Data Preprocessing.

[0083] S2.1 Normalization Processing: The Zscore normalization method is adopted to uniformly map time, space, and spectrum data with different ranges and dimensions to the [-1,1] interval, thereby eliminating the influence of dimensions.

[0084] For the voltage, current, and power factor data in the time information, let the original data be X = {x1, x2, ..., x...} n}, its mean is The standard deviation is σ x The Z-score normalization formula is:

[0085]

[0086] Among them, z i This represents the normalized data value. This formula standardizes the original data, giving it zero mean and unit variance. For z... i The transformation is performed to map the data to the interval [-1, 1]. The transformation formula is as follows:

[0087]

[0088] Here, min(z) and max(z) are the minimum and maximum values ​​of the data z after initial Z-score normalization, respectively. iThis refers to the time information data ultimately mapped to the [-1, 1] interval. This processing method not only eliminates the influence of dimensions on time information data but also places different data on the same comparable scale, facilitating subsequent analysis. For example, voltage data may vary between hundreds of volts and thousands of volts, while current data may vary between a few amperes and hundreds of amperes. After this normalization process, they become numerically comparable and can better participate in subsequent model calculations.

[0089] Similar processing is applied to temperature, vibration, and partial discharge parameter data in spatial information. Taking temperature data as an example, let the original temperature data set be T={t1,t2,…,t…} m First, calculate its mean. and standard deviation σ t :

[0090]

[0091] Then, preliminary Z-score normalization is performed:

[0092]

[0093] Then map it to the interval [-1, 1]:

[0094]

[0095] in, This is the temperature value after preliminary Z-score normalization, min(t) z ) and max(t z These are the minimum and maximum values ​​of the temperature data after preliminary normalization. This results in temperature data that is ultimately normalized to the [-1,1] interval. By unifying spatially distributed temperature data to the same scale in this way, temperature data from different locations can be directly compared and analyzed, which is beneficial for uncovering the spatial characteristics of the temperature data.

[0096] For spectral information, suppose the wavelet packet coefficient data for a certain frequency band obtained after variational mode decomposition and adaptive wavelet packet transform is C = {c1, c2, ..., c k First, perform Z-score normalization on it:

[0097]

[0098] in, It is the mean of the wavelet packet coefficients, σ c The standard deviation is calculated using the following formulas:

[0099]

[0100] Next, to map it to the interval [-1, 1], a transformation is used:

[0101]

[0102] here, These are the wavelet packet coefficient values ​​after initial Z-score normalization, min(c z ) and max(c z These are the minimum and maximum values ​​of the wavelet packet coefficients after initial normalization. This results in spectral information data that is ultimately normalized to the [-1,1] interval. This processing allows spectral characteristic data from different frequency bands to be expressed on the same scale, enhancing the comparability between spectral features and facilitating a more accurate analysis of the relationship between spectral information and transformer faults.

[0103] Through the detailed Z-score normalization and mapping process described above, the dimensional differences between time, space, and spectral data can be effectively eliminated, providing a standardized and comparable data foundation for subsequent transformer fault prediction models based on multi-source information fusion, thereby improving the model's performance and the accuracy of fault prediction.

[0104] S2.2 Denoising Operation: A multi-scale morphological filtering algorithm is used to remove noise and interference from the data, thereby improving the purity of the data.

[0105] Multi-scale morphological filtering, based on the principles of mathematical morphology, processes data using structuring elements at different scales to effectively remove noise at various scales. Mathematical morphology primarily employs two basic operations: erosion and dilation, from which opening and closing operations are derived.

[0106] For a discrete one-dimensional data sequence x(n), n = 1, 2, ..., N (taking voltage data from time information as an example; other data processing methods are similar), the formula for the erosion operation is:

[0107]

[0108] Where b(m) is the structuring element and Z is the set of integers. The erosion operation can remove tiny bumps in the signal and suppress positive impulse noise.

[0109] The expansion formula is:

[0110]

[0111] The expansion operation, the opposite of erosion, is used to fill tiny voids in a signal and suppress negative impulse noise.

[0112] The opening operation involves erosion followed by dilation, and its formula is:

[0113]

[0114] Opening operations can remove isolated noise points and minor interference from a signal, smooth the signal profile, and preserve the main characteristics of the signal.

[0115] The closing operation involves dilation followed by erosion, and the formula is as follows:

[0116]

[0117] Closing operations can fill tiny gaps in a signal and connect broken parts, thus smoothing the signal.

[0118] In multi-scale morphological filtering, multiple structuring elements of different scales are used to process the data. Let the scale of the structuring element be s, and for structuring elements b of different scales... s (m), performs multi-scale opening and closing operations on the data x(n). For example, first perform small-scale opening operations on the data:

[0119]

[0120] Where s1 represents a smaller scale. This is the data after undergoing small-scale opening. Small-scale opening can remove small-scale noise and detail interference, while preserving the basic characteristics of the signal.

[0121] Then perform a large-scale closing operation:

[0122]

[0123] Here, s2 > s1, where s2 represents a larger scale. This is data after large-scale closing operations. Large-scale closing operations can further smooth the signal, remove larger-scale noise and interference, while preserving the overall trend of the signal.

[0124] By combining multi-scale opening and closing operations, noise at different scales in the data can be comprehensively processed. Small-scale operations focus on removing detailed noise, while large-scale operations are used to handle more macroscopic interference, thereby improving the purity of the data.

[0125] For two-dimensional spatial information data, such as temperature distribution data T(x) at different locations on a transformer... i ,y i The structural element is transformed into a two-dimensional form b(x). m ,y m The corresponding formulas for corrosion, expansion, opening, and closing operations are extended in two-dimensional space. For example, the two-dimensional corrosion operation formula is:

[0126]

[0127] Here, S is the region covered by the structuring element b. The formulas for two-dimensional expansion, opening, and closing operations are also extended according to similar logic to adapt to the processing needs of two-dimensional data:

[0128]

[0129] By processing spatial data using two-dimensional structural elements of different scales, noise and outliers in the spatial distribution are removed, thus better showcasing the spatial characteristics of the data.

[0130] For spectral information, since its data characteristics differ from those of temporal and spatial information, when performing multi-scale morphological filtering, the spectral data is treated as a discrete frequency sequence X(f). j For j = 1, 2, ..., M, the erosion, expansion, opening, and closing operations formulas described above are also applied, but at this time the structuring element b(f) k The design of the structuring element needs to be adjusted according to the frequency characteristics of the spectrum. For example, a Gaussian structuring element can be designed, whose expression in the frequency domain is:

[0131]

[0132] In this context, σ controls the scale of the structuring element; the larger σ is, the wider the coverage of the structuring element in the frequency domain. By using structuring elements of different scales to perform multi-scale morphological filtering on spectral data, noise bands in the spectrum can be removed, highlighting useful spectral features and providing more accurate data for subsequent spectral feature-based analysis.

[0133] By using the multi-scale morphological filtering algorithm described above, time, space, and spectral information data are processed in a targeted manner using structural elements of different scales. This effectively removes various types of noise and interference, improves data quality, and provides more reliable data support for subsequent transformer fault prediction.

[0134] S3, Feature Extraction.

[0135] S3.1 Time Feature Extraction: Extract statistical features, correlation features, and frequency domain features from time series data. Here, statistical features include mean, variance, kurtosis, and skewness; correlation features include autocorrelation coefficient and cross-correlation coefficient; and frequency domain features include power spectral density, thereby obtaining the temporal variation pattern of the data.

[0136] 1) Statistical characteristics

[0137] <1> Mean: For a voltage sequence V(n), n=1,2,…,N, its mean is... The calculation formula is:

[0138]

[0139] The average voltage reflects the average voltage level during transformer operation and can be used as a benchmark value for measuring voltage stability. For example, if the average voltage deviates significantly from the rated voltage over a period of time, it may indicate an abnormality in transformer operation.

[0140] For a current sequence I(n), the mean for:

[0141]

[0142] The average current reflects the average magnitude of the transformer's load current, which helps to understand the transformer's load conditions.

[0143] Mean of the power factor sequence PF(n)

[0144]

[0145] The average power factor reflects the average level of power utilization efficiency of a transformer and is of great significance for evaluating the economic efficiency of transformer operation.

[0146] <2> Variance: Voltage variance The formula used to measure the degree of voltage fluctuation around its mean is:

[0147]

[0148] The larger the voltage variance, the more severe the voltage fluctuation, which may cause greater impact on the insulation and other components of the transformer, increasing the risk of failure.

[0149] Current variance

[0150]

[0151] Current variance reflects the stability of the current. Abnormal current fluctuations may be related to factors such as internal winding faults in the transformer or sudden load changes.

[0152] Power due to variance

[0153]

[0154] The variance of the power factor reflects its stability. An unstable power factor can affect the power transmission efficiency of the power grid and may be a signal of problems in the internal or external circuits of the transformer.

[0155] <3> Kurtosis: Voltage kurtosis K VThe formula used to describe the peaks in the voltage distribution is:

[0156]

[0157] If the voltage kurtosis value deviates significantly from the kurtosis value of the normal distribution (usually 3), it may indicate the presence of abnormal spikes or pulse signals in the voltage. This may be due to harmonic interference in the power grid, partial discharge inside the transformer, or other reasons.

[0158] Current kurtosis K I :

[0159]

[0160] Abnormal changes in current kurtosis can reflect abnormalities in the current signal, such as sudden current changes or an increase in harmonic components, which may be related to transformer faults.

[0161] Power factor kurtosis K PF :

[0162]

[0163] Changes in power factor kurtosis can help determine the stability of the power factor and whether there are abnormal fluctuations. Abnormal kurtosis values ​​may indicate abnormal changes in the power conversion efficiency of the transformer.

[0164] <4> Deflection: Voltage deflection S V The formula used to measure the degree of asymmetry in voltage distribution is:

[0165]

[0166] Voltage skewness can reflect whether there are asymmetrical fluctuations in the voltage signal. For example, if the voltage skewness is not zero, it may indicate that there is a certain DC component or other asymmetrical interference factors in the voltage.

[0167] Current deflection S I :

[0168]

[0169] Current skewness helps to identify asymmetric characteristics in current distribution, which may be related to factors such as transformer winding asymmetry and load imbalance.

[0170] Power factor skewness S PF :

[0171]

[0172] Changes in power factor skewness can reflect the asymmetry of power factor distribution, which is of reference value for analyzing the power conversion characteristics of transformers and potential faults.

[0173] 2) Correlation characteristics

[0174] <1> Autocorrelation coefficient: Voltage autocorrelation coefficient r V (k) is used to measure the correlation between the voltage sequence and itself at different time delays k, and the formula is:

[0175]

[0176] By analyzing the change of voltage autocorrelation coefficient with time delay k, we can understand the periodicity and memory of voltage fluctuations. For example, if the voltage autocorrelation coefficient is high at a specific time delay k, it indicates that the voltage has a certain similarity and periodicity within that time interval, which may be related to the periodic load changes of the power grid or the inherent characteristics of the transformer itself.

[0177] Current autocorrelation coefficient r I (k):

[0178]

[0179] The autocorrelation coefficient of current can reflect the time correlation of current signals and help to determine the pattern and stability of current fluctuations.

[0180] Power factor autocorrelation coefficient r PF (k):

[0181]

[0182] The autocorrelation coefficient of the power factor helps to analyze the changing trend and stability of the power factor, and has certain reference value for predicting future changes in the power factor.

[0183] <2> Cross-correlation coefficient: the cross-correlation coefficient r between voltage and current VI (k) is used to measure the correlation between voltage V(n) and current I(n), and the formula is:

[0184]

[0185] Cross-correlation coefficients can help analyze the phase relationship and mutual influence between voltage and current. For example, under normal operating conditions, there is a certain phase difference between voltage and current, and cross-correlation coefficients can quantify this relationship. When cross-correlation coefficients show abnormal changes, it may indicate that the load characteristics of the transformer have changed or that there is a fault.

[0186] Cross-correlation coefficient r between voltage and power factor VPF (k):

[0187]

[0188] This cross-correlation coefficient can reflect the impact of voltage changes on the power factor and the potential relationship between the two, which is of great significance for analyzing the power conversion efficiency and operating status of transformers.

[0189] Cross-correlation coefficient r between current and power factor IPF (k):

[0190]

[0191] By analyzing this cross-correlation coefficient, we can understand the relationship between current change and power factor, and determine whether the load characteristics and power conversion efficiency of the transformer are normal.

[0192] 3) Frequency domain characteristics - power spectral density: voltage power spectral density P V (f) This describes the power distribution of the voltage signal at different frequencies, calculated using the Welch method, with the following formula:

[0193]

[0194] Where M is the number of segments, V i (f) is the Fourier transform of the i-th segment of the voltage sequence V(n). The voltage power spectral density can clearly show the energy distribution of the voltage in different frequency bands. If an abnormal peak appears in the power spectral density in a specific frequency band, it may indicate the existence of a voltage problem related to that frequency, such as harmonic interference.

[0195] Current power spectral density P I (f):

[0196]

[0197] Among them I i (f) is the Fourier transform of the i-th segment of the current sequence I(n). The current power spectral density helps to discover the harmonic components and frequency characteristics of the current, and plays an important role in judging the load characteristics and operating status of the transformer.

[0198] Power factor and power spectral density P PF (f):

[0199]

[0200] Where PF i(f) is the Fourier transform of the i-th segment of the power factor sequence PF(n). The power factor power spectral density can reflect the energy distribution of the power factor at different frequencies, helping to analyze the frequency characteristics and stability of the power factor, and providing valuable information for transformer operation monitoring and fault prediction.

[0201] By using the specific calculation formulas for voltage, current, and power factor data, we can more comprehensively and deeply explore the time variation patterns in time series data, providing strong time-dimensional information support for transformer fault prediction.

[0202] S3.2 Spatial Feature Extraction: Based on the spatially distributed sensor data, calculate the spatial differences in temperature gradient and vibration intensity, and the spatial location of partial discharge, thereby obtaining the spatial distribution characteristics that reflect the internal state of the transformer.

[0203] 1) Temperature gradient calculation: The temperature gradient is used to measure the drastic degree of temperature change between different locations on the transformer, reflecting the direction and rate of heat transfer. For transformers located at position (x... i ,y i ) and (x j ,y j The two temperature sensors have temperature values ​​T(x) and x(x) respectively. i ,y i ) and T(x j ,y j Temperature gradient vector The formula for calculating components in two-dimensional space is:

[0204] In the x-direction: (when x j -x i (≠0 o'clock)

[0205] In the y-direction: (when y j -y i (≠0 o'clock)

[0206] Combining the x and y directions, the magnitude of the temperature gradient is:

[0207] Temperature gradient direction: (The gradient direction is determined by the arctangent function)

[0208] Calculating temperature gradients helps identify areas of concentrated heat and heat propagation paths within a transformer. For example, a large temperature gradient in a transformer winding indicates drastic temperature changes in that area, potentially suggesting poor heat dissipation or localized overheating. This could accelerate insulation aging and increase the risk of failure.

[0209] 2) Calculation of Spatial Difference in Vibration Intensity: The spatial difference in vibration intensity can be measured by calculating the standard deviation of the vibration amplitude at different locations. Let the vibration sensor be located at n different locations (x... m ,y m The collected vibration amplitude is A(x) m ,y m ), m = 1, 2, ..., n. First, calculate the mean value of the vibration amplitude at all positions.

[0210]

[0211] Then calculate the spatial variation (standard deviation) σ of the vibration intensity. A :

[0212]

[0213] Significant spatial differences in vibration intensity mean that the vibration conditions vary considerably across different locations of the transformer, potentially indicating internal structural inhomogeneity or localized defects. For example, if the spatial difference in vibration intensity in a certain area is significantly higher than in other areas, it may indicate the presence of loose components or structural damage in that area. This can affect the mechanical stability of the transformer, and long-term accumulation could lead to more serious failures.

[0214] 3) Spatial localization calculation of partial discharge: Spatial localization can be performed using multiple partial discharge sensors, employing the hyperbolic localization method. Assume there are three partial discharge sensors located at (x1, y1), (x2, y2), and (x3, y3), and they detect partial discharge signals at times t1, t2, and t3, respectively. Given the propagation speed of the partial discharge signal as v, the propagation distances are r1 = v(t1 - t0), r2 = v(t2 - t0), and r3 = v(t3 - t0) (where t0 is the onset time of the partial discharge, which can be estimated using signal processing methods).

[0215] According to the definition of a hyperbola, for two sensors, such as sensor 1 and sensor 2, we have:

[0216] |r2-r1|=2a 12 (2a 12 (where the length of the hyperbola is the real axis)

[0217] |r3-r1|=2a 13 (2a 13 (The real axis length of another hyperbola)

[0218] By establishing the hyperbola equation:

[0219]

[0220] Solving the two hyperbolic equations simultaneously yields the coordinates (x, y) of the partial discharge source. In practical calculations, numerical methods such as iterative algorithms or least squares methods can be used to improve the location accuracy. Accurate spatial location of partial discharge is crucial for assessing the internal insulation condition of transformers and determining the fault location, helping to take timely and targeted measures to prevent further expansion of the fault.

[0221] By calculating the spatial differences in temperature gradient, vibration intensity, and spatial location of partial discharge, we can comprehensively obtain the spatial distribution characteristics of the transformer's internal state, providing key spatial dimension information for transformer fault prediction.

[0222] S3.3 Spectral Feature Extraction: Extract characteristic frequencies, band energy, center frequency, bandwidth, and spectral entropy from the spectral data to analyze the relationship between the spectral characteristics of the signal and the operating state of the transformer.

[0223] 1) Feature Frequency Extraction: Feature frequencies refer to frequency components in the spectrum that have special significance and can characterize the operating state of the transformer. By analyzing the spectrum of the transformer neutral point current, winding current, and partial discharge signal, a peak search algorithm is used to determine the feature frequencies. Let X(f) be the spectrum data of a certain frequency band obtained after variational mode decomposition and adaptive wavelet packet transform. Within the frequency range [f...], ... min ,f max Within ], define the characteristic frequency f ch ar The search criteria are:

[0224] and And X(f)>τ

[0225] in, and These are the first and second derivatives of the spectral data X(f) with respect to frequency f, used to determine the extreme points of the spectral curve; τ is a pre-set threshold used to filter out peak frequencies with large amplitudes and practical significance. The frequency f that satisfies the above conditions is the characteristic frequency f. ch ar .

[0226] For example, when a transformer core fault occurs, harmonics of specific frequencies are generated. These harmonic frequencies appear as high-amplitude peaks in the frequency spectrum. The aforementioned characteristic frequency extraction method can accurately identify these fault-related characteristic frequencies, providing an important basis for fault diagnosis.

[0227] 2) Bandwidth Energy Calculation: Bandwidth energy reflects the amount of energy contained in a signal within a specific frequency range, and is of great significance for analyzing the frequency component distribution of transformer signals. The entire spectrum is divided into multiple frequency bands [f] i ,f i+1], i = 1, 2, ..., N, frequency band energy E i The calculation formula is:

[0228]

[0229] In practical calculations, since spectral data is usually discrete, the trapezoidal integral method is used for approximate calculations.

[0230]

[0231] Where, k i and k i+1 These are the frequency bands [f] i ,f i+1 The corresponding discrete frequency point index is Δf, where Δf is the frequency resolution.

[0232] Different types of faults can cause changes in the energy distribution of signals across different frequency bands. For example, a winding short-circuit fault may significantly increase the energy in certain frequency bands. By monitoring changes in frequency band energy, abnormalities in the transformer's operating status can be detected in a timely manner.

[0233] 3) Center Frequency Calculation: The center frequency describes the central location of the frequency distribution within a frequency band and can more accurately reflect the frequency characteristics of the band. For each frequency band [f i ,f i+1 Center frequency The calculation formula is:

[0234]

[0235] Similarly, in the case of discrete data, numerical calculation methods are used:

[0236]

[0237] Changes in center frequency reflect shifts in signal frequency components. When a fault occurs inside a transformer, the frequency components of the signal change, leading to a change in the center frequency. By monitoring the trend of center frequency changes, it is possible to effectively determine whether the transformer is operating normally.

[0238] 4) Bandwidth Calculation: Bandwidth measures the breadth of a frequency band and reflects the degree of dispersion of signal frequency components. For a frequency band [f i ,f i+1 ], bandwidth B i The calculation can be performed using the half-power bandwidth method, which involves finding two frequencies f within the bandwidth where the amplitude is half of the maximum power amplitude. low and f high (satisfy f low ≤f highIf the bandwidth B is... i for:

[0239] B i =f high -f low

[0240] Variations in bandwidth are closely related to the type and severity of transformer faults. For example, partial discharge faults can broaden the signal bandwidth because partial discharge generates multiple frequency components, resulting in a more dispersed frequency distribution. Monitoring these bandwidth variations can provide valuable information for fault diagnosis.

[0241] 5) Calculation of Spectral Entropy: Spectral entropy is an indicator that measures spectral complexity and uncertainty, reflecting the uniformity of the distribution of signal frequency components. For spectral data X(f), the formula for calculating the spectral entropy S is based on information entropy theory and uses the following formula:

[0242]

[0243] Where M is the number of discrete points in the spectrum data. This represents the proportion of energy at the i-th frequency point to the total energy.

[0244] Under normal operating conditions, the spectral entropy value of a transformer is relatively stable. When a transformer malfunctions, the spectral composition becomes more complex, causing changes in the spectral entropy value. For example, during insulation aging, the spectral entropy value may gradually increase. By monitoring changes in the spectral entropy value, potential fault hazards can be detected in advance.

[0245] By extracting the aforementioned spectral features, namely characteristic frequency, band energy, center frequency, bandwidth, and spectral entropy, and combining them with temporal and spatial information for comprehensive analysis, a more comprehensive and in-depth understanding of the transformer's operating status can be achieved, providing strong support for transformer fault prediction.

[0246] S4. Model building and training.

[0247] S4.1 Constructing a Deep Belief Network (DBN): This consists of multiple Restricted Boltzmann Machine (RBM) layers stacked together, plus an output layer. The number of input layer nodes is determined by the number of extracted features, and the number of output layer nodes corresponds to different fault types or fault probability levels of the transformer.

[0248] A Restricted Boltzmann Machine (RBM) is an energy-based undirected graph model consisting of a visible layer and a hidden layer. Nodes within a layer are unconnected, while connections exist only between nodes in different layers. For a given array of nodes with size n... v visible units and n h The energy function of an RBM with hidden units is defined as:

[0249]

[0250] Among them, v i h represents the state of the i-th unit in the visible layer (usually 0 or 1). j This represents the state of the j-th unit in the hidden layer (its value is usually 0 or 1), a i It is the bias of visible layer unit i, b j It is the bias of hidden layer unit j, w ij It is the weight connecting visible layer unit i and hidden layer unit j.

[0251] Based on the energy function, the joint probability distribution of RBMs is:

[0252]

[0253] Where Z = ∑ v ∑ h exp(-E(v,h)) is the partition function, used to normalize the probability distribution. Since the exact calculation of the partition function is very difficult in high-dimensional cases, the contrastive divergence (CD) algorithm is usually used for approximate training.

[0254] In DBN, the hidden layer output of the previous RBM is used as the visible layer input of the next RBM. For the l-th RBM (l = 1, 2, ..., L-1, where L is the total number of RBM layers), let its visible layer input be v. (l) The hidden layer output is h (l) .

[0255] The activation probability from the visible layer to the hidden layer is calculated as follows:

[0256]

[0257] in, It is the sigmoid activation function. It is the number of nodes in the visible layer of the l-th RBM. and These are the bias of hidden layer unit j in the l-th RBM and the weight connecting visible layer unit i and hidden layer unit j, respectively. This formula represents the bias of hidden layer unit j in the l-th RBM and the weight connecting visible layer unit i and hidden layer unit j, respectively. (l) When the probability of hidden layer unit j being activated (with a value of 1) depends on the bias of the hidden layer unit, the connection weights, and the state of the visible layer unit.

[0258] The reconstruction probability from the hidden layer to the visible layer is calculated as follows:

[0259]

[0260] in, It is the bias of the l-th visible layer cell i of the RBM. This is the number of nodes in the l-th RBM hidden layer. This formula is used when the hidden layer state h is known. (l) In the case of , calculate the probability that visible layer cell i is reconstructed as 1.

[0261] During training, the parameters of the RBM (weights w and biases a, b) are adjusted by minimizing the reconstruction error. Taking the l-th RBM as an example, the contrastive divergence algorithm is used, and its goal is to minimize the following reconstruction error function:

[0262]

[0263] Where V is the set of all visible layer samples in the training dataset. The contrastive divergence algorithm approximates the gradient by performing k steps (usually k=1 or k=2) of Gibbs sampling on the training samples, and then updates the parameters. The specific update formula is as follows:

[0264]

[0265] Where ∈ is the learning rate, which controls the step size for parameter updates; <·> data This represents the expectation calculated based on the training data, <·> recon This represents the expectation calculated based on the reconstructed data. By continuously iterating and updating these parameters, the RBM can better learn the feature representation of the data.

[0266] After training all the RBMs layer by layer, the hidden layer output of the last RBM is connected to the output layer. Assuming the output layer is a fully connected layer, for multi-class problems (corresponding to different transformer fault types), the output layer uses the softmax activation function, calculated as follows:

[0267]

[0268] Among them, o k It is the output value of the kth node in the output layer, representing the probability that the sample belongs to the kth type of fault; These are the weights connecting the last hidden layer unit j of the RBM to the output layer unit k. is the bias of output layer unit k, and C is the total number of fault types.

[0269] For regression problems (corresponding to fault probability levels), the output layer can use a linear activation function, with the output value directly representing the predicted fault probability level. Through this network structure and training method, DBN can effectively learn the complex relationship between input data (data fused from temporal, spatial, and spectral features) and transformer fault types or fault probability levels, providing strong model support for fault prediction.

[0270] S4.2 employs the Whale Optimization Algorithm (WOA) to optimize the DBN weight parameters; the preprocessed and feature-extracted data is divided into training and testing sets. The DBN is trained using the training set, and during training, the Whale Optimization Algorithm (WOA) is used to optimize the DBN weight parameters. The specific steps are as follows:

[0271] 1) Initialization: Randomly generate a set of individual whales (solutions), each representing a set of weight parameters for DBN. Set the algorithm parameters, including population size, maximum number of iterations, and search space range.

[0272] Let the size of the whale population be N. In this patent, to ensure that the Whale Optimization Algorithm (WOA) can search in a diverse initial solution space and increase the probability of finding the global optimum, the weights and biases of each layer of the Restricted Boltzmann Machine (RBM) in the DBN are randomly initialized.

[0273] For the l-th layer RBM in DBN (l=1,2,…,L-1, where L is the total number of RBM layers), the weight matrix W connecting the visible layer and the hidden layer (l) elements ( The number of visible nodes in the l-th layer of the RBM. (where the number of hidden layer nodes is ), the initialization formula is:

[0274]

[0275] in, The reason for setting the search space range in this way is that the magnitude of the weight values ​​has a crucial impact on the training process and performance of the neural network. Smaller weight values ​​help prevent overfitting in the early stages of training, making neuron updates more stable and facilitating robust learning of data features. Conversely, larger weight values ​​may lead to excessively large neuron outputs, making it difficult for the model to converge during training. Setting the weight value range to [-1, 1] ensures sufficient diversity of weights, allowing the model to explore different parameter combinations, while also limiting the magnitude of weight changes to avoid adverse effects on model training from excessively large or small weight values. `rand()` is a random number generator function that uniformly distributes random numbers within the interval [0, 1]. The random numbers generated by this function randomly initialize the weights of each individual whale within the specified range, thus ensuring the diversity of the initial population and providing a rich initial solution space for WOA to find the global optimum in subsequent search processes.

[0276] For the bias vector a of the visible layer of the l-th RBM (l) elements and the bias vector b of the hidden layer (l) elements The initialization formulas are as follows:

[0277]

[0278] in, Bias in neural networks provide neurons with additional learnable parameters, helping the model better fit the data. Setting the bias search space to [-1, 1] avoids excessive influence from overly large or small bias values ​​on neuron output, ensuring that the model can reasonably adjust the bias during training, thereby optimizing model performance. Randomly generating bias values ​​within the appropriate range using the rand() function further increases the diversity of initial solutions, allowing the algorithm to search from different starting points in the early stages of training, increasing the probability of finding the global optimum.

[0279] The maximum number of iterations is set to T = 300. The maximum number of iterations is a key parameter in WOA (World of Arithmetic Analysis), determining the number of iterations the algorithm performs in the search for the optimal solution. If the maximum number of iterations is set too small, the algorithm may not be able to fully explore the solution space, making it difficult to find the global optimum or a parameter combination close to the global optimum, resulting in low model prediction accuracy. Because transformer fault prediction involves multi-source information fusion and complex data features, DBN (Deep Transformer Network) contains multiple RBM (Resource Block Model) layers, and the parameters of each layer need fine-tuning. Too few iterations will prevent the algorithm from fully exploring the potential patterns in the data. For example, if the number of iterations is only a few dozen, WOA may stop searching before fully exploring the complex relationship between data features and transformer faults, causing the trained DBN model to not fit the data well, thus affecting fault prediction.

[0280] 2) Fitness evaluation: For each individual whale, its corresponding weight parameters are applied to the DBN, training is performed using training data, and the mean squared error on the validation data is calculated. The reciprocal of the mean squared error is used as the fitness value. The higher the fitness value, the better the DBN performance corresponding to that set of weight parameters.

[0281] Let the training set be Where x m It is an input sample that incorporates temporal, spatial, and spectral features. m These are the corresponding real labels (which can be fault type codes or fault probability levels). The validation set is...

[0282] When a set of weight parameters representing individual whales is applied to DBN, DBN applies the weights to the input sample x. n The predicted output is Mean squared error (MSE) measures the degree of error between the predicted and the true value. Its calculation formula is as follows:

[0283]

[0284] This formula means that the difference between the predicted value and the true value for each sample in the validation set is squared, and then the average of these squared differences is calculated. The squaring operation is used to prevent positive and negative errors from canceling each other out, making the error measurement more reasonable. The smaller the mean squared error, the closer the model's predicted values ​​are to the true values, and the better the model's predictive performance.

[0285] The reciprocal of the mean squared error is used as the fitness value, i.e.:

[0286]

[0287] The reciprocal of the mean squared error (MSE) is used as the fitness value because, during optimization, the goal is to find the combination of weight parameters that minimizes the model's prediction error. The smaller the MSE, the larger its reciprocal; in the whale optimization algorithm, a larger fitness value represents a better solution. This method quantifies the model's predictive performance on the validation set into a fitness value, allowing the whale optimization algorithm to judge the quality of each individual whale (i.e., each set of weight parameters) based on the fitness value during iteration. This guides the algorithm towards finding better weight parameters, continuously improving the prediction accuracy of the DBN model and thus better serving the transformer fault prediction task.

[0288] 3) Hunting Behavior Simulation: Simulate the humpback whale's behavior of surrounding prey, spiraling upwards, and random searching, and update the position of individual whales (i.e., the weight parameters of DBN). By adjusting the position, the individual whales are gradually brought closer to the optimal solution.

[0289] In the whale optimization algorithm, humpback whales surround their prey when searching for it. Let X be the current position vector of the individual whale (corresponding to the weight parameter vector of DBN). i The optimal position vector (corresponding to the optimal combination of weight parameters) found in the current iteration is X. * First, we introduce the parameter 'a', which decreases linearly from 2 to 0 with the number of iterations 't', as shown in the formula:

[0290]

[0291] Where t is the current iteration number, and T is the set maximum iteration number. 'a' plays a crucial role in controlling the search range and convergence speed in the algorithm. As iterations proceed, 'a' gradually decreases, causing the search range of the whale to shrink, and the algorithm to gradually converge towards the optimal solution.

[0292] The formulas for calculating the coefficient vectors A and C are as follows:

[0293] A = 2a·r1-a

[0294] C = 2·r²

[0295] Here, r1 and r2 are random vectors within the interval [0,1], and are regenerated for each calculation. A and C are used to determine the method and direction of updating the whale's individual position. Vector A, combined with the value of a, allows the whale to search a wider range of values ​​and explore different solution spaces in the early stages of the algorithm, as a is larger. In the later stages, as a decreases, the range of A also shrinks, and the whale's search area gradually focuses on the vicinity of the optimal solution. Vector C, through its randomized value, increases the randomness of the search process, preventing the algorithm from getting trapped in local optima.

[0296] Based on the relationship between the absolute value of A and 1, whales exhibit two distinct behavioral patterns. When |A| < 1, the whale engages in a predator-surrounding behavior, and the formula for updating its position is:

[0297] X i (t+1)=X * (t)-A·|C·X * (t)-X i (t)|

[0298] This formula indicates that an individual whale moves closer to its current optimal position, and the degree and direction of this approach are determined by A, C, and the difference between the current individual's position and the optimal position. Using the current optimal position as a reference point, the whale adjusts its position based on the calculation results of A and C, gradually reducing the distance to the optimal solution.

[0299] When |A|≥1, the whale performs a random search, as shown in the formula:

[0300] X i (t+1)=X rand -A·|C·X rand -X i (t)|

[0301] Among them, X rand The position is randomly selected from the population of whales. In this case, the whale no longer searches around its current optimal position, but instead randomly selects a position in the solution space as a reference for position updating. This helps the algorithm to explore the solution space extensively in the early stages of the search, increasing the likelihood of finding the global optimum and preventing the algorithm from converging to a local optimum too early.

[0302] In addition to the two behaviors mentioned above, humpback whales also exhibit a spiral ascent approach to their prey. Let's introduce a parameter l, which follows a normal distribution. Used to simulate a spiral shape. The position update formula for the spiral ascent behavior is:

[0303] X i (t+1)=X * (t)·e bl ·cos(2πl)+X * (t)

[0304] Here, b is a constant used to control the shape of the spiral. In this patent, based on actual testing and experience, b = 1 is set. This formula indicates that the individual whale gradually approaches the current optimal position in a spiral manner. This behavior pattern combines the characteristics of global search and local search, maintaining a fine search of the area near the current optimal solution while also being able to break out of local optima to a certain extent and explore better solutions.

[0305] In the actual algorithm execution, for each individual whale, a certain probability p (p = 0.5 in this patent) determines whether to perform spiral ascent behavior or the aforementioned prey encirclement or random search behavior. Through this comprehensive hunting behavior simulation, the position of the individual whales (i.e., the weight parameters of DBN) is continuously updated, allowing the whale group to continuously explore and optimize in the solution space, gradually approaching the global optimum. This leads to finding a more suitable combination of DBN weight parameters for transformer fault prediction tasks, thereby improving the model's predictive performance.

[0306] 4) Iterative optimization: Repeat the fitness evaluation and hunting behavior simulation steps until the maximum number of iterations is reached or the termination condition is met to obtain the optimal weight parameters. During iterative optimization, in addition to repeating the fitness evaluation and hunting behavior simulation steps, a dynamic learning rate adjustment strategy and a population diversity maintenance mechanism are introduced to further improve algorithm performance.

[0307] In each iteration, dynamically adjusting the learning rate helps balance the algorithm's convergence speed and search accuracy. An excessively large learning rate may cause the algorithm to skip the optimal solution during the search process, while an excessively small learning rate will slow down the convergence speed. Therefore, the following dynamic learning rate adjustment formula is used:

[0308]

[0309] Where η(t) represents the learning rate at the t-th iteration, η0 is the initial learning rate, λ is a parameter controlling the rate of decay of the learning rate, and T is the maximum number of iterations. As the number of iterations t increases, the learning rate η(t) gradually decreases. A larger learning rate is maintained in the early stages of the algorithm to speed up the search, while the learning rate is reduced as the algorithm approaches the optimal solution, enabling it to converge to the optimal solution more accurately.

[0310] To maintain population diversity and prevent the algorithm from prematurely getting trapped in local optima, a diversity assessment method based on population entropy is introduced. The formula for calculating population entropy H is:

[0311]

[0312] Where N is the population size, p i It represents the proportion of the i-th individual whale in the population. When the population entropy H falls below a certain threshold H... thresh When the population diversity is insufficient, random perturbation is applied to some individual whales to increase population diversity. The formula for random perturbation is:

[0313]

[0314] in, X represents the position of the individual whale after perturbation (corresponding to the weight parameter vector of DBN). iThe original position is α, the disturbance intensity coefficient is X. max and X min These are the maximum and minimum values ​​in the search space, respectively, and randn() is a random number that follows a standard normal distribution.

[0315] During the iteration process, the iteration stops when the maximum number of iterations T is reached or a specific termination condition is met (such as the change of the optimal solution after multiple consecutive iterations being less than a certain minimum value ∈), and the optimal weight parameters are obtained.

[0316] S5. Fault Prediction: The test set data is input into the trained DBN model to obtain the transformer fault prediction results, thereby determining whether the transformer has a fault and the type or severity of the fault. The following prediction and judgment schemes are adopted for different types of faults:

[0317] (1) Determination of winding short circuit fault

[0318] 1) Probability based on DBN model output: When the DBN model outputs the probability P corresponding to a winding short-circuit fault. w Exceeding the set threshold P w-th resh At that time, it was initially determined that the transformer might have a winding short-circuit fault. w-th resh These are empirical values ​​derived from statistical analysis of a large amount of historical fault data and normal operation data, used to distinguish between normal operation and winding short-circuit fault states. For example, through analysis of long-term monitoring data from multiple transformers, it was found that when P... w When P is greater than 0.7, the probability of a winding short-circuit fault actually occurring is relatively high. In this case, P can be adjusted. w-th resh Set it to 0.7.

[0319] 2) Combining time characteristics: The time-varying autoregressive moving average (TV-ARMA) model is used to analyze the changes in the time-varying coefficients of voltage V(t) and current I(t). If, within a certain time period, the voltage autoregressive coefficient a... i (t) and current autoregressive coefficient c i The rate of change of (t) exceeds their respective thresholds S. a-th resh and S c-th resh This increases the basis for judging the possibility of winding short-circuit faults. The formulas for calculating the rate of change are as follows:

[0320]

[0321] Where Δt is the time interval. S a-th resh and S c-th resh It is obtained by statistical analysis of the rate of change of time-varying coefficients under normal operation and known winding short-circuit fault conditions, and is used to measure the degree of abnormality of coefficient changes.

[0322] 3) Combining spatial characteristics: Calculate the temperature gradient in the transformer winding region. like This indicates that the temperature in this area is fluctuating drastically, possibly due to localized overheating caused by a winding short circuit. g-th resh This threshold is determined based on historical data of the temperature gradient in the winding area during normal transformer operation; exceeding this value indicates an abnormal temperature distribution. Simultaneously, if the spatial difference (standard deviation) of vibration intensity σ... A Greater than the vibration threshold σ A-th resh This further supports the judgment of winding short-circuit faults. σ A-th resh It is based on statistical analysis of a large amount of transformer vibration data under normal operation and fault conditions, and is used to define the normal range of vibration differences.

[0323] 4) Combining spectral characteristics: Analyze the winding current spectrum. If the frequency band energy E of a specific frequency range (such as f1-f2) is... w Satisfy E w >E w-th resh And the center frequency f of this frequency band c Deviation from the normal range exceeding Δf c-th resh This indicates an abnormal winding current spectrum, related to a winding short-circuit fault. The formula for calculating bandwidth energy is: The formula for calculating the center frequency is:

[0324]

[0325] E w-th resh , Δf c-th resh The threshold is determined by spectral analysis of a large number of winding short-circuit fault samples and normal operation samples, and is used to judge the degree of abnormality of frequency band energy and center frequency.

[0326] (2) Determination of iron core grounding fault

[0327] 1) Probability output based on DBN model: If the probability P of core grounding fault in the DBN model prediction result is... c Greater than threshold P c-th resh Therefore, it is preliminarily determined that the transformer has a risk of core grounding fault. c-th resh It is determined based on historical fault data and normal operation data. For example, through the analysis of a large amount of transformer data, it was found that when P c When the value is greater than 0.6, the probability of an actual core grounding fault occurring is relatively high. Therefore, P... c-thresh Set it to 0.6.

[0328] 2) Combining time characteristics: Calculate the skewness S of the power factor PF(t). PF and kurtosis K PF If |S PF |>S PF-th resh And |KPF |>K PF-th resh This indicates an abnormal power factor distribution, which may be related to a core grounding fault. PF-thresh and K PF-th resh It is a threshold determined based on statistical data of the skewness and kurtosis of the power factor during normal operation, and is used to measure the degree of abnormality in the power factor distribution.

[0329] 3) Combining spatial characteristics: Analyze the spatial correlation of temperature in the transformer core region using the spatial autocorrelation function I(T). If in the core region, I(T) is greater than the threshold I0. T-thresh This indicates that the temperature correlation in this area is too strong, which may be due to localized heat accumulation caused by the grounding of the iron core. T-thresh It is obtained through statistical analysis of the spatial autocorrelation function values ​​of the core area temperature during normal operation and core grounding faults, and is used to determine abnormal conditions of temperature correlation.

[0330] 4) Combining spectral characteristics: If a characteristic frequency f is detected in the neutral point current spectrum of the transformer... ch ar-c (For example, around 100Hz, which is a characteristic frequency that often occurs in iron core grounding faults), and the amplitude A corresponding to this characteristic frequency. ch ar-c Greater than threshold A c-th resh Meanwhile, the spectral entropy value S at that frequency c The spectral entropy value S under normal conditions c-normal The difference exceeds ΔS c-th resh This further confirms the possibility of a core grounding fault. c-th resh ΔS c-th resh The threshold is determined by comparing and analyzing a large amount of neutral point current spectrum data during core grounding faults and normal operation, and is used to judge the degree of abnormality of characteristic frequency amplitude and spectral entropy value.

[0331] (3) Insulation aging fault judgment

[0332] 1) Based on DBN model output probability: When the insulation aging failure probability P output by the DBN model... a Greater than threshold P a-th resh At that time, it was considered that the transformer might have an insulation aging fault. a-th resh It is determined based on historical data and experience related to transformer insulation aging. For example, through long-term monitoring and analysis of multiple transformers, it was found that when P... a When the value is greater than 0.5, the probability of the transformer showing signs of insulation aging is relatively high. P can be... a-thresh Set it to 0.5.

[0333] 2) Combining time characteristics: Analyze the power spectral density P of voltage V(t) and current I(t). V (f), P I(f). If within a specific frequency range (such as f3 - f4), the rate of change of the voltage power spectral density and the rate of change of the current power spectral density respectively exceed the thresholds and then it indicates that the frequency characteristics of the voltage and current have changed, which may be related to insulation aging.

[0334] f3 and f4 are boundary parameters for defining the frequency range when analyzing the rate of change of the voltage and current power spectral densities. f3 is the starting frequency of this frequency range, f4 is the ending frequency, and f3 < f4. The values of these two parameters are related to the electrical characteristics during the normal operation of the transformer and the frequency change characteristics under fault conditions. They are determined through spectral analysis of a large amount of transformer operation data and combined with theoretical research on frequency changes during transformer faults. In practical applications, faults such as transformer insulation aging will cause changes in the frequency components of the voltage and current signals, and these changes will be concentrated in certain specific frequency ranges. For example, in some research and actual monitoring, it has been found that when a transformer has an insulation aging fault, within the frequency range of 1000Hz - 5000Hz, the changes in the power spectral densities of the voltage and current are more obvious. At this time, f3 can be set to 1000Hz and f4 to 5000Hz. For different types of transformers and different operating environments, the values of f3 and f4 may vary, and need to be adjusted and optimized according to the characteristics and historical data of the specific transformer to ensure that frequency change information related to faults can be effectively captured within this frequency range and provide an accurate basis for fault determination.

[0335] The calculation formulas for the rate of change of the power spectral density are respectively:

[0336]

[0337] and are obtained through statistical analysis of the rate of change of the voltage and current power spectral densities during normal operation and insulation aging, and are used to measure the abnormality of the change in the power spectral density.

[0338] 3) Combining spatial characteristics: Analyze the probability density function p(Q) of partial discharge. If in some areas inside the transformer, the probability density p(Q) of the partial discharge quantity is greater than the threshold p Q-th resh , and the spatial autocorrelation function I(Q) of the partial discharge differs from that in the normal state by more than ΔI Q-th resh , then it indicates that the partial discharge activity is abnormal, which may be caused by a decrease in insulation performance due to insulation aging. p Q-th resh , ΔI Q-th reshIt is a threshold determined by statistical analysis of partial discharge data during normal operation and insulation aging, used to judge the degree of abnormality of partial discharge probability density and spatial autocorrelation function.

[0339] 4) Calculate the spectral entropy value S of the transformer partial discharge signal based on spectral characteristics. p If S p >S p-th resh This indicates increased spectral complexity and a more uneven distribution of frequency components, consistent with the characteristics of partial discharge spectrum changes during insulation aging. p-th resh It is obtained through statistical analysis of the spectral entropy values ​​of partial discharge signals during normal operation and insulation aging, and is used to define the normal range of spectral entropy values.

[0340] By combining the output results of the DBN model with the comprehensive judgment scheme of multi-source information features of time, space, and spectrum and corresponding thresholds, we can more accurately predict transformer winding short-circuit faults, core grounding faults, and insulation aging faults, providing strong support for transformer operation and maintenance, timely detection of potential fault hazards, and ensuring the safe and stable operation of the power system.

Claims

1. A transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion, characterized in that, The method includes the following steps: S1. Data acquisition, which includes: time information acquisition, spatial information acquisition and spectrum information acquisition; S2. Data preprocessing, which includes normalization and noise reduction operations; S3. Feature extraction, which includes: temporal feature extraction, spatial feature extraction and spectral feature extraction; S4. Model building and training, which includes: building a deep belief network (DBN) and using the whale optimization algorithm (WOA) to optimize the weight parameters of the DBN; S5. Fault prediction, the fault prediction includes: inputting test data into the trained DBN model, outputting the fault probability, and combining the threshold judgment rules of time, space and spectrum characteristics to comprehensively judge winding short circuit, core grounding and insulation aging faults. The time information acquisition process is as follows: Voltage signals during transformer operation are acquired in real time using a high-precision sensor. Current signal Power factor The voltage and current signals are dynamically modeled using a time-varying autoregressive moving average (TV-ARMA) model, and the power factor is calculated using an improved power factor formula. The spatial information acquisition process is as follows: sensors are arranged in the transformer windings, core, and tank walls to acquire spatial distribution data of temperature, vibration, and partial discharge parameters; and spatial autocorrelation function is used to obtain spatial distribution data of temperature, vibration, and partial discharge parameters. and variogram Analyze the spatial distribution characteristics; by using the probability density function, we can understand the spatial distribution probability of partial discharge and identify areas with high incidence of partial discharge. The process of acquiring the spectrum information is as follows: Spectral analysis is performed on the neutral point current, winding current, and partial discharge signal. Variational Mode Decomposition (VMD) combined with Adaptive Wavelet Packet Transform (AWPT) algorithm is used to extract spectral features. The extracted spectral features include: frequency band energy ratio. and center frequency ; The fault determination rules in the fault prediction include: Winding short circuit fault: DBN output probability and temperature gradient Frequency band energy ; Core grounding fault: DBN output probability And power factor skewness Spectral characteristic frequency amplitude ; Insulation aging fault: DBN output probability And the entropy value of the partial discharge spectrum Power spectral density change rate .

2. The transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion according to claim 1, characterized in that, The spatial autocorrelation function The calculation formula is: ； in, It is the average of all temperature data; spatial weight matrix , It is the sensor spacing. This is the weight decay parameter; The TV-ARMA model solves for the time-varying coefficients using the alternating direction multiplier method (ADMM). and It is updated in real time to reflect dynamic changes caused by load mutations.

3. The transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion according to claim 1, characterized in that, The normalization process is as follows: Z-score normalization is performed on the time, space, and spectral data respectively, and the data is mapped to the [-1,1] interval; The denoising process is as follows: a multi-scale morphological filtering algorithm is used to perform opening and closing operations on the data using structuring elements of different scales to remove noise.

4. The transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion according to claim 3, characterized in that, The multi-scale morphological filtering uses two-dimensional structuring elements to process spatial data, and the erosion and dilation calculation formulas are as follows: ； 。 5. The transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion according to claim 1, characterized in that, The process of extracting the time features is as follows: extracting the statistical features, autocorrelation coefficient, cross-correlation coefficient, and power spectral density of voltage, current, and power factor; the statistical features include: mean, variance, kurtosis, and skewness; The spatial feature extraction process is as follows: calculating the temperature gradient. Spatial differences in vibration intensity Partial discharge spatial positioning coordinates ; The process of extracting spectral features is as follows: extracting feature frequencies. Frequency band energy Center frequency ,bandwidth and spectral entropy .

6. The transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion according to claim 5, characterized in that, Characteristic frequency The search criteria are: ； in, and These are spectrum data Regarding frequency The first and second derivatives, It is a pre-set threshold; Bandwidth energy The calculation formula is: ； in, and They are frequency bands The corresponding discrete frequency point index, It is frequency resolution; Center frequency The calculation formula is: ； bandwidth The calculation formula is: ; Spectral entropy The calculation formula is: ； in, It is the number of discrete points in the spectrum data. Indicates the first The proportion of energy at each frequency point to the total energy.

7. The transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion according to claim 1, characterized in that, The process of constructing a deep belief network (DBN) is as follows: it consists of multiple stacked Restricted Boltzmann Machine (RBM) layers, with the number of input layer nodes matching the number of features, and the output layer corresponding to the fault type or probability; The process of optimizing DBN weight parameters using the Whale Optimization Algorithm (WOA) includes: initializing the population, fitness assessment, hunting behavior simulation, and iterative optimization until convergence.

8. The transformer fault prediction method based on spatiotemporal-spectral multi-source information fusion according to claim 7, characterized in that, The output layer of the DBN uses the softmax function for multi-class classification, with the following formula: ； in, It is the output layer. The output value of the nth node indicates that the sample belongs to the nth node. The probability of a type of failure; , It connects to the last RBM hidden layer unit. With output layer unit The weight, It is an output layer unit The bias, This represents the total number of fault types. The parameter update formula of the Whale Optimization Algorithm (WOA) includes: Encirclement behavior: ; Spiral behavior: ; Random search behavior: .

Images