A Fault Diagnosis Method for On-Load Tap Changer Based on Multi-Feature Fusion

Abstract

A fault diagnosis method for on-load tap changers based on multi-feature fusion, which relates to the technical field of on-load tap changer fault diagnosis and is used to solve the problem that the existing vibration signal analysis methods usually focus on the extraction of single features and cannot comprehensively capture the complex information in the signals. It includes data acquisition, feature extraction, and fault diagnosis. The present invention combines advanced feature extraction techniques of complementary ensemble empirical mode decomposition (CEEMD), quality factor adjustable wavelet transform (QFAT), and short-time Fourier transform (STFT) to comprehensively extract fault information in vibration signals. By optimizing the core parameters of the outlier support vector machine (OOSVM) classifier and combining it with the predator optimization algorithm to enhance the global search ability and local fine search ability of the algorithm, the accuracy, robustness, and timeliness of fault diagnosis are further improved, realizing the requirements of high precision, real-time performance, and adaptability for power equipment fault diagnosis.

Description

Technical Field

[0001] The present invention relates to the technical field of on-load tap-changer fault diagnosis, and specifically, it is a method for on-load tap-changer fault diagnosis based on multi-feature fusion. Background Art

[0002] The on-load tap-changer (OLTC) is an important component in a power transformer, undertaking the key task of regulating the output voltage. With the increase in the service life of power equipment and the influence of a complex operating environment, the OLTC is prone to faults such as mechanical wear and poor contact. If these faults are not detected in time, they will seriously affect the normal operation of the transformer and the stability of the power system. Traditional fault diagnosis methods rely on manual inspection and regular maintenance, unable to achieve real-time monitoring and with low diagnostic accuracy. Especially when facing complex fault patterns and noise interference, the adaptability and timeliness of traditional methods are poor, and it is difficult to meet the requirements of modern power equipment fault diagnosis.

[0003] Traditional feature extraction methods, such as time-domain analysis, frequency-domain analysis, and empirical mode decomposition (EMD), usually can only extract signal features for specific types of faults, lacking a comprehensive analysis of the multi-dimensional and dynamic changes of signals, and it is difficult to cope with diverse fault patterns in a complex electrical equipment operating environment. In addition, existing vibration signal analysis methods usually focus on the extraction of single features, unable to comprehensively capture the complex information in the signal, lacking the ability to comprehensively analyze vibration signals from multiple dimensions, such as time-frequency characteristics, frequency components, and non-linear characteristics, resulting in a low accuracy of fault diagnosis. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for on-load tap-changer fault diagnosis based on multi-feature fusion, which is used to solve the problem that existing vibration signal analysis methods usually focus on the extraction of single features, unable to comprehensively capture the complex information in the signal, and lacking the ability to comprehensively analyze vibration signals from multiple dimensions.

[0005] The technical solution adopted by the present invention to solve its technical problems is: a method for on-load tap-changer fault diagnosis based on multi-feature fusion, including the following steps.

[0006] S1 Data acquisition.

[0007] Collect the original vibration signal X i (t) of the on-load tap-changer under four working conditions: normal state, gear jamming, contact looseness, and fastening screw looseness.

[0008] S2 Feature extraction.

[0009] S2.1 Feature extraction of the on-load tap-changer vibration signal.

[0010] S2.1.1 Decompose the original vibration signal X i (t) into several intrinsic mode functions IMF by using complementary ensemble empirical mode decomposition CEEMD. Add a positive and a negative white noise m i (t) to the original vibration signal X i (t) respectively to obtain two different new signals P i and N i .

[0011] P i = X i (t)+ m i (t).

[0012] N i = X i (t)- m i (t).

[0013] Perform EMD decomposition on the above two different new signals respectively, obtain j IMF components and calculate their average value to get the final IMF component C j (t).

[0014] .

[0015] C ij + (t) represents the j-th intrinsic mode function IMF obtained after performing EMD decomposition on the signal P i ; C ij - (t) represents the j-th intrinsic mode function IMF obtained after performing EMD decomposition on the signal N i ; m is the number of white noises added during EMD decomposition.

[0016] Repeat the above steps, add a new normal distribution white noise sequence each time, and take the IMF with a larger correlation coefficient each time as the final result; find the 5th order component C i (t) with a larger correlation coefficient, and obtain the energy E i of each component, .

[0017] Construct the total energy E ∑ of the vibration signal according to the energy of each component: .

[0018] Find the energy entropy H EN of the IMF component: .

[0019] S2.1.2 Extract the subsequence energy E j of the on-load tap-changer vibration signal by using the tunable quality factor wavelet transform TQWT, The specific steps are as follows: First, determine the low-pass scale α and high-pass scale β of the filter to obtain a wavelet basis that conforms to the signal characteristics. The calculation formulas are as follows: ; where Q is the quality factor; r is the oversampling rate of the TQWT-transformed signal, which is calculated by dividing the sum of the wavelet coefficients after decomposition by the signal length; the low-pass output of each filter bank of the TQWT is used as the input of the continuous filter bank, and the maximum decomposition level Jmax of the multi-layer filter bank is: .

[0020] where N is the length of the analytic signal, is the floor function.

[0021] Then perform the decomposition of the TQWT; the frequency response function of the low-pass filter of each layer of the wavelet filter is set as H0(ω), and the frequency response function of the high-pass filter of each layer of the wavelet filter is set as H1(ω), w1, w2, …, w J+1 are the wavelet coefficients from the 1st to the (J + 1)th layer, and the calculation formula is as follows.

[0022] .

[0023] .

[0024] In the formula, θ(ω) is the functional function used to construct the frequency response functions H0(ω) and H1(ω), .

[0025] Obtain the spectrum X(ω) of the signal through Fourier transform and input it as the input signal into the filter bank of the TQWT; the input signal of the (j - 1)th layer is decomposed into two components, denoted as w (j) (ω) and v (j) (ω), w (j) (ω) is the high-quality factor component corresponding to the jth layer filter bank, and v (j) (ω) is the low-quality factor component corresponding to the jth layer filter bank; use v (j) (ω) as the input quantity of the next layer filter, perform the inverse Fourier transform on w (j) (ω), and the obtained time-domain component w (j) (t) is used as the wavelet decomposition coefficient of the jth layer of the TQWT; store the wavelet decomposition coefficients corresponding to the 1st to the Jth layers into the corresponding rows of the decomposition coefficient matrix W to obtain the J-layer decomposition coefficient matrix W, and the coefficient of the (J + 1)th layer is the decomposition residue.

[0026] After decomposition, the oscillation energy characteristics of each layer subsequence of the OLTC original vibration signal correspond to the filter of that layer and can characterize the vibration energy characteristics of the original vibration signal. For the energy {E1, E2, …, E J} Perform statistics to obtain the time-frequency domain characteristic quantity sequence of the original vibration signal; among them, the energy of the j-th layer subsequence is: .

[0027] S2.1.3 Adopt the vibration frequency component amplitude entropy VFCAE based on short-time Fourier transform to describe the irregularity of the vibration signal frequency components of the on-load tap-changer at different time periods.

[0028] First, convert the multi-channel signal into an STFT matrix, and then extract VFCAE to construct a feature vector.

[0029] For the discrete time series s = { s1, s2,.., s N}, its STFT is calculated as follows: .

[0030] Among them, h(n-k) is the Gaussian window function, and the value of S(f, k) is a complex number; s(n) is the discrete time series of the original signal, and e -j2πfn / N is the complex exponential part of the Fourier transform.

[0031] Calculate the Euclidean norm of the elements in S as follows: .

[0032] In the formula: ‖ ‖ is the Euclidean norm.

[0033] Based on the maximum value, divide the amplitude of the STFT matrix into L intervals, and count the distribution D(p) of the amplitude intervals, where p is an integer between 0 and L; then, add weights to the calculated entropy values. ; among them, ω(p) is the weight of the p-th interval; finally, calculate the average entropy value within a time period as a characteristic quantity, and take the expected feature vector dimension as y, then the characteristic quantity is: ; where k represents calculating the average entropy of the k-th feature dimension; ; represents rounding down.

[0034] S2.2 Feature fusion.

[0035] Perform standardization processing on the above three characteristic quantities F(k), E j and H EN .

[0036] .

[0037] Among them, μ F、 μ E、 μ H is the mean value of the characteristic quantity; σ F、 σ E、 σ HIt is the standard deviation of the characteristic quantity.

[0038] The eigenvectors are combined using the weighted summation method.

[0039] ; where ω1, ω2, ω3 are the set weights.

[0040] S3 Fault diagnosis.

[0041] The improved viper optimization algorithm-based optimization operator support vector machine IVOA-OOSVM is used to diagnose faults in the features.

[0042] The IOWA operator is defined as follows: .

[0043] where b j is the a i with the j-th largest u i in IOWA for (u i ), a i value, u i is the order-induced variable, a i is the independent variable; ω j represents the weight of the IOWA operator.

[0044] The expression of Shannon entropy is as follows: .

[0045] First, sort the slack variable ξ according to the vector u: .

[0046] In the formula, B u is the ordered parameter vector of the slack variable ξ, that is, b j is the slack variable corresponding to the IOWA with the j-th largest value in u; different weights are assigned according to the sorted positions, and the sorted slack variable ξ and the corresponding weight vector ω are defined; the weight vector ω is determined in advance according to the definition of the IOWA operator and satisfies and , where n is the number of data points; ω T is the transpose of the weight vector ω.

[0047] Calculate the weighted slack variable sum L(ξ), which is expressed by the mathematical formula as: .

[0048] where v(i) is the sorted index, indicating the position of the i-th largest slack variable ξ in the original vector.

[0049] The positive sample optimization equation of the optimization operator support vector machine is: .

[0050] The negative sample optimization equation for the optimized operator support vector machine is as follows: .

[0051] In the formula, matrix A represents the positive class samples, B represents the negative class samples; K is the kernel function to be determined; b1 and b2 are the biases; ξ and η are the slack variables; W1 and W2 are the normal vectors; e1 and e2 are the unit column vectors, the number of rows of e1 is the same as that of the kernel function and the number of rows of e2 is the same as that of the kernel function ; c1 and c2 are the penalty factors.

[0052] Use the improved viper optimization algorithm IVOA to optimize the core parameters of the support vector machine, namely the penalty factor C and the Gaussian kernel function parameter g.

[0053] Furthermore, by arranging sensors on the surface of the fuel tank, the original vibration signal X i (t) of the on-load tap-changer is collected.

[0054] Furthermore, the mathematical model of the improved viper optimization algorithm IVOA is as follows: The first stage: The prey recognition process.

[0055] The prey is freely distributed in the environment, and its behavior simulation expression is: .

[0056] In the formula, P i is the location of the target object of the i-th viper; X k is the state of the viper; k is a natural number in [1, N]; N is the viper population size.

[0057] The viper randomly selects a prey to attack, and the viper behavior simulation expression is as follows: .

[0058] In the formula, F ri is the ideal fitness value; is the new state of the i-th viper in the v-th dimension; r is an arbitrary value in [1, N], used to generate the irregular behavior of the viper; I takes the value of 1 or 2; x i，v is the original state of the i-th viper in the v-th dimension; p i，v is the position of the target of the i-th viper in the v-th dimension, and F i is the fitness value of the current position of the i-th viper.

[0059] .

[0060] In the formula, is the new state of the i-th viper; is the fitness value corresponding to the new state of the viper in the v-th dimension; calculate the current optimal solution : .

[0061] is the fitness value corresponding to each viper individual .

[0062] Initialize the population using backpropagation initialization to accelerate the search process: .

[0063] where is the new position of the i-th individual in the initialized population; is the adjustment factor; is the distance between the current solution and the optimal solution, representing the initialization direction of the population individuals: .

[0064] Introduce the wild horse algorithm search mechanism to improve the parameter r, and calculate the individual angular coordinate θ1(i): .

[0065] x i is the current iterative individual; x best is the optimal individual in the population; o is the origin; the value range of i is (1, N).

[0066] Divide the iterative region into 4 parts, s is the region serial number, and calculate the number of individuals in each region; assign a probability P(s) to the region according to the number of individuals, and the probability calculation expression is as follows: .

[0067] In the formula, num(s) is the number of individuals in the s-th region.

[0068] Adopt the fitness weight to search the region and calculate the region fitness weight A(s): .

[0069] Select the search region s by the way of region selection, and select the angle θ k in the s-th region, and its generation expression is as follows: .

[0070] where α is the adjustment factor that controls the intensity of the adjustment amount, α ∈ [0, 1]; d i is the Euclidean distance between the current individual and the optimal solution: ; d max is the possible maximum distance in the entire search space: , that is, the distance between the farthest individual and the optimal solution.

[0071] After generating θk Substitute into the following formula for iteration: .

[0072] In the formula, is the new state of the i-th viper in the v-th dimension after iteration.

[0073] Second stage: Chase and escape process.

[0074] After the viper catches the prey, the prey tries to escape; during the chase process, assuming that this hunting is close to an attack range with a radius of R, the mathematical expression of this process is as follows: , .

[0075] In the formula, is the new state of the i-th viper in the v-th dimension in the second stage; t is the current iteration number; T is the maximum iteration number.

[0076] .

[0077] In the formula, is the new state of the i-th viper in the second stage; is the fitness value in the new state.

[0078] Introduce the particle swarm optimization algorithm PSO to update the position of the viper, and the update expression is as follows: .

[0079] In the formula, is the new state of the i-th viper in the v-th dimension after introducing the PSO algorithm in the second stage; c1, c2 are learning factors that control the guiding degree of the individual and the group; r1, r2 are random numbers in the range of [0,1]; X best,v is the historical optimal position of the individual, X global,v is the global optimal position of the individual.

[0080] The beneficial effects of the present invention are as follows: The present invention combines advanced feature extraction technologies of complementary ensemble empirical mode decomposition CEEMD, quality factor adjustable wavelet transform QFAT, and short-time Fourier transform STFT to comprehensively extract fault information in vibration signals. By optimizing the core parameters of the OOSVM classifier and combining the predator optimization algorithm to improve the global search ability and local fine search ability of the algorithm, the accuracy, robustness, and timeliness of fault diagnosis are further improved, and the requirements of high precision, real-time performance, and adaptability for power equipment fault diagnosis are realized. Brief Description of the Drawings

[0081] Figure 1 is the flow chart of the diagnostic method of the present invention.

[0082] Figure 2 Figure showing gear jamming caused by the looseness of the transmission shaft

[0083] Figure 3 Figure showing the looseness of the contact

[0084] Figure 4 Figure showing the looseness of the fastening bolt

[0085] Figure 5 Flow chart of the optimization operator support vector machine based on the improved viper optimization algorithm Detailed implementation method

[0086] The following describes in detail a on-load tap-changer fault diagnosis method based on multi-feature fusion of the present invention with reference to the accompanying drawings

[0087] As Figure 1 shown, a on-load tap-changer fault diagnosis method based on multi-feature fusion includes the following steps

[0088] S1 Data acquisition

[0089] Collect the vibration signal of the on-load tap-changer by arranging sensors on the surface of the oil tank, simulate and collect the vibration signals of the on-load tap-changer under four working conditions of normal state, gear jamming, contact looseness, and fastening screw looseness, and conduct multiple groups of tests under the normal and mechanical fault states of the on-load tap-changer to collect multiple groups of data. As Figure 2 shown, it is a figure of gear jamming; as Figure 3 shown, it is a figure of contact looseness; as Figure 4 shown, it is a figure of fastening bolt looseness

[0090] S2 Feature extraction

[0091] S2.1 Feature extraction of the vibration signal of the on-load tap-changer

[0092] S2.1.1 Use complementary ensemble empirical mode decomposition (CEEMD) to decompose the original vibration signal into several intrinsic mode functions (IMFs). Add a positive and a negative white noise m i (t) to the original vibration signal X i (t) respectively to obtain two different new signals P i and N i .

[0093] P i = X i (t)+ m i (t).

[0094] N i = X i (t)- m i (t).

[0095] Perform EMD decomposition on these two different new signals respectively, obtain j IMF components and calculate their average values to get the final IMF component C j (t).

[0096] .

[0097] C ij + (t) represents the jth intrinsic mode component IMF obtained after performing EMD decomposition on the signal P i ; C ij - (t) represents the jth intrinsic mode component IMF obtained after performing EMD decomposition on the signal N i ; m is the number of white noises added during EMD decomposition.

[0098] Repeat the above steps, add a new white noise sequence with a normal distribution each time, and take the IMF with a relatively large correlation coefficient obtained each time as the final result. Find the 5th order component C i (t) with a relatively large correlation coefficient, and obtain the energy E of each component i : .

[0099] Construct the total energy E of the vibration signal according to the energy of each component ∑ : .

[0100] Find the energy entropy H of the IMF component EN : .

[0101] S2.1.2 Extract the subsequence energy E of the on-load tap-changer vibration signal by using the tunable quality factor wavelet transform TQWT j . The quality factor Q, redundancy parameter r, and decomposition level J are the main parameters of TQWT. Among them, the quality factor Q is a technical index for evaluating the energy concentration degree of the signal, which is the ratio of the central frequency corresponding to the maximum value in the signal frequency domain to the 3 dB bandwidth.

[0102] The specific steps are as follows: First, determine the low-pass scale α and high-pass scale β of the filter to obtain a wavelet basis that conforms to the signal characteristics. The calculation formulas are respectively:[[]] .

[0103] Among them, r is the oversampling rate of the TQWT-transformed signal, which is calculated by dividing the sum of the wavelet coefficients after decomposition by the signal length. The low-pass output of each filter bank of TQWT is used as the input of the continuous filter bank, and the maximum decomposition level Jmax of the multi-layer filter bank is:[[]] .

[0104] Where N is the length of the analytical signal, is the floor function.

[0105] Next, the decomposition of TQWT is performed. The frequency response function of the low-pass filter of each layer of the wavelet filter is set as H0(ω), and the frequency response function of the high-pass filter of each layer of the wavelet filter is set as H1(ω), w1, w2, …, w J+1 are the wavelet coefficients from the 1st to the J + 1st layer, and the calculation formula is as follows.

[0106] .

[0107] .

[0108] In the formula, θ(ω) is the functional function used to construct the frequency response functions H0 (1) (ω) and H1(ω), . ω is the angular frequency.

[0109] The spectrum X(ω) of the signal is obtained through Fourier transform and input into the filter bank of TQWT as the input signal. The input signal of the (j - 1)th layer is decomposed into two components, denoted as w (j) (ω) and v (j) (ω), which are the high- and low-quality factor components corresponding to the jth layer filter bank respectively. Take v (j) (ω) as the input quantity of the next layer filter, and perform the inverse Fourier transform on w (j) (ω). The obtained time-domain component w (j) (t) is used as the wavelet decomposition coefficient of the jth layer of TQWT. The wavelet decomposition coefficients corresponding to the 1st to the Jth layers are stored in the corresponding rows of the decomposition coefficient matrix W to obtain the J-layer decomposition coefficient matrix W, and the coefficient of the (J + 1)th layer is the decomposition residue.

[0110] After decomposition, the oscillation energy characteristics of each layer of subsequences of the OLTC vibration signal correspond to the filter of that layer and can characterize the vibration energy characteristics of the original signal. Therefore, the energies {E1, E2, …, E J} of the J-layer subsequences can be statistically analyzed as the time-frequency domain characteristic quantity sequence of the original vibration signal. Among them, the energy of the jth layer subsequence is: .

[0111] S2.1.3 Combining the advantages of time-frequency domain analysis and information entropy technology in strong anti-noise feature extraction, the vibration frequency component amplitude entropy VFCAE based on short-time Fourier transform is used to describe the irregularity of the vibration signal frequency components of the on-load tap-changer at different time periods. Among various signal processing technologies, the short-time Fourier transform STFT is applicable to long vibration signals to display the signal frequency components changing with time. First, the multi-channel signals are transformed into an STFT matrix, and then the VFCAE is extracted to construct the feature vector.

[0112] For the discrete time series s = { s1, s2,.., s N}, its STFT is calculated as follows: .

[0113] where h(n - k) is the Gaussian window function, and the value of S(f, k) is a complex number. s(n) is the discrete time series of the original signal, and e -j2πfn / N is the complex exponential part of the Fourier transform.

[0114] To convert it into a real matrix form, the Euclidean norm of the elements in S is calculated as follows: .

[0115] In the formula: ‖ ‖ is the Euclidean norm.

[0116] First, based on the maximum value, the amplitude of the STFT matrix is divided into L intervals, and the distribution D(p) of the amplitude intervals is counted, where p is an integer between 0 and L. Then, in order to highlight the influence of the frequency components with high amplitudes, weights are added to the calculated entropy values. The higher the amplitude, the greater the weight.

[0117] .

[0118] where ω(p) is the weight of the p-th interval. Finally, in order to reduce the size of the entropy feature, the average entropy value within a time period is calculated as the feature quantity. Assuming the expected feature vector dimension is y, then the feature quantity is: .

[0119] where k represents the index, indicating that the average entropy of the k-th feature dimension is being calculated; ; represents rounding down.

[0120] S2.2 Feature fusion.

[0121] Since the dimensions of the feature quantities are different, the above three feature quantities F(k), E j and H EN are normalized: .

[0122] Among them, μ and σ are the mean and standard deviation of the characteristic quantity respectively.

[0123] The feature vectors are combined using the weighted summation method: .

[0124] Among them, ω1, ω2, ω3 are the set weights.

[0125] S3 Fault diagnosis.

[0126] Such as Figure 5 shown, the improved viper optimization algorithm-based optimization operator support vector machine IVOA - OOSVM is used to diagnose the faults of the features.

[0127] The optimization operator support vector machine is improved on the basis of the twin vector machine. Its basic idea is to find two non-parallel hyperplanes by solving two quadratic programming problems, and use the IOWA operator to replace the hinge loss function generalization in order to make different trade-offs for the slack variable ξ. The IOWA operator is an average aggregation operator, which generalizes the OWA operator by using a reordering process based on induced variables. The n-dimensional IOWA operator is a mapping IOWA: R n →R, which has an n-dimensional weight vector ω associated with and . The IOWA operator is defined as follows: .

[0128] Among them, b j is the value of a i with the j-th largest u i in IOWA for (u i , a i ). u i is the order-induced variable, and a i is the independent variable. ω j represents the weight of the IOWA operator.

[0129] The expression of Shannon entropy is as follows: .

[0130] First, sort the slack variable ξ. According to the previous definition, sort the slack variable ζ according to the vector u: .

[0131] In the formula, B u is the ordered parameter vector of the slack variable ξ, that is, b jIt is the slack variable corresponding to the j-th largest IOWA among the u medians. Different weights are assigned according to the sorted positions, and these weights will be used to calculate the sum of weighted slack variables, thereby replacing ξ in the original hinge loss function. Define the sorted slack variable ξ and the corresponding weight vector ω. ω T is the transpose of the weight vector ω. The weight vector ω is predetermined according to the definition of the IOWA operator and satisfies and , where n is the number of data points.

[0132] Calculate the sum of weighted slack variables, which is expressed by the mathematical formula: .

[0133] where v(i) is the sorted index, representing the position of the i-th largest slack variable ξ in the original vector.

[0134] From this, the positive sample optimization equation of the optimized operator support vector machine can be obtained as follows: .

[0135] The negative sample optimization form of the optimized operator support vector machine is similar: .

[0136] In the formula, matrix A represents positive class samples, B represents negative class samples; K is the kernel function to be determined; b1 and b2 are biases; ξ and η are slack variables; W1 and W2 are normal vectors; e1 and e2 are unit column vectors, and the number of rows of e1 is the same as that of the kernel function , and the number of rows of e2 is the same as that of the kernel function ; c1 and c2 are penalty factors.

[0137] From the above analysis, it can be seen that the penalty factors c1, c2 and the Gaussian kernel function parameter g are two very important hyperparameters, which have a significant impact on the performance, generalization ability and computational complexity of the model. Use the improved viper optimization algorithm IVOA to optimize the core parameters of the support vector machine, namely the penalty factor C and the Gaussian kernel function parameter g.

[0138] The algorithm establishes a mathematical model based on different hunting stages of the assumed viper as follows: The first stage: the prey recognition process.

[0139] Prey are freely distributed in the environment, and their behavioral simulation expression is: .

[0140] In the formula, P i is the location of the target object of the i-th viper; X kFor the state of the viper; k is a natural number in [1, N]; N is the number of the viper population. The viper randomly selects a prey to attack, and the viper behavior simulation expression is as follows: .

[0141] In the formula, F ri is the ideal fitness value; is the new state of the i-th viper; is the new state of the i-th viper in the v-th dimension; r is an arbitrary value in [1, N], used to generate the irregular behavior of the viper; I takes the value of 1 or 2. x i，v is the original state of the i-th viper in the v-th dimension, p i，v is the position of the target of the i-th viper in the v-th dimension, that is, the position of the prey; F i is the fitness value of the current position of the i-th viper.

[0142] .

[0143] In the formula, is the fitness value corresponding to the new state of the viper in the v-th dimension. Calculate the current optimal solution : .

[0144] is the fitness value corresponding to each viper individual . In each iteration of the algorithm, the population distribution is dynamically adjusted according to the current optimal solution using backpropagation, and the updated optimal solution each time is used as a reference for the next population initialization. The position of the population is adjusted through error propagation, so that the population approaches the current optimal solution. If , then update the current optimal solution and its corresponding fitness.

[0145] Use backpropagation initialization to initialize the population to accelerate the search process: .

[0146] Among them, is the new position of the i-th individual in the initialized population, is the adjustment factor, is the distance between the current solution and the optimal solution, representing the initialization direction of the population individuals: .

[0147] Introduce the search mechanism of the wild horse algorithm to improve the parameter r. Calculate the individual angular coordinate θ1(i): .

[0148] x i is the current iteration individual; x bestis the optimal individual of the population; o is the origin; the value range of i is (1, N).

[0149] Divide the iterative region into 4 parts, s is the region serial number, and calculate the number of individuals in each region. Assign the probability P(s) to the region according to the number of individuals. The more individuals, the greater the probability. The probability calculation expression is as follows: .

[0150] In the formula, num(s) is the number of individuals in the s-th region. Use the fitness weight to search the region and calculate the region fitness weight A(s): .

[0151] Select the search region s by the way of region selection, and select the angle θ in the s-th region k , and its generation expression is as follows: .

[0152] Among them, α is the adjustment factor that controls the intensity of adjustment, α ∈ [0,1]; d i is the Euclidean distance between the current individual and the optimal solution: ; d max is the possible maximum distance in the entire search space: , that is, the distance between the farthest individual and the optimal solution.

[0153] After generating θ k Substitute it into the following formula for iteration: .

[0154] In the formula, is the new state of the i-th viper in the v-th dimension after iteration.

[0155] The second stage: the chasing and escaping process.

[0156] After the viper catches the prey, the prey tries to escape. During the chasing process, assume that this hunting is close to an attack range with a radius of R. The mathematical expression of this process is as follows: .

[0157] .

[0158] In the formula, is the new state of the i-th viper in the v-th dimension in the second stage; t is the current iteration number; T is the maximum iteration number.

[0159] .

[0160] In the formula, is the new state of the i-th viper in the second stage; is the fitness value in the new state.

[0161] In the second stage of the algorithm, the prey tries to escape. Since the iteration parameter R is defined within a limited range, the particle swarm optimization (PSO) algorithm is introduced to update the viper's position to avoid the algorithm falling into a local optimum. The update expression is as follows: .

[0162] In the formula, is the new state of the i-th viper in the v-th dimension after introducing the PSO algorithm in the second stage; c1 and c2 are learning factors that control the guiding degree of the individual and the group; r1 and r2 are random numbers within the range [0, 1]; X best,v and X global,v are the historical optimal position of the individual and the global optimal position respectively.

[0163] The model is trained using known fault samples. The extracted fault features are input into the trained fault diagnosis model to obtain the fault diagnosis result. The optimized operator support vector machine based on the improved viper optimization algorithm has better global search ability and can conduct global exploration in a complex environment. By mimicking the search behavior of predators in nature, the algorithm can find the global optimal solution in a broader solution space. Compared with other optimization algorithms, this algorithm has better iterative convergence ability and faster optimization speed for both multimodal functions and unimodal functions, demonstrating the superiority of the optimized operator support vector machine based on the improved viper optimization algorithm.

Claims

1. A fault diagnosis method for on-load tap changer based on multi-feature fusion, characterized in that: The following steps are involved: S1 Data Collection Collect the original vibration signals of the on-load tap changer under four working conditions: normal state, gear jamming, contact loosening, and fastening screw loosening. i (t); S2 feature extraction S2.1 On-load tapchanger vibration signal feature extraction S2.1.1 The complementary ensemble empirical mode decomposition (CEEMD) is used to transform the original vibration signal X i (t) is decomposed into several intrinsic mode components IMF, and the energy entropy H of the intrinsic mode components IMF is calculated. EN ; S2.1.2 Using quality factor adjustable wavelet transform TQWT to extract the subsequence energy E of the on-load tap changer vibration signal j ; S2.1.3 uses the short-time Fourier transform vibration frequency component amplitude entropy VFCAE to describe the irregularity of the frequency component of the on-load tapchanger vibration signal in different time periods; calculate the average entropy value F(k) within a time period; S2.2 Feature Fusion For the above three feature quantities F(k), E j and H EN To standardize: Among them, μ F , μ E , μ H is the mean of the characteristic quantity; σ F , σ E , σ H is the standard deviation of the characteristic quantity; The eigenvectors are combined using the weighted summation method: F fused =ω1·F norm +ω2·E j，norm +ω3·H EN，norm ; Among them, ω1, ω2, ω3 are the set weights; S3 Troubleshooting The optimization operator support vector machine IVOA-OOSVM based on the improved Viper optimization algorithm is used to perform fault diagnosis on the features; The specific content of step S2.1.3 is: first, convert the multi-channel signal into an STFT matrix, and then extract the VFCAE to construct a feature vector; For a discrete time series s = {s1, s2, .., s N }, its STFT is calculated as follows: Among them, h(nk) is the Gaussian window function, the value of S(f,k) is a complex number; s(n) is the discrete time series of the original signal, e -j2πfn / N is the complex exponential part of the Fourier transform; The Euclidean norm of the elements in S is calculated as follows: Where: ‖‖ is the Euclidean norm; The amplitude of the STFT matrix is ​​divided into L intervals based on the maximum value, and the distribution of the amplitude intervals is counted D(p), where p is an integer between 0 and L; then, a weight is added to the calculated entropy value; Among them, ω(p) is the weight of the pth interval; finally, the average entropy value in a time period is calculated as the feature value, and the expected feature vector dimension is taken as y, then the feature value is: Among them, k means that the average entropy of the kth feature dimension is being calculated; Indicates rounding down; The mathematical model of the improved Viper optimization algorithm IVOA is: Stage 1: Prey Identification Process The prey is freely distributed in the environment, and its behavior simulation expression is: : Where P i is the location of the target object of the i-th viper; X k is the state of the viper; k is a natural number [1, N]; N is the number of viper populations; The viper randomly selects a prey to attack. The viper behavior simulation expression is as follows: In the formula, F ri is the ideal fitness value; is the new state of the ith viper in the vth dimension; r is an arbitrary value in [1,N], which is used to generate irregular behavior of the viper; the value of I is 1 or 2; x i，v is the original state of the ith viper in the vth dimension; p i，v is the position of the target of the ith viper in the vth dimension, F i is the fitness value of the current position of the i-th viper; In the formula, is the new state of the i-th viper; is the fitness value corresponding to the new state of the viper in the vth dimension; calculate the current optimal solution X best : F(X i ) for each viper individual X i The corresponding fitness value; Initialize the population using backpropagation initialization to speed up the search process: in, is the new position of the i-th individual in the initial population; α is the adjustment factor; ΔX i is the distance between the current solution and the optimal solution, indicating the initialization direction of the population individuals: The wild horse algorithm search mechanism is introduced to improve the parameter r and calculate the individual angle coordinate θ1(i): x i is the current iteration individual; x best is the optimal individual in the population; o is the origin; the value range of i is (1, N); The iteration area is divided into 4 parts, s is the area number, and the number of individuals in each area is calculated; the probability P(s) is assigned to the area according to the number of individuals, and the probability calculation expression is as follows: Where num(s) is the number of individuals in the sth region; Use fitness weight to search the area and calculate the regional fitness weight A(s): Select the search area s by area selection, and select the angle θ in the sth area k , and its generating expression is as follows: Among them, α is the adjustment factor that controls the intensity of the adjustment, α∈[0,1]; d i is the Euclidean distance between the current individual and the optimal solution: =X-Xpesel; d max is the maximum possible distance in the entire search space: d max =max||X i -X bese ||, that is, the distance between the farthest individual and the optimal solution; After the generation k Substitute into the following formula and iterate: In the formula, is the new state of the i-th viper in the v-th dimension after iteration; Stage 2: Chase and Escape After the viper catches its prey, the prey tries to escape; during the chase, assuming that the hunting is close to an attack range with a radius of R, the mathematical expression of this process is as follows: In the formula, is the new state of the ith viper in the vth dimension in the second stage; t is the current number of iterations; T is the maximum number of iterations; In the formula, is the new state of the i-th viper in the second stage; is the fitness value in the new state; The particle swarm optimization algorithm PSO is introduced to update the position of the viper. The update expression is as follows: In the formula, is the new state of the i-th viper in the v-th dimension after the introduction of the PSO algorithm in the second stage; c1, c2 are learning factors that control the degree of guidance of individuals and groups; r1, r2 are random numbers in the range of [0, 1]; X best,v is the individual's historical optimal position, X global,v is the global optimal position of an individual.

2. The on-load tap changer fault diagnosis method based on multi-feature fusion according to claim 1, characterized in that: The original vibration signal X of the on-load tap changer is collected by placing sensors on the surface of the oil tank. i (t).

3. The on-load tap changer fault diagnosis method based on multi-feature fusion according to claim 2 is characterized in that: The specific steps of step S2.1.1 are: To the original vibration signal X i (t) adds a positive and a negative white noise m respectively i (t), we get two different new signals P i and N i : P i =X i (t)+m i (t); N i =X i (t)-m i (t); Perform EMD decomposition on the above two different new signals respectively, obtain j IMF components and calculate their average values ​​to obtain the final IMF component C j (t); C ij + (t) represents the signal P i The jth intrinsic mode component IMF obtained after EMD decomposition; C ij - (t) represents the signal N i The jth intrinsic mode component IMF obtained after EMD decomposition; m is the amount of white noise added during EMD decomposition; Repeat the above steps, adding a new normally distributed white noise sequence each time, and taking the IMF with the largest correlation coefficient each time as the final result; calculate the fifth-order component C with the largest correlation coefficient i (t), and find the energy E of each component i :E i =∫|C i (t)| 2 dt,i=1,2,…,n; Construct the total energy E of the vibration signal based on the energy of each component ∑ : Find the IMF component energy entropy H EN :

4. The on-load tap changer fault diagnosis method based on multi-feature fusion according to claim 1, characterized in that: The specific steps of step S2.1.2 are as follows: First, determine the low-pass scale α and high-pass scale β of the filter to obtain the wavelet basis that meets the signal characteristics. The calculation formulas are: Where Q is the quality factor; r is the oversampling rate of the TQWT transformed signal, which is calculated by dividing the sum of the decomposed wavelet coefficients by the signal length; the low-pass output of each filter group of TQWT is used as the input of the continuous filter group, and the maximum number of decomposition layers Jmax of the multi-layer filter group is: Where N is the length of the analytical signal, To round down; Then, the TQWT is decomposed; the low-pass filter frequency response function of each layer of wavelet filter is set to H0(ω), and the high-pass filter frequency response function of each layer of wavelet filter is set to H1(ω), w1, w2, …, w J+1 is the wavelet coefficient from the 1st to the J+1th layer, and the calculation formula is: Where θ(ω) is the performance function, which is used to construct the frequency response functions H0(ω) and H1(ω). The spectrum X(ω) of the signal is obtained by Fourier transform and is input into the filter bank of TQWT as the input signal; the input signal of the (j-1)th layer is decomposed into two components, denoted as w (j) (ω) and v (j) (ω), w (j) (ω) is the high quality factor component corresponding to the j-th filter bank, v (j) (ω) is the low quality factor component corresponding to the j-th filter bank; v (j) (ω) is used as the input of the next layer of filters, and w (j) (ω) is transformed by inverse Fourier transform, and the time domain component w is obtained (j) (t) is used as the j-th wavelet decomposition coefficient of TQWT; the wavelet decomposition coefficients corresponding to layers 1 to J are stored in the corresponding rows of the decomposition coefficient matrix W to obtain the J-th layer decomposition coefficient matrix W, and the (J+1)-th layer coefficient is the decomposition residue; The oscillation energy characteristics of each subsequence of the decomposed OLTC original vibration signal correspond to the filter of this layer, which can characterize the vibration energy characteristics of the original vibration signal. J } to perform statistics as a time-frequency domain feature sequence of the original vibration signal; Among them, the energy of the j-th layer subsequence is:

5. The on-load tap changer fault diagnosis method based on multi-feature fusion according to claim 1, characterized in that: The specific content of step S3 is: The IOWA operator is defined as follows: where b j is IOWA for (u i , a i ) has the jth largest u i a i Value, u i is the order-inducing variable, a i is the independent variable; ω j Represents the weight of the IOWA operator; The expression of Shannon entropy is as follows: First, sort the slack variables ξ according to the vector u: In the formula, B u is the ordered parameter vector of the slack variable ξ, i.e. b j is the slack variable corresponding to the jth largest IOWA in u; different weights are assigned according to the sorted positions, and the sorted slack variables ξ and the corresponding weight vector ω are defined; the weight vector ω is predetermined according to the definition of the IOWA operator, satisfying ω i ≥0 and where n is the number of data points; ω T is the transpose of the weight vector ω; Calculate the weighted slack variables and L(ξ), expressed mathematically as: Among them, v(i) is the sorted index, indicating the position of the i-th largest slack variable ξ in the original vector; The positive sample optimization equation of the optimization operator support vector machine is: The negative sample optimization equation of the optimization operator support vector machine is: In the formula, matrix A represents positive samples, B represents negative samples; K is the kernel function to be determined; b1 and b2 are biases; ξ and η are slack variables; W1 and W2 are normal vectors; e1 and e2 are unit column vectors, and the number of rows in e1 is the same as the kernel function K(A, [A, B] T ), the number of rows in e2 is the same as that of the kernel function K(B, [A, B] T ) are the same; c1 and c2 are penalty factors; The improved Viper optimization algorithm IVOA is used to optimize the core parameter penalty factor C and Gaussian kernel function parameter g of the support vector machine.

Images