Transformer fault diagnosis method based on multi-sensor multi-domain feature fusion

Abstract

The invention discloses a transformer fault diagnosis method based on multi-sensor multi-domain feature fusion, and the method comprises the following steps: 1, obtaining a vibration signal data set of a transformer box, and carrying out the decentralization of a collected signal; 2, decomposing the decentralized transformer vibration signal by adopting a dung beetle optimized ICEEMD method, calculating a Pearson correlation coefficient between an IMF component and an original signal, and selecting a strongly correlated component to reconstruct the transformer vibration signal, so as to filter and denoise the transformer vibration signal; step 3, extracting time domain, frequency domain and space domain features in the transformer vibration signals by adopting an MVAR model, VMD-BIS decomposition and Riemannian manifold, and constructing a joint feature vector; 4, performing fusion dimension reduction by adopting a kernel principal component analysis algorithm to obtain a fusion dimension reduction feature vector; and a fifth step of constructing a kernel extreme learning machine model, training parameters of the kernel extreme learning machine model in combination with the time-frequency-space fusion features, and optimizing hyper-parameters of the kernel extreme learning machine model in combination with a Harris eagle optimization algorithm so as to obtain an optimal kernel extreme learning machine classification model. According to the method, the problems of misjudgment, missed judgment and low diagnosis rate caused by unbalanced distribution due to nonlinear characteristics of vibration signals in transformer fault diagnosis are effectively solved.

Description

Technical Field

[0001] The present invention relates to the research field of transformer fault diagnosis technology, and in particular to a transformer fault diagnosis method based on multi-sensor multi-domain feature fusion. Background Art

[0002] As one of the important components of my country's power system, the health status of the transformer directly affects the stable operation of the power system. In order to ensure the safe operation of the power system and prevent economic losses caused by power system shutdown, an efficient and accurate transformer fault diagnosis model must be established.

[0003] With the rise of neural networks and intelligent technologies, transformer fault diagnosis methods continue to emerge, such as BP neural networks, support vector machines, random forests, Bayesian methods, expert systems, etc., which have been widely used and have greatly improved accuracy. However, although the above methods are flexible and have achieved good application results, the model may misjudge or miss faults due to the nonlinear characteristics of the vibration signal at the fault location. In addition, many actual operating scenarios are outdoors. Due to the harsh working conditions of the transformer, the strong noise environment increases the difficulty of collecting vibration signals and extracting effective sample features, resulting in poor model performance and reduced accuracy. Summary of the invention

[0004] In view of the above problems, the present invention proposes a transformer fault diagnosis method with multi-sensor and multi-domain feature fusion. Based on time-frequency-space domain feature fusion, the multi-core extreme learning machine optimized by Harris Eagle is adopted to analyze its power system characteristics, so as to effectively solve the problems of fault misjudgment, missed judgment, poor model performance and low accuracy caused by the nonlinear characteristics of vibration signals.

[0005] The technical solution of the present invention is as follows: A transformer fault diagnosis method based on multi-domain feature fusion of sensors comprises the following steps:

[0006] The first step is to install sensors at key locations of the transformer to collect vibration signals from the transformer box and perform decentralized processing on the vibration signals.

[0007] In the second step, the ICEEMD algorithm of Dung Beetle Optimization (DBO) is used to decompose the transformer box vibration signal to obtain multiple modal components IMF; the Pearson correlation coefficient of each modal component is calculated, the strongly correlated components are screened and the signal is reconstructed to complete the transformer box vibration signal data denoising;

[0008] In the third step, the MVAR model is used to extract the time domain features of the reconstructed transformer vibration signal, the VMD-BIS algorithm is used to extract the frequency domain features of the reconstructed transformer vibration signal, and the Riemann manifold is used to extract the spatial domain features of the transformer vibration signal.

[0009] The fourth step is to concatenate the time-domain, frequency-domain and space-domain feature vectors head to tail to construct a time-frequency-space joint feature vector, extract the joint feature vectors of all samples of various fault types, and construct a joint feature matrix.

[0010] The fifth step is to perform kernel principal component analysis on the joint feature matrix, select the number of principal components according to the principle that the cumulative contribution rate is greater than 85%, and project the sample data into the new principal component space to obtain the fused feature vector after dimensionality reduction.

[0011] The sixth step is to construct a kernel extreme learning machine model, train its parameters by combining the time-frequency-space fusion features, and optimize its hyperparameters by combining the Harris Eagle optimization algorithm to obtain the optimal kernel extreme learning machine classification model.

[0012] Preferably, the specific steps of the first step are as follows:

[0013] Step 1: Install sensors at key locations of the transformer to collect transformer box vibration signal sample data matrix

[0014] X is the signal matrix constructed from the original sampling signal. Each row in the matrix is ​​the collected signal of a sensor, and each column is a sampling point, that is, m is the number of sensors and n is the number of sampling points.

[0015] Step 2: Specifically, the vibration signals of the transformer box in three fault states, namely, loose transformer winding, loose core, and winding deformation, and in a normal working state are collected.

[0016] Step 3: Use the common average reference to decentralize the transformer box vibration signal to obtain a decentralized matrix

[0017] Preferably, the specific steps of the second step are as follows:

[0018] Step 1: Initialize the parameters of the dung beetle population, including the position of the individual dung beetle population vector [white noise amplitude weight (Nstd) and the number of noise additions (NE)], the number of iterations, etc., and divide the dung beetle population into rolling dung beetles, egg balls, small dung beetles, and stealing dung beetles in a ratio of 6:6:7:11;

[0019] Step 2: The rolling ball dung beetle simulates the rolling ball behavior and moves along the given direction in the entire search space. During the rolling process, the position of the rolling ball dung beetle is updated. When the dung beetle encounters an obstacle and cannot move forward, it needs to reposition itself by dancing in order to obtain a new route.

[0020] Step 3: Simulate the boundary selection strategy of the dung beetle's egg-laying area. Once the egg-laying area is determined, the female dung beetle will lay eggs, and one egg will be laid in each iteration;

[0021] Step 4: The grown-up dung beetles come out to forage for food. Establish the best dung beetle area and guide the dung beetles to forage for food.

[0022] Step 5: Thieves will steal food from other dung beetles and locate the best food source;

[0023] Step 6: Repeat the above process, update the position information of the dung beetle, determine whether the target agent is beyond the boundary, update the current optimal solution and its fitness value, until the maximum number of iterations is reached, and obtain the optimal parameter combination of the ICEEMD algorithm [white noise amplitude weight (Nstd) and number of noise additions (NE)];

[0024] Step 7: Bring the optimal parameter combination of white noise amplitude weight (Nstd) and noise addition times (NE) into the ICEEMD algorithm to decompose the signal and obtain n modal component IMFs;

[0025] Step 8: Calculate the Pearson correlation coefficient of each modal component IMF;

[0026] Step 9: Filter out the components with strong correlation and reconstruct the signal to obtain the reconstructed signal Y i (t);

[0027] Preferably, the specific steps of the third step are as follows:

[0028] Step 1: The MVAR model is a generalized form of the univariate AR model. The MVAR model is used to extract the time domain characteristics of the transformer vibration signal. The p-order MVAR model can be described as follows:

[0029]

[0030]

[0031] Step 2: Construct the time domain feature vector V of the transformer vibration signal corresponding to m sensors AR , the expression is:

[0032]

[0033] Step 3: Y i (t) Perform variational mode decomposition and decompose it into K IMF components; variational mode decomposition is to obtain the optimal solution of the constrained variational model through iterative search, automatically obtain the center frequency and bandwidth of the intrinsic mode function, Y i The constrained variational model of (t) is as follows:

[0034]

[0035]

[0036] In the formula {u k}={u 1 ,u 2 ,…u k} and {ω k}={ω 1 ,ω 2 ,…ω k} are the IMF component set and the center frequency set corresponding to the IMF component set, and the alternating direction multiplier method is used to iteratively update u k ,ω k , solve the saddle point of the augmented Lagrangian function and continuously update u k ,ω k , the iteration terminates if the given judgment expression is satisfied.

[0037] Step 4: Screen the number of IMF components based on the center frequency, and determine the number of IMF components K according to the principle of no over-decomposition. Specifically, when decomposed into K components, if the difference between the center frequency of the K+1th component and the center frequency of the Kth component and the difference between the center frequency of the Kth component and the center frequency of the K-1th component are in the same order of magnitude, then over-decomposition occurs when decomposing into K+1 components, and the decomposition into K IMF components is selected;

[0038] Step 5: Perform bispectral analysis on the K IMF components obtained by variational mode decomposition. The analysis process is as follows:

[0039] Let {u k}={u 1 ,u 2 ,…u k} is the IMF component decomposed by variational mode. According to the discrete Fourier transform of the signal, the bispectrum is defined as:

[0040]

[0041] Step 6: Where ω j is the frequency variable, U i (ω i ) is {u k The discrete Fourier transform of each IMF component in} is used, E[] is the mathematical expectation, and the indirect method is used for bispectrum estimation. First, the third-order cumulant of the obtained IF component is estimated, and then the cumulant is subjected to a two-dimensional discrete Fourier transform, thus obtaining the bispectral estimation of the IMF component.

[0042] Step 7: Calculate all B in the main domain of the bispectral plane x (ω1 ,ω 2 ) is taken as the frequency domain feature, and the mathematical expression formula of the bispectral feature is as follows:

[0043]

[0044] Step 8: Use variational mode decomposition-bispectral analysis to extract the frequency domain features of the selected m channels, and perform signal smoothing and normalization processing on the extracted frequency domain features. The normalization processing formula is as follows:

[0045]

[0046] Step 9: Where B is the bispectral feature after signal smoothing and normalization. is the bispectral feature vector extracted from the selected m channels, M represents The data length, min and max are operators for calculating the maximum and minimum values; the constructed frequency domain feature vector V VMD-BIS The form is as follows:

[0047] V VMD-BIS =[B 1.1 ,B 1.2 ,...,B 1.k ]

[0048] Step 10: Use Riemann manifold to extract spatial features. This study uses the covariance matrix of the transformer vibration signal as the statistical feature of the data. The spatial covariance matrix of the EEG signal of a single participant is:

[0049]

[0050] Step 11: Since the covariance matrix is ​​symmetric and the off-diagonal items are calculated twice during the calculation process, only the upper triangular part of the matrix is ​​taken as the spatial domain feature, so its off-diagonal elements need to be scaled to maintain the standard. Vectorize it into a spatial domain feature vector:

[0051]

[0052] Preferably, the specific steps of the fourth step are as follows:

[0053] Step 1: Construct a joint feature vector. Connect the time-frequency domain feature vectors in sequence to construct the following joint feature vector:

[0054]

[0055] The joint feature vector f is the union of the features in each domain extracted from the same sample data using different signal processing methods, and is a t-dimensional feature vector.

[0056] Step 2: Construct a combined feature matrix, and construct an n*t combined feature matrix T from the selected n sample data

[0057] T = [f 1 , f 2 , …, f n T

[0058] Each row of the matrix represents the combined feature vector of a sample.

[0059] Preferably, the specific steps of the fifth step are as follows:

[0060] Step 1: Project the combined feature matrix T into a linearly separable high-dimensional space through the kernel function φ for kernel principal component analysis to obtain the principal component coefficient matrix K:

[0061]

[0062] During principal component analysis, the sample covariance matrix can be obtained. It consists of the eigenvalues forming a column vector, and the eigenvalues are arranged in descending order, and the eigenvalues λ 1 ≥ λ 2 ≥ … ≥ λ t and the corresponding eigenvectors α κ (k = 1, …, t);

[0063] Step 2: The ratio reflects the contribution rate of the k-th principal component in the generated feature space. According to the feature

[0064] value variance cumulative contribution rate formula, determine the number S of kernel principal components. The calculation formula is as follows:

[0065]

[0066] Step 3: Select the cumulative contribution rate of the first s (s < t) principal components to exceed 85%. Therefore, the principal component coefficient matrix becomes:

[0067]

[0068] Project the combined feature vector f constructed by each sample onto the new principal components, and the fused feature vector f ∈ R after dimensionality reduction is obtained. The projection method is as follows: 1×s , and the projection method is as follows:

[0069]

[0070] Preferably, the specific steps of the sixth step are as follows:

[0071] ​Step 1: Build a kernel extreme machine model, select a mixed weighted combination of radial basis kernel function and polynomial kernel function as the kernel function of the kernel extreme learning machine, and achieve dimensionality reduction by multiplying the original high-dimensional data input by the weight matrix between the hidden layer and the output layer:

[0072]

[0073] Among them, K(y i ,y j ) is the hybrid kernel function of the kernel extreme learning machine, η, ξ are the parameters of the mixing kernel function;

[0074] Step 2: Use the Harris Eagle optimization algorithm to optimize the parameters in step 1. First, construct the objective function of parameter optimization through the validation set recognition accuracy. The objective function is defined as follows:

[0075]

[0076] Among them, ACC represents the overall recognition accuracy of the validation set, NC i is the number of correct classifications of the i-th category, NE i is the number of misclassifications of the i-th category, and l is the total number of categories;

[0077] Secondly, the Harris Eagle optimization algorithm is applied to the above objective function to obtain the optimal values ​​of the parameters of the hybrid kernel extreme learning machine through the global search phase, the search to development phase, and the development phase;

[0078] The global search phase updates the value of the parameter vector (Harris Hawk) through two strategies:

[0079]

[0080] Among them, H(t) is the current position vector of the eagle, H rabbit (t) is the current position of the individual with the best fitness, H rand (t) is the position of the eagle randomly selected in the current population, H(t+1) is the position vector of the eagle in the next iteration, r 1 , r 2 , r 3 , r 4 and q are random numbers between (0,1), LB and UB are the lower and upper limits of the parameters to be optimized, respectively. m (t) is the average position of the population;

[0081] The transition between the global search phase and the development phase is realized according to the value of |E|. When |E| < 1, the development phase is entered. The expression of E is:

[0082]

[0083] Where E represents the escape energy of the prey, E 0 represents the initial state of energy, T represents the maximum number of iterations;

[0084] In the development phase, four attack methods are selected according to the energy E of the prey and the escape probability gl of the prey to update the position of the parameter vector, where:

[0085] When |E|≥0.5 and gl≥0.5, a soft encirclement attack is used;

[0086] When |E|≤0.5 and gl≥0.5, a hard encirclement attack is used;

[0087] When |E|≥0.5 and gl<0.5, a soft encirclement attack method of progressive rapid dive is adopted;

[0088] When |E|≤0.5 and gl<0.5, a hard encirclement attack method with progressive rapid dive is adopted.

[0089] The present invention is based on the fusion of time-frequency-space domain features and adopts the multi-core extreme learning machine optimized by Harris Eagle to analyze the characteristics of its power system and effectively solve the problems of fault misjudgment, missed judgment, poor model performance and low accuracy caused by the nonlinear characteristics of vibration signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0090] Figure 1 is a flowchart of the method of the present invention;

[0091] Figure 2 Flowchart of the optimization algorithm for Harris Hawk. DETAILED DESCRIPTION

[0092] In order to make the purpose, technical solution and advantages of the present invention clearer, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. The specific embodiments described here are only used to explain the present invention and are not used to limit the present invention.

[0093] The first step is to collect the vibration signal of the transformer box and decentralize the vibration signal. The specific steps are as follows:

[0094] S101: Install sensors at key locations of the transformer to collect transformer box vibration signal sample data matrix

[0095] Specifically, four sensors are installed at key positions of the transformer to collect 25 groups of transformer box vibration signal data in three fault states, namely loose winding, loose core, and deformed winding, and in normal working state, with a data length of 2048 data points.

[0096] S102: Decentralize the transformer box vibration signal using a common average reference. The specific formula is as follows:

[0097]

[0098] Where Y i (t) is the original transformer vibration signal of the i-th sensor.

[0099] In the second step, the Dung Beetle Optimization (DBO) algorithm is used to optimize the optimal ICEEMD parameters of the vibration signal under each state, and the vibration signal is decomposed by ICEEMD to obtain multiple modal components IMF; the Pearson correlation coefficient of each modal component is calculated, and the components with strong correlation and above are selected and the signal is reconstructed to complete the denoising of the transformer fault signal data;

[0100] S201: Initializing the population parameters of dung beetles, including the position of the individual population vector of dung beetles [white noise amplitude weight (Nstd) and the number of noise additions (NE)], the number of iterations, etc., and dividing the dung beetle population into rolling ball dung beetles, egg balls, small dung beetles, and stealing dung beetles in a ratio of 6:6:7:11;

[0101] S202: The rolling dung beetle simulates rolling behavior according to the light intensity and moves along a given direction in the entire search space. During the rolling process, the position of the rolling dung beetle is updated, and the formula is as follows:

[0102] Y i (t+1)=Y i (t)+α*k*Y i (t-1)+b*ΔY

[0103] ΔY=|Y i (t)-Y w |

[0104] Where t is the number of iterations, Yi(t) is the position information of the i-th dung beetle at the t-th iteration, k∈(0, 0.2] is a constant, denoted as the deflection coefficient, b∈(0, 1) is a fixed value, α is 1 or -1, -1 means deviation from the original direction, 1 means no deviation from the original direction, Y w Recorded as the global worst position, ΔY simulates the light intensity change, where k is set to 0.1 and b is set to 0.3;

[0105] When a dung beetle encounters an obstacle and cannot move forward, it needs to reposition itself by dancing in order to obtain a new route. The formula is as follows:

[0106] Y i (t+1)=Y i (t)+tanθ|Y i (t)-Y i (t-1)|

[0107] In the formula, θ∈[0,π] is recorded as the deflection angle of the direction. When θ=0, When θ = π, the dung beetle does not update its position;

[0108] S203: Simulate the boundary selection strategy of the area where the female dung beetle lays eggs. Once the egg-laying area is determined, the female dung beetle will lay eggs, and lay eggs once in each iteration. The formula is as follows:

[0109] R=1-t / Tmax

[0110] Lb * =max(X * ×(1-R), Lb)

[0111] Ub * =min(X * ×(1+R),Ub)

[0112] B i (t+1)=X * +b 1 ×(B i (t)-Lb * )+b 2 *(B i (t)-Ub * )

[0113] Where, X * is recorded as the global optimal position, Lb is recorded as the lower bound of the optimization problem, Ub is recorded as the upper bound of the optimization problem, and Lb * Recorded as the lower bound of the spawning area, Ub * The upper limit of the spawning area is B i (t) is recorded as the position information of the i-th egg ball at the t-th iteration, b 1 With b 2 Denoted as 1×D independent random vectors, D is the dimension of the optimization problem;

[0114] S204: The grown-up dung beetles come out to forage for food. An optimal dung beetle area is established to guide the dung beetles to forage for food. The formula is as follows:

[0115] Lb b =max(C b ×(1-R), Lb)

[0116] Ub b =min(X 0 ×(1+R),Ub)

[0117] Y i (t+1)=Y i (t)+C 1 ×(Yi (t)-Lb b )+C 2 *(Y i (t)-Ub b )

[0118] Where Y i (t) is recorded as the position information of the i-th dung beetle at the t-th iteration, Lb b Denote as the lower bound of the optimal feeding area, Ub b Recorded as the upper bound of the optimal feeding area;

[0119] S205: Stealing dung beetles will steal food from other dung beetles. The formula for determining the location of the optimal food source is as follows:

[0120] Y i (t+1)=X b +s×g×(|Y i (t)-X * |+|Y i (t)-X b )

[0121] Yi(t) is denoted as the position information of the i-th stealing dung beetle at the t-th iteration, g is denoted as a random vector of size 1×D, serving the normal distribution, and S is a constant value.

[0122] S206: Repeat the above process, update the position information of the dung beetle, determine whether the target agent exceeds the boundary, update the current optimal solution and its fitness value, until the maximum number of iterations is reached, and obtain the best parameter combination of the ICEEMD algorithm [white noise amplitude weight (Nstd) and noise addition times (NE)];

[0123] S207: Bring the optimal parameter combination [white noise amplitude weight (Nstd) and noise addition times (NE)] into the ICEEMD algorithm to decompose the signal and obtain 10 to 12 modal component IMFs;

[0124] S208: Calculate the Pearson correlation coefficient of each modal component IMF;

[0125] S209: Filter out the components with strong correlation or above to reconstruct the signal, and obtain the reconstructed signal Yi(t);

[0126] In the third step, the AR model is used to extract the time domain features of the reconstructed transformer vibration signal, and the VMD-BIS algorithm is used to extract the frequency domain features of the reconstructed transformer vibration signal.

[0127] S301: The MVAR model is a generalized form of the univariate AR model. The MVAR model is used to extract the time domain features of the transformer vibration signal. The p-order MVAR model can be described as follows:

[0128]

[0129]

[0130] S302: Constructing the time domain feature vector V of the transformer vibration signal corresponding to m sensors AR , the expression is:

[0131]

[0132] S303: For y i (t) Perform variational mode decomposition and decompose it into K IMF components; variational mode decomposition is to obtain the optimal solution of the constrained variational model through iterative search, and automatically obtain the center frequency and bandwidth of the intrinsic mode function, y i The constrained variational model of (t) is as follows:

[0133]

[0134]

[0135] In the formula {u k}={u 1 ,u 2 ,…u k} and {ω k}={ω 1 ,ω 2 ,…ω k} are the IMF component set and the center frequency set corresponding to the IMF component set, and the alternating direction multiplier method is used to iteratively update u k ,ω k , solve the saddle point of the augmented Lagrangian function and continuously update u k ,ω k , the iteration terminates if the given judgment expression is satisfied.

[0136] S304: Screening the number of IMF components based on the center frequency, and determining the number of IMF components K according to the principle of not over-decomposing, specifically, when decomposing into K components, when the difference between the center frequency of the K+1th component and the center frequency of the Kth component and the difference between the center frequency of the Kth component and the center frequency of the K-1th component are in the same order of magnitude, over-decomposition occurs in decomposing into K+1 components, and decomposition into K IMF components is selected;

[0137] S305: Perform bispectral analysis on the K IMF components obtained by variational mode decomposition. The analysis process is as follows:

[0138] Let {u k}={u 1 ,u 2 ,…u k} is the IMF component decomposed by variational mode. According to the discrete Fourier transform of the signal, the bispectrum is defined as:

[0139]

[0140] S306: where ω j is the frequency variable, U i (ω i ) is {u k The discrete Fourier transform of each IMF component in} is used, E[] is the mathematical expectation, and the indirect method is used for bispectrum estimation. First, the third-order cumulant of the obtained IF component is estimated, and then the cumulant is subjected to a two-dimensional discrete Fourier transform, thus obtaining the bispectral estimation of the IMF component.

[0141] S307: Calculate all B in the main domain of the bispectral plane x (ω 1 ,ω 2 ) is taken as the frequency domain feature, and the mathematical expression formula of the bispectral feature is as follows:

[0142]

[0143] S308: Using variational mode decomposition-bispectral analysis method to extract frequency domain features of the selected m channels, and performing signal smoothing and normalization processing on the extracted frequency domain features. The normalization processing formula is as follows:

[0144]

[0145] S309: Wherein B is the bispectral feature after signal smoothing and normalization, is the bispectral feature vector extracted from the selected m channels, M represents The data length, min and max are operators for calculating the maximum and minimum values; the constructed frequency domain feature vector V VMD-BIS The form is as follows:

[0146] V VMD-BIS =[B 1.1 ,B 1.2 ,...,B 1.k ]

[0147] S310: Riemann manifold is used to extract spatial features. This study uses the covariance matrix of the transformer vibration signal as the statistical feature of the data. The spatial covariance matrix of the EEG signal of a single participant is:

[0148]

[0149] S311: Since the covariance matrix is ​​symmetric and the off-diagonal items are calculated twice during the calculation process, only the upper triangular part of the matrix is ​​taken as the spatial domain feature, so its off-diagonal elements need to be scaled to maintain the standard. Vectorize it into a spatial domain feature vector:

[0150]

[0151] The fourth step is to concatenate the time domain and frequency domain feature vectors to construct a time-frequency joint feature vector, extract the time-frequency joint feature vectors of all samples of various fault types, and construct a time-frequency joint feature matrix.

[0152] S401: Construct a joint feature vector, connect the time-frequency domain feature vectors in sequence, and construct the following joint feature vector:

[0153]

[0154] The joint feature vector f is the union of the features in each domain extracted from the same sample data using different signal processing methods, and is a t-dimensional feature vector.

[0155] S402: Construct a joint feature matrix, construct an n*t joint feature matrix T from the selected n sample data

[0156] T=[f 1 ,f 2 ,…,f n ] T

[0157] Each row of the matrix represents the joint feature vector of a sample.

[0158] The fifth step is to perform kernel principal component analysis on the time-frequency joint feature matrix, select the number of principal components according to the principle that the cumulative contribution rate is greater than 85%, and project the sample data into the new principal component space to obtain the fused feature vector after dimensionality reduction.

[0159] S501: Project the joint feature matrix T into a linearly separable high-dimensional space through the kernel function φ to perform kernel principal component analysis to obtain the principal component coefficient matrix K:

[0160]

[0161] The sample covariance matrix can be obtained while performing principal component analysis. It consists of eigenvectors formed by eigenvalues, and the eigenvalues are arranged in descending order, obtaining eigenvalues λ 1 ≥λ 2 ≥…≥λ t and the corresponding eigenvectors α κ (k = 1,…,t);

[0162] S502: The ratio reflects the contribution rate of the k-th principal component in the generated feature space. According to the

[0163] cumulative contribution rate formula of eigenvalue variance, determine the number of kernel principal components S. The calculation formula is as follows:

[0164]

[0165] S503: Select the first s (s < t) principal components whose cumulative contribution rate exceeds 85%. Therefore, the principal component coefficient matrix becomes:

[0166]

[0167] Project the combined feature vector f constructed by each sample onto the new principal components to obtain the fused feature vector after dimensionality reduction The projection method is as follows:

[0168]

[0169] Sixth step, construct a kernel extreme learning machine model, train its parameters by combining time-frequency fused features, and optimize its hyperparameters by combining the Harris hawk optimization algorithm to obtain the optimal kernel extreme learning machine classification model:

[0170] S601: Construct a kernel extreme machine model, select a mixed weighted radial basis kernel function and polynomial kernel function as the kernel function of the kernel extreme learning machine, and realize dimensionality reduction by multiplying the original high-dimensional data input by the weight matrix between the hidden layer and the output layer:

[0171]

[0172] Among them, K(y i ,y j ) is the mixed kernel function of the kernel extreme learning machine, η, ξ are the parameters of the mixed kernel function;

[0173] S602: Use the Harris hawk optimization algorithm to optimize the parameters in S601. First, construct the objective function for parameter optimization through the recognition accuracy of the validation set. The objective function is defined as follows:

[0174]

[0175] Among them, ACC represents the overall recognition accuracy of the validation set, NC i is the number of correct classifications of the i-th category, NE i is the number of misclassifications of the i-th category, and l is the total number of categories;

[0176] Secondly, the Harris Eagle optimization algorithm is applied to the above objective function to obtain the optimal values ​​of the parameters of the hybrid kernel extreme learning machine through the global search phase, the search to development phase, and the development phase;

[0177] The global search phase updates the value of the parameter vector (Harris Hawk) through two strategies:

[0178]

[0179] Among them, H(t) is the current position vector of the eagle, H rabbit (t) is the current position of the individual with the best fitness, H rand (t) is the position of the eagle randomly selected in the current population, H(t+1) is the position vector of the eagle in the next iteration, r 1 , r 2 , r 3 , r 4 and q are random numbers between (0,1), LB and UB are the lower and upper limits of the parameters to be optimized, respectively. m (t) is the average position of the population;

[0180] The transition between the global search phase and the development phase is realized according to the value of |E|. When |E| < 1, the development phase is entered. The expression of E is:

[0181]

[0182] Where E represents the escape energy of the prey, E 0 represents the initial state of energy, T represents the maximum number of iterations;

[0183] In the development phase, four attack methods are selected according to the energy E of the prey and the escape probability gl of the prey to update the position of the parameter vector, where:

[0184] When |E|≥0.5 and gl≥0.5, a soft encirclement attack is used;

[0185] When |E|≤0.5 and gl≥0.5, a hard encirclement attack is used;

[0186] When |E|≥0.5 and gl<0.5, a soft encirclement attack method of progressive rapid dive is adopted;

[0187] When |E|≤0.5 and gl<0.5, a hard encirclement attack method with progressive rapid dive is adopted.

[0188] Example:

[0189] The present invention first installs 4 sensors at key positions of the sample transformer to collect 25 groups of transformer box vibration signal data in three fault states of transformer winding looseness, core looseness, winding deformation and normal working state, with a data length of 2048 data points.

[0190] The present invention uses a dung beetle optimization algorithm to optimize the best ICEEMD parameters of vibration signals under various states, the dung beetle optimization algorithm initialization parameter is set to 30, wherein the dung beetle population is divided into rolling ball dung beetles, egg balls, small dung beetles, and stealing dung beetles in a ratio of 6:6:7:11, wherein the optimization problem dimension D is set to 2, the maximum number of iterations is 20, the white noise amplitude weight (Nstd) is 0.2 and the number of noise additions (NE) is 50, which are brought into ICEEMD for modal decomposition, and finally 10 to 12 IMF components are obtained, the Pearson correlation coefficients of the 10 to 12 IMF components are calculated respectively, the IMF components with correlation coefficients ≥0.6 are screened out, and finally 5 IMF components that meet the requirements are screened to complete signal reconstruction and perform data denoising;

[0191] The 4th-order MVAR model is used to extract the MVAR model coefficients of the transformer vibration signal and construct the time domain feature vector V AR , extract the frequency domain features of the transformer vibration signal to obtain the 5th order frequency domain feature vector V VMD-BIS , triangular Riemann space characteristic V mr , concatenate the dimensional feature vector, take 94 representative samples, construct a 100*94 feature matrix M, and then input it into the Harris Eagle optimized multi-core extreme learning machine for classification.

[0192] Based on the background of harsh working conditions of transformers, the patent of this invention proposes a transformer fault diagnosis method with multi-sensor multi-domain feature fusion, and combines time-frequency-space feature fusion with Harris Eagle optimized multi-core extreme learning machine to perform transformer fault diagnosis. Compared with the existing technology, the advantage of this patent is that it introduces the optimization algorithm based on multi-domain feature fusion and combines it with the machine learning framework to complete fault diagnosis, thereby achieving the effect of solving the problem of low diagnostic accuracy of transformer fault diagnosis and the problem of misjudgment and missed judgment of faults by the transformer fault diagnosis model caused by the nonlinear characteristics of vibration signals.

[0193] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, it should be understood by those skilled in the art that the technical solutions described in the aforementioned embodiments can still be modified, or some or all of the technical features therein can be replaced by equivalents. Therefore, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope defined by the claims of the present invention.

Claims

1. A transformer fault diagnosis method based on multi-sensor multi-domain feature fusion, characterized in that: The following steps are involved: The first step is to install multiple sensors at key locations of the transformer to collect vibration signals from the transformer box and perform decentralized processing on the vibration signals. In the second step, the ICEEMD algorithm of Dung Beetle Optimization (DBO) is used to decompose the transformer box vibration signal to obtain multiple modal components IMF; the Pearson correlation coefficient of each modal component is calculated, the strongly correlated components are screened and the signal is reconstructed to complete the transformer box vibration signal data denoising; In the third step, the MVAR model is used to extract the time domain features of the reconstructed transformer vibration signal, the VMD-BIS algorithm is used to extract the frequency domain features of the reconstructed transformer vibration signal, and the Riemann manifold is used to extract the spatial domain features of the transformer vibration signal; The fourth step is to concatenate the time-domain, frequency-domain and space-domain feature vectors head to tail to construct a time-frequency-space joint feature vector, extract the joint feature vectors of all samples of various faults, and construct a joint feature matrix; The fifth step is to perform kernel principal component analysis on the joint feature matrix, select the number of principal components according to the principle that the cumulative contribution rate is greater than 85%, and project the sample data into the new principal component space to obtain the fusion feature vector after dimensionality reduction; The sixth step is to construct a kernel extreme learning machine model, train its parameters by combining the time-frequency-space fusion features, and optimize its hyperparameters by combining the Harris Eagle optimization algorithm to obtain the optimal kernel extreme learning machine classification model.

2. The transformer fault diagnosis method based on multi-domain feature fusion of sensors according to claim 1 is characterized in that: The specific steps of the first step are as follows: Step 1: Install multiple sensors at key locations of the transformer to collect transformer box vibration signal sample data matrix X is the signal matrix constructed from the original sampling signal. Each row in the matrix is ​​the collected signal of a sensor, and each column is a sampling point, that is, m is the number of sensors, and n is the number of sampling points; Step 2: Specifically, collecting transformer box vibration signals in three fault states of transformer winding looseness, iron core looseness, winding deformation and normal working state; Step 3: Use the common average reference to decentralize the transformer box vibration signal to obtain a decentralized matrix 3. The transformer fault diagnosis method based on multi-domain feature fusion of sensors according to claim 1 is characterized in that: The specific steps of the second step are as follows: Step 1: Initialize the parameters of the dung beetle population, including the position of the individual dung beetle population vector [white noise amplitude weight (Nstd) and the number of noise additions (NE)], the number of iterations, etc., and divide the dung beetle population into rolling dung beetles, egg balls, small dung beetles, and stealing dung beetles in a ratio of 6:6:7:11; Step 2: The rolling ball dung beetle simulates the rolling ball behavior and moves along the given direction in the entire search space. During the rolling process, the position of the rolling ball dung beetle is updated. When the dung beetle encounters an obstacle and cannot move forward, it needs to reposition itself by dancing in order to obtain a new route. Step 3: Simulation The boundary selection strategy of the dung beetle's egg-laying area. Once the egg-laying area is determined, the female dung beetle will lay eggs, one egg per iteration. Step 4: The grown-up dung beetles come out to forage for food. Establish the best dung beetle area and guide the dung beetles to forage for food. Step 5: Thieves will steal food from other dung beetles and locate the best food source; Step 6: Repeat the above process, update the position information of the dung beetle, determine whether the target agent is beyond the boundary, update the current optimal solution and its fitness value, until the maximum number of iterations is reached, and obtain the optimal parameter combination of the ICEEMD algorithm [white noise amplitude weight (Nstd) and number of noise additions (NE)]; Step 7: Bring the optimal parameter combination of white noise amplitude weight (Nstd) and noise addition times (NE) into the ICEEMD algorithm to decompose the signal and obtain n modal component IMFs; Step 8: Calculate the Pearson correlation coefficients of each modal component IMF; Step 9: Filter out the components with strong correlation and reconstruct the signal to obtain the reconstructed signal Y i (t).

4. The transformer fault diagnosis method based on multi-domain feature fusion of sensors according to claim 1 is characterized in that: The specific steps of the third step are as follows: Step 1: The MVAR model is a generalized form of the univariate AR model. The time-domain characteristics of the transformer vibration signal are extracted using the MVAR model. The MVAR model of order p can be described as follows: Step 2: Construct the time domain feature vector V of the transformer vibration signal corresponding to m sensors AR , the expression is: Step 3: Y i (t) Perform variational mode decomposition and decompose it into K IMF components; variational mode decomposition is to obtain the optimal solution of the constrained variational model through iterative search, automatically obtain the center frequency and bandwidth of the intrinsic mode function, Y i The constrained variational model of (t) is as follows: In the formula {u k }={u1,u2,…u k } and {ω k }={ω1,ω2,…ω k } are the IMF component set and the center frequency set corresponding to the IMF component set, and the alternating direction multiplier method is used to iteratively update u k ,ω k , solve the saddle point of the augmented Lagrangian function and continuously update u k ,ω k , if the given judgment expression is satisfied, the iteration is terminated; Step 4: Screen the number of IMF components based on the central frequency. According to the principle of no over-decomposition, determine the number of IMF components K. Specifically, when decomposed into K components, the difference between the central frequency of the (K + 1)th component and the central frequency of the Kth component, and the difference between the central frequency of the Kth component and the central frequency of the (K - 1)th component are in the same order of magnitude, then decomposing into K + 1 components results in over-decomposition, and choose to decompose into K IMF components; Step 5: Perform bispectrum analysis on the K IMF components decomposed by variational mode decomposition. The analysis process is as follows: Let {u k }={u1,u2,…u k } is the IMF component decomposed by variational mode. According to the discrete Fourier transform of the signal, the bispectrum is defined as: Step 6: Where ω j is the frequency variable, U i (ω i ) is {u k The discrete Fourier transform of each IMF component in}, E[] is the mathematical expectation, and the indirect method is used for bispectral estimation. First, the third-order cumulant of the obtained IF component is estimated, and then the cumulant is subjected to a two-dimensional discrete Fourier transform, thus obtaining the bispectral estimation of the IMF component; Step 7: Calculate all B in the main domain of the bispectral plane x The sum of the absolute logarithms of (ω1, ω2) is used as the frequency domain feature, and the mathematical expression formula of the bispectral feature is as follows: Step 8: Use the variational mode decomposition-bispectrum analysis method to extract the frequency-domain characteristics of the selected m channels, and perform signal smoothing and normalization processing on the extracted frequency-domain characteristics. The normalization processing formula is as follows: Step 9: Where B is the bispectral feature after signal smoothing and normalization. is the bispectral feature vector extracted from the selected m channels, M represents The data length, min and max are operators for calculating the maximum and minimum values; the constructed frequency domain feature vector V VMD-BIS The form is as follows: V VMD-BIS =[B 1.1 ,B 1.2 ,...,B 1.k ] Step 10: Use Riemannian manifold to extract spatial domain characteristics. In this study, the covariance matrix of the transformer vibration signal is used as the statistical feature of the data. The spatial covariance matrix of the EEG signal of a single participant is: Step 11: Since the covariance matrix is symmetric and the off-diagonal terms are calculated twice during the calculation, only the upper triangular part of the matrix is taken as the spatial domain feature. Therefore, its off-diagonal elements need to be scaled to maintain the standard; vectorize it into a spatial domain feature vector:

5. The transformer fault diagnosis method based on multi-domain feature fusion of sensors according to claim 1 is characterized in that: The specific steps of the fourth step are as follows: Step 1: Construct a joint feature vector by connecting the time-frequency domain feature vectors in sequence at the head and tail to construct the following joint feature vector: The joint feature vector f is the combination of the features in each domain extracted from the same sample data using different signal processing methods, and is a t-dimensional feature vector; Step 2: Construct a joint feature matrix by constructing an n*t joint feature matrix T from the selected n sample data T=[f1,f2,…,f n ] T Each row of the matrix represents the joint feature vector of a sample.

6. The transformer fault diagnosis method based on multi-domain feature fusion of sensors according to claim 1 is characterized in that: The specific steps of the fifth step are as follows: Step 1: Project the joint feature matrix T into a linearly separable high-dimensional space through the kernel function φ to perform kernel principal component analysis to obtain the principal component coefficient matrix K: The sample covariance matrix can be obtained at the same time as the principal component analysis. It consists of eigenvalues ​​that form column vectors, and the eigenvalues ​​are arranged in descending order. The eigenvalues ​​λ1≥λ2≥…≥λ t and the corresponding eigenvector α κ (k=1,…,t); Step 2: Ratio It reflects the contribution rate of the kth principal element in the generated feature space. Determine the number of kernel principal components S according to the value variance cumulative contribution rate formula. The calculation formula is as follows: Step 3: Select the cumulative contribution rate of the first s (s < t) principal components to exceed 85%. Therefore, the principal component coefficient matrix becomes: Project the joint feature vector f constructed by each sample onto the new principal element to obtain the fusion feature vector after dimensionality reduction The projection method is as follows:

7. The transformer fault diagnosis method based on multi-domain feature fusion of sensors according to claim 1 is characterized in that: The specific steps of the sixth step are as follows: Step 1, construct a kernel extreme learning machine model. Select the mixed weighting of the radial basis kernel function and the polynomial kernel function as the kernel function of the kernel extreme learning machine, and realize dimensionality reduction by multiplying the input of the original high-dimensional data by the weight matrix between the hidden layer and the output layer: Among them, K(y i ,y j ) is the hybrid kernel function of the kernel extreme learning machine, η, ξ are the parameters of the hybrid kernel function; Step 2, use the Harris hawk optimization algorithm to optimize the parameters in Step 1. First, construct the objective function for parameter optimization through the recognition accuracy of the validation set. The objective function is defined as follows: Among them, ACC represents the overall recognition accuracy of the validation set, NC i is the number of correct classifications of the i-th category, NE i is the number of misclassifications of the i-th category, and l is the total number of categories; Secondly, the Harris Eagle optimization algorithm is applied to the above objective function to obtain the optimal values ​​of the parameters of the hybrid kernel extreme learning machine through the global search phase, the search to development phase, and the development phase; The global search phase updates the value of the parameter vector (Harris Hawk) through two strategies: Among them, H(t) is the current position vector of the eagle, H rabbit (t) is the current position of the individual with the best fitness, H rand (t) is the position of the eagle randomly selected in the current population, H(t+1) is the position vector of the eagle in the next iteration, r1, r2, r3, r4 and q are all random numbers between (0,1), LB and UB are the lower and upper limits of the parameters to be optimized, respectively, and H m (t) is the average position of the population; The transition between the global search phase and the development phase is realized according to the value of |E|. When |E| < 1, the development phase is entered. The expression of E is: Among them, E represents the escape energy of the prey, E0 represents the initial state of energy, and T represents the maximum number of iterations; In the development phase, four attack methods are selected according to the energy E of the prey and the escape probability gl of the prey to update the position of the parameter vector, where: When |E|≥0.5 and gl≥0.5, a soft encirclement attack is used; When |E|≤0.5 and gl≥0.5, a hard encirclement attack is used; When |E|≥0.5 and gl<0.5, a soft encirclement attack method of progressive rapid dive is adopted; When |E|≤0.5 and gl<0.5, a hard encirclement attack method with progressive rapid dive is adopted.

Images