Breaker vibration signal composite feature space construction method based on multi-channel fusion

Abstract

The invention discloses a circuit breaker vibration signal composite feature space construction method based on multi-channel fusion, which comprises the following steps of: 1) acquiring a multi-channel circuit breaker operating mechanism vibration signal, and constructing a noise-free reference signal and a standard signal-to-noise ratio test signal; 2) performing CFA-MVMD decomposition on the multi-channel circuit breaker vibration signal to obtain a fusion mode component and a frequency spectrum characteristic thereof; 3) performing segmentation operation on the obtained fusion mode component, performing segmentation time-frequency feature extraction, and constructing a segmentation time-frequency feature matrix; 4) extracting global time-frequency features, and constructing a global time-frequency feature matrix; combining segmentation and global time-frequency feature extraction to construct a multi-channel fused circuit breaker vibration signal composite feature space; and 5) inputting the multi-channel vibration signal composite characteristics of the circuit breaker into the CFA-KELM fault diagnosis model to obtain an early fault state of the circuit breaker. According to the invention, a high-precision and low-cost solution is provided for real-time monitoring of early mechanical faults of the circuit breaker.

Description

Technical Field

[0001] The present invention relates to a method for constructing a composite feature space of breaker vibration signals based on multi-channel fusion, and belongs to the technical field of breaker vibration signal processing. Background Art

[0002] With the development of the power system in recent years, China is forming a power supply mode of distributed power sources and centralized power grids. The demand for breakers to disconnect large currents and large capacities is gradually increasing. However, with the increase in the number of switching operations, the performance of breakers will gradually deteriorate, and faults may occur during long-term operation. The losses caused by breaker faults far exceed their own costs and even threaten the safety of the entire power supply line. According to the worldwide statistics of high-voltage breaker faults at the 2014 International Conference on Large High Voltage Electric Systems, the faults of breaker operating mechanisms and their auxiliary circuits account for 59.9% of all faults and show an increasing trend year by year. In the power system, preventive maintenance can reduce 30% of the maintenance cost, avoid 35%-45% of accidents, and reduce 75% of the downtime.

[0003] The vibration signal waveform of the circuit breaker operating mechanism is sensitive to the changes in its mechanical structure, with sufficient characteristic expression and is convenient for non-intrusive measurement. In recent years, a large number of scholars have carried out research on the state diagnosis of circuit breakers based on vibration signals combined with feature extraction and fault classification algorithms. The main research method is to collect multi-channel vibration signals and then perform necessary preprocessing. Then, various signal processing methods are used to extract the characteristic parameters contained in the preprocessed signals. Finally, the mechanical faults of the operating mechanism are diagnosed by combining data fusion diagnosis technology and state recognition technology. Common feature extraction methods include time-frequency analysis methods such as wavelet transform, modal decomposition methods such as EMD and VMD, and other methods such as phase space reconstruction and singular value decomposition. (1) Detection method of closing and opening coil current signals: The changes in the closing and opening coil current waveforms can intuitively reflect the changes in the states of the electromagnet and the operating mechanism during the operation of the circuit breaker. The closing and opening current waveforms of the circuit breaker contain the working state information of the circuit breaker during this process. By monitoring the closing and opening currents, various fault types in the circuit breaker control loop can be judged, such as coil core jamming, too low power supply voltage, too long core idle stroke, etc. At the same time, the easy acquisition of the coil current makes it very suitable for fault diagnosis of circuit breakers. (2) Vibration signal detection method: During the closing and opening processes of the circuit breaker, the energy changes rapidly, and the operating mechanism bears a huge impact force instantaneously. The position changes of the spring energy storage, rocker arm, cam, pawl, and pull rod, as well as the movement and stillness of the contacts, will all generate significant vibration signals. These signals are carriers with very rich information, containing a large amount of equipment state information, which is sufficient to reflect the key mechanical characteristic parameters such as the closing and opening speeds of the circuit breaker. The vibration signals have small attenuation and strong anti-interference ability, and can fully reflect the closing and opening states of the circuit breaker. (3) Contact stroke signal detection method: The contact stroke signal contains motion characteristic parameters such as closing and opening times, initial closing and opening speeds, stroke, and overtravel. Obtain the detailed motion trajectory and stroke-time curve of the transmission components during the operation of the circuit breaker. Displacement, main shaft rotation process, and angle parameters of the transmission components can be obtained from the curve. Since the motion between the transmission components and the moving contact is rigid, the closing and opening times of the moving contact, the stroke of the moving contact, and the motion speed can be directly calculated through the parameters on the curve, accurately depicting the motion trajectory of the moving contact of the circuit breaker during closing and opening.

[0004] With the development of informatization and intelligent technologies, artificial intelligence algorithms provide a feasible solution for the early mechanical fault state diagnosis of circuit breakers. The early mechanical fault diagnosis based on machine learning and deep learning has gradually attracted the attention of domestic and foreign researchers. The process of such fault diagnosis methods is usually to denoise and decompose signals such as vibration, sound, and current, further extract features, analyze the features through a classification model, and finally classify the fault data and output the results, and put forward opinions on maintenance. However, there are still problems such as insufficient feature extraction and low diagnostic accuracy in the application of current artificial intelligence technologies in the field of circuit breaker fault diagnosis. The main defects existing in the current fault diagnosis of circuit breaker operating mechanisms are:

[0005] (1) For minor mechanical defects such as loose screws and fatigue springs, traditional models rely on a single feature space such as time-domain statistical features or frequency-domain energy spectra for classification, which are easily affected by equipment individual differences and environmental fluctuations, and the identification ability of traditional models is insufficient in scenarios with weak feature differences.

[0006] (2) For the vibration monitoring of circuit breakers, multiple sensors are usually deployed. Traditional methods mostly use single-channel independent analysis or simple linear weighted fusion, ignoring the phase coupling and modal correlation of signals in different channels in the time-frequency domain, and failing to effectively solve the frequency-domain alignment and non-linear fusion problems of multi-source signals, resulting in the lack of fault feature dimensions and incomplete input information for the classification model. Summary of the Invention

[0007] The purpose of the present invention is to provide a method for constructing a composite feature space of circuit breaker vibration signals based on multi-channel fusion. First, the CFA-MVMD multi-channel fusion decomposition method is used to introduce the correlation between different channels to obtain the decomposed modal components of multi-channel information fusion. Then, the segmented time-frequency feature extraction is performed on the decomposed modal components. Finally, the global time-frequency features are combined to construct a composite feature space that can fully reflect the early faults of the circuit breaker.

[0008] The purpose of the present invention is realized through the following technical solutions:

[0009] A method for constructing a composite feature space of circuit breaker vibration signals based on multi-channel fusion, comprising the following steps:

[0010] 1) Collect the vibration signals of the multi-channel circuit breaker operating mechanism, and construct a noise-free reference signal and a standard signal-to-noise ratio test signal;

[0011] 2) Perform CFA-MVMD decomposition on the multi-channel circuit breaker vibration signals to obtain the fusion modal components and their spectral characteristics;

[0012] 3) Perform a segmentation operation on the obtained fusion modal components, perform segmented time-frequency feature extraction, and construct a segmented time-frequency feature matrix;

[0013] 4) Extract the global time-frequency features and construct a global time-frequency feature matrix; combine the segmented and global time-frequency feature extractions to construct a composite feature space of multi-channel fusion circuit breaker vibration signals;

[0014] 5) Input the composite features of the multi-channel vibration signals of the circuit breaker into the CFA-KELM fault diagnosis model to obtain the early fault state of the circuit breaker.

[0015] Furthermore, in step 1), the dual-tree complex wavelet transform algorithm is used to denoise the vibration signals. The steps of the dual-tree complex wavelet decomposition denoising algorithm include:

[0016] ψ(t) = ψ h (t) + iψ g (t)

[0017] where ψ(t) is a complex wavelet, ψ h (t) and ψ g (t) are two real wavelets, and i is the imaginary unit.

[0018] The wavelet coefficients and scaling coefficients of the real part tree transform above the dashed line are calculated by Eqs. (2) and (3)

[0019]

[0020] The wavelet coefficients and scaling coefficients of the imaginary part tree transform below the dashed line are calculated by Eqs. (4) and (5)

[0021]

[0022] Therefore, the wavelet coefficients and scaling coefficients of the dual-tree complex wavelet transform are

[0023]

[0024] The wavelet coefficients and scaling coefficients of the dual-tree complex wavelet transform are reconstructed by Eqs. (8) and (9).

[0025]

[0026] The reconstructed signal of the dual-tree complex wavelet transform can be expressed as

[0027]

[0028] Furthermore, the specific method of step 2) is as follows:

[0029] (1) Let the multi-channel signal vector be X(t), expressed as X(t) = [x1(t), x2(t),..., x c (t), ], and extract the predefined K multi-variable modulated oscillations from the multi-channel signal as follows:

[0030]

[0031] where u k is the k-th multi-channel joint modal component, representing the component with specific frequency characteristics extracted from the multi-channel signal. Each component is strongly coupled with its center frequency, and the total bandwidth of each mode is kept to a minimum;

[0032] (2) Construct the analytic signal of the modal component

[0033]

[0034] Denote the analytic signal of the \(k\)th mode component processed by Hilbert transform, Denote the original signal \(u\) k (t) of the Hilbert transform;

[0035] (3) The square of the second norm after frequency shifting is used as the estimated value of the bandwidth. Add the bandwidth estimates of all channels and modes to obtain the total cost function \(g\), and the formula is

[0036]

[0037] The cost function \(g\) represents the sum of the bandwidths of all mode components in all channels and is the object to be minimized in the MVMD algorithm; \(k\) is the mode component number, \(c\) is the channel number, and \(\omega\) k Denote the center frequency of the \(k\)th mode component, which is used to shift the spectrum of the analytic signal to the baseband;

[0038] (4) Represent the constrained optimization problem as the following constrained variational model:

[0039]

[0040] Construct an augmented function to improve the reconstruction accuracy:

[0041]

[0042] \(\lambda\) is the Lagrange operator; \(\beta\) is the penalty parameter;

[0043] The update methods of the mode component, center frequency, and Lagrange operator are as follows:

[0044]

[0045] Is the updated value of the mode component of the \(k\)th mode and the \(c\)th channel in the \((n + 1)\)th iteration in the frequency domain; Is the frequency domain representation of the original signal of the \(c\)th channel; Is the frequency domain representation of the Lagrange operator \(\lambda\) in the \(c\)th channel; \(\alpha\) is the penalty parameter, which controls the frequency domain filtering characteristics; \(\tau\) is the relaxation parameter, which controls the update step size;

[0046] The iteration stop condition in MVMD becomes:

[0047]

[0048] \(\varepsilon\) is the preset convergence threshold to judge whether the iteration stops;

[0049] Select the minimum envelope entropy as the index for optimizing the MVMD model parameters, introduce the chaotic firefly algorithm CFA to adaptively select the model mode number \(K\) and the penalty factor \(\alpha\), and obtain the optimal values of each parameter;

[0050] (1) Generate a chaotic firefly population

[0051] Adopt the Logistic chaotic mapping strategy. The expression of the Logistic mapping is:

[0052] x n+1 = l c x n (1 - x n )(22)

[0053] where x n ∈ [0, 1] represents the numerical solution of the nth iteration of the chaotic mapping, x0 is the initial iteration value, and l c is the control parameter of the chaotic mapping;

[0054] Assume that the required scale of the firefly population is m, and the spatial dimension of the firefly variable is d. The steps of the chaotic mapping are as follows:

[0055] (a) First, randomly generate a firefly individual y1 = (y 1,1 , y 1,2 ,..., y 1,d ) within the interval [0, 1];

[0056] (b) Let the firefly numbers be i = 1, 2,..., m, and the dimensions of the firefly variables be j = 1, 2,..., n. Generate the remaining m - 1 firefly individuals according to the Logistic mapping iteration to obtain the firefly chaotic sequence Y = (y1, y2,... y m ). The iteration calculation formula is:

[0057] y i+1,j = 4y i,j (1 - y i,j )(23)

[0058] (c) Map the generated firefly population chaotic sequence Y to the firefly search range [x min , x max according to the following formula to obtain the chaotic variable X = (x1, x2,..., x m ) of the firefly population.

[0059]

[0060] (2) Population brightness mechanism and attraction mechanism

[0061] (a) Brightness mechanism

[0062] The absolute brightness I of a firefly is proportional to the objective function f(x). The higher the brightness, the better the objective function. The absolute brightness of the ith firefly individual is Ii = f(x i ); The relative brightness refers to the brightness perceived on other individuals. The relative brightness is inversely proportional to the distance between two individuals. The farther the distance, the lower the brightness. The relative brightness formula is defined as:

[0063]

[0064] I0 is the initial light intensity of the firefly, γ is the light absorption coefficient, l i,j is the Euclidean distance between two firefly individuals, l i,j = ||x i - x j ||;

[0065] (b) Gravitational mechanism

[0066] Any two individuals among the randomly generated firefly populations will attract each other. The attraction between individuals is directly proportional to the relative brightness and inversely proportional to the distance. The definition of attraction is:

[0067]

[0068] where β0 represents the attraction without distance;

[0069] (3) Position update strategy

[0070] When there are N firefly individuals in the population in space, the fireflies will move according to the attraction, and their position update method is:

[0071]

[0072] In the formula, α is the step size, rand is a random number between [0,1], and β i,j is the attraction of firefly i to firefly j. After each iteration, the brightest firefly moves randomly, as shown in the formula:

[0073]

[0074] Furthermore, in step 3), when performing segmented time-frequency feature extraction, the method for constructing the segmented time-frequency feature matrix is as follows:

[0075] (1) Fusion of modal components segmentation:

[0076] The segmentation range is selected within the time t after the mutation point and divided into 6 segments;

[0077] (2) Segmented time-frequency feature extraction

[0078] (a) Energy entropy

[0079] The segmented energy entropy feature can characterize the energy complexity of each IMF in the time domain after the circuit breaker action starts. From this, the change of the complexity of each IMF over time can be analyzed. The segmented energy entropy E energy_entropy is calculated as follows:

[0080]

[0081] where x i represents the vibration amplitude, p i represents the ratio of the current position energy to the total energy, and beg and end represent the starting and ending serial numbers of the segmented signal;

[0082] (b) Kurtosis

[0083] Kurtosis is a measure of the fourth-order standard moment, used to describe the spikiness of a signal, and can be used to capture non-periodic mutations in the signal. In fault diagnosis, it is used to detect instantaneous faults in the signal. The formula for calculating kurtosis K is as follows:

[0084]

[0085] where N represents the number of data sequence samples in the segmented signal;

[0086] (c) Envelope entropy

[0087] Calculate the envelope of the segmented IMF obtained, and calculate the segmented envelope entropy for the segmented envelope. The formula for calculating the segmented envelope entropy is as follows:

[0088]

[0089] where p i is the probability distribution of the envelope signal;

[0090] (d) Spectral kurtosis

[0091] Spectral kurtosis represents the spikiness of the frequency spectrum. The formula for calculating spectral kurtosis is as follows:

[0092]

[0093] where N represents the number of data sequence samples in the segmented signal;

[0094] In the formula, f i represents the frequency, psd(f i ) represents the power spectral density at the current frequency, represents the average power spectral density;

[0095] (3) Construct a segmented time-frequency feature matrix: Define the energy entropy of a single modal component as e i , i = 1, 2,..., 6, and define the kurtosis as k i, where \(i = 1, 2, \cdots, 6\), the envelope entropy is defined as \(v\) i , where \(i = 1, 2, \cdots, 6\), the spectral kurtosis is defined as \(sk\) i , where \(i = 1, 2, \cdots, 6\), thus the vector form of the feature extracted from the \(t\)

[0096] -th modal component can be expressed as:

[0097] \(v\) IMFt = [ \(e\) t1 , \(\cdots\), \(e\) t6 , \(k\) t1 , \(\cdots\), \(k\) t6 , \(v\) t1 , \(\cdots\), \(v\) t6 , \(sk\) t1 , \(\cdots\), \(sk\) t6 T (33)

[0098] Suppose there are \(N\) vibration sampling channels, the matrix form of the segmented time-frequency features is:

[0099] \(V = [v\) IMF1 , \(v\) IMF2 , \(\cdots\), \(v\) IMFN T (34)

[0100] The dimension of the segmented time-frequency feature matrix is \([N\times24]\).

[0101] Furthermore, step 4) specifically includes:

[0102] (1) Global time-frequency feature extraction

[0103] Extract global time-domain features: standard deviation, root mean square, skewness, impulse factor, shape factor, and frequency-domain features: center frequency, root mean square frequency, frequency variance, frequency standard deviation;

[0104] (2) Construct the global time-frequency feature matrix

[0105] Suppose there are \(N\) vibration sampling channels, then the global time-frequency feature \(M\) of the \(k\)-th channel k = [ \(T\) k1 , \(T\) k2 , \(\cdots\), \(T\) k5 , \(F\) k1 , \(\cdots\), \(F\) k4 T , and the global time-frequency feature matrix is

[0106] \(M = [M_1, M_2, \cdots, M\) N T (35)

[0107] The dimension of the global time-frequency feature matrix is [N×9].

[0108] Furthermore, step 5) specifically includes:

[0109] Let the connection weights between the ELM hidden layer and the input layer and the output layer be ω and β respectively, set the bias of the hidden layer as b, and the activation function as f(x). The ELM model is expressed as:

[0110]

[0111] In the formula, ω i = [ω i1 , ω i2 ,..., ω in , x j = [x 1j , x 2j ,..., x n,j , T , t j = [t 1j , t 2j ,..., t m,j , T , N represents the number of output nodes, L represents the number of neurons in the hidden layer, m represents the dimension of the output t j , n represents the dimension of the input sample x j ; Let H be the output of the hidden layer, and the output of the network is T = [t1, t2,..., t n . The matrix form of ELM is expressed as:

[0112] Hβ = T T (37)

[0113] The output weight β is obtained through the least squares solution, and its solution is expressed as:

[0114] β = H T (HH T ) -1 T(38)

[0115] To avoid matrix inversion operations, the generalized inverse matrix H + of H is introduced, and the calculation formula is simplified to:

[0116]

[0117] The kernel extreme learning machine realizes the non-linear mapping from input to output by introducing a kernel function as the activation function in ELM; the selection of the kernel function needs to satisfy Mercer's theorem. Let K(x i , x j ) be expressed as the kernel function, and the Gaussian kernel function is used. Its formula is as follows:

[0118]

[0119] where σ is the bandwidth parameter of the Gaussian kernel;

[0120] Taking the F1-score metric F_score as the fitness function for CFA optimization to optimize the parameter selection of KELM; According to the true class of the samples and the predicted result classes of the fault diagnosis model, it is divided into four cases: true positive (TP), false positive (FP), true negative (TN), and false negative (FN). The F1-score metric F_score is the harmonic mean of precision and recall, and the calculation formula is as follows:

[0121]

[0122] where P is precision and R is recall;

[0123] The calculation formula for precision is as follows:

[0124]

[0125] The calculation formula for recall is as follows:

[0126]

[0127] Compared with the prior art, the beneficial effects of the present invention are: A diagnosis framework based on a composite feature space and an optimized classification model is proposed. The vibration signal is denoised by dual-tree complex wavelet transform, and its robustness in low signal-to-noise ratio scenarios is verified, effectively retaining the weak features of the fault signal; Further combining multivariate variational mode decomposition for frequency-domain alignment and fusion of multi-channel vibration-acoustic signals, solving the defects of mode mixing and insufficient utilization of multi-source information in traditional methods, and significantly improving the discrimination of fault features through the construction of a composite space of segmental time-frequency features and global features. The superiority of this method under complex working conditions is verified, providing a high-precision and low-cost solution for the real-time monitoring of early mechanical faults of circuit breakers. Description of the Drawings

[0128] Figure 1 is a schematic diagram of dual-tree complex wavelet decomposition;

[0129] Figure 2 is the vibration simulation signal with different signal-to-noise ratios;

[0130] Figure 3 is the noise reduction effect diagram of dual-tree complex wavelet decomposition for signals with different signal-to-noise ratios;

[0131] Figure 4 is the modal component diagram obtained by multivariate variational mode decomposition;

[0132] Figure 5 is the segmented diagram of the fused modal components;

[0133] Figure 6 is the segmented energy entropy statistical chart;

[0134] Figure 7 is the segmented kurtosis statistical chart;

[0135] Figure 8 is the segmented envelope entropy statistical chart;

[0136] Figure 9 is the spectral kurtosis statistical chart;

[0137] Figure 10 is the segmented time-frequency feature extraction flow chart;

[0138] Figure 11 is the fault diagnosis confusion matrix based on the composite feature space;

[0139] Figure 12 is the fault diagnosis flow chart based on the composite feature space of the breaker vibration signal. Specific implementation manners

[0140] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

[0141] First, collect the vibration signals of the multi-channel breaker operating mechanism, and construct a noise-free reference signal and a standard signal-to-noise ratio test signal, as shown in the appendix Figure 1 .

[0142] Regard the breaker vibration as a superposition of a series of exponentially decaying sine waves, and the formula is described as:

[0143]

[0144] where A represents the amplitude of the i-th exponentially decaying sine wave component; μ represents the decay factor; t i represents the starting time of the i-th exponentially decaying sine wave component; f i represents the vibration frequency of the i-th exponentially decaying sine wave component; ε(t - t i ) represents the step function.

[0145] Inject full-band random noise into the original signal to make its signal-to-noise ratio (SNR) 10 and 5, and simulate the noise contained in the measured signal. SNR represents the energy ratio of the signal to the mixed-in noise, and its calculation formula is:

[0146]

[0147] where P s , P n represent the signal power and the noise power respectively, and the calculation formulas are:

[0148]

[0149] In the formula, x(i) and s(i) respectively represent the noise-free signal of the reference and the noisy signal.

[0150] The vibration signal of the circuit breaker has strong non-periodic and non-stationary characteristics, and the characteristics are concentrated in the operating stage. However, the complex working conditions make the collected signal often contain noise, which interferes with the extraction of effective characteristics. Therefore, it is necessary to perform noise reduction preprocessing on the data before feature extraction to prevent noise from affecting the subsequent diagnosis accuracy. The dual-tree complex wavelet transform algorithm is selected to denoise the vibration signal. The schematic diagram of the dual-tree complex wavelet decomposition is as attached Figure 2 .

[0151] The steps of the dual-tree complex wavelet decomposition denoising algorithm include:

[0152] ψ(t) = ψ h (t) + iψ g (t)

[0153] where ψ(t) is the complex wavelet, and ψ h (t) and ψ g (t) are two real wavelets, and i is the imaginary unit.

[0154] The wavelet coefficients and scale coefficients of the real part tree transform above the dotted line are calculated by Equations (2) and (3)

[0155]

[0156] The wavelet coefficients and scale coefficients of the imaginary part tree transform below the dotted line are calculated by Equations (4) and (5)

[0157]

[0158] Therefore, the wavelet coefficients and scale coefficients of the dual-tree complex wavelet transform are

[0159]

[0160] The wavelet coefficients and scale coefficients of the dual-tree complex wavelet transform are reconstructed by Equations (8) and (9).

[0161]

[0162] The signal reconstructed by the dual-tree complex wavelet transform can be expressed as

[0163]

[0164] The noise reduction effect of the dual-tree complex wavelet decomposition on signals with different signal-to-noise ratios is as attached Figure 3 .

[0165] To comprehensively extract the fault characteristics of the circuit breaker operating mechanism, multi-channel vibration and sound signals are collected. For the decomposition of multi-channel homologous vibration signals, the multivariate variational mode decomposition (MVMD) algorithm is introduced to achieve frequency alignment of multi-channel signals and complete the fusion of multi-channel frequency-domain characteristics.

[0166] The multivariate variational mode decomposition model is obtained by optimizing the variational mode decomposition algorithm through multi-channel joint variational model construction, frequency consistency constraint, and alternating direction method of multipliers (ADMM) optimization solution;

[0167] Let the multi-channel signal vector be X(t), expressed as X(t) = [x1(t), x2(t),..., x c (t),], and the predefined K multi-variable modulated oscillations are extracted from the multi-channel signals as follows:

[0168]

[0169] where u k is the k-th multi-channel joint modal component, representing the component with specific frequency characteristics extracted from the multi-channel signals. Each component is strongly coupled with its center frequency, and the total bandwidth of each mode is kept to a minimum.

[0170] (2) Construct the analytic signal of the modal component

[0171]

[0172] represents the analytic signal of the k-th modal component processed by the Hilbert transform, represents the Hilbert transform of the original signal u k (t).

[0173] (3) The square of the two-norm after frequency shifting is the estimated value of the bandwidth. The estimated values of the bandwidths of all channels and modes are added to obtain the total cost function g, and the formula is

[0174]

[0175] The cost function g represents the sum of the bandwidths of all modal components in all channels and is the object to be minimized in the MVMD algorithm; k is the modal component number, c is the channel number, and ω k represents the center frequency of the k-th modal component, which is used to shift the spectrum of the analytic signal to the baseband.

[0176] (4) The constrained optimization problem is expressed as the following constrained variational model:

[0177]

[0178] Construct an augmented function to improve the reconstruction accuracy:

[0179]

[0180] λ is the Lagrange operator; β is the penalty parameter.

[0181] The update methods of the modal component, center frequency, and Lagrange operator are as follows:

[0182]

[0183] is the updated value of the modal component of the k-th mode and the c-th channel in the (n + 1)-th iteration in the frequency domain; is the frequency-domain representation of the original signal of the c-th channel; is the frequency-domain representation of the Lagrange operator λ in the c-th channel; α is the penalty parameter, controlling the frequency-domain filtering characteristics; τ is the relaxation parameter, controlling the update step size.

[0184] The iteration stop condition in MVMD becomes:

[0185]

[0186] ε is the preset convergence threshold, used to judge whether the iteration stops.

[0187] Select the minimum envelope entropy as the index for optimizing the MVMD model parameters, introduce the chaotic firefly algorithm CFA to adaptively select the model modal number K and the penalty factor α, and obtain the optimal values of each parameter.

[0188] (1) Generate a chaotic firefly population

[0189] Adopt the Logistic chaotic mapping strategy. The expression of the Logistic mapping is:

[0190] x n+1 = l c x n (1 - x n ) (22)

[0191] where x n ∈ [0, 1] represents the numerical solution of the n-th iteration of the chaotic mapping, x0 is the initial iteration value, and l c is the control parameter of the chaotic mapping. In this paper, its value is set to 4.

[0192] Assume that the required firefly population size is m, and the spatial dimension of the firefly variable is d. The chaotic mapping steps are as follows:

[0193] (a) First, randomly generate a firefly individual y1 = (y 1,1 , y 1,2 ,..., y1,d )

[0194] (b) Let the firefly numbers be \(i = 1, 2, \cdots, m\), and the dimension of the firefly variables be \(j = 1, 2, \cdots, n\). Generate the remaining \(m - 1\) firefly individuals by iterating the Logistic map to obtain the firefly chaotic sequence \(Y=(y_1,y_2,\cdots,y m ), and the iterative calculation formula is:

[0195] y i+1,j = 4y i,j (1 - y i,j ) (23)

[0196] (c) Map the generated firefly population chaotic sequence \(Y\) to the firefly search range \([x min , x max \) according to the following formula to obtain the chaotic variables \(X=(x_1,x_2,\cdots,x m ) of the firefly population.

[0197]

[0198] (2) Population brightness mechanism and attraction mechanism

[0199] (a) Brightness mechanism

[0200] The absolute brightness \(I\) of a firefly is proportional to the objective function \(f(x)\). The higher the brightness, the better the objective function. The absolute brightness of the \(i\)-th firefly individual is \(I i = f(x i ); The relative brightness refers to the brightness felt on other individuals. The relative brightness is inversely proportional to the distance between two individuals. The farther the distance, the lower the brightness. The relative brightness formula is defined as:

[0201]

[0202] \(I_0\) is the initial light intensity of this firefly, \(\gamma\) is the light absorption coefficient, \(l i,j is the Euclidean distance between two firefly individuals, \(l i,j = ||x i - x j ||.

[0203] (b) Gravitational mechanism

[0204] Any two individuals among the randomly generated firefly populations will attract each other. The attraction between individuals is proportional to the relative brightness and inversely proportional to the distance. The definition of attraction is:

[0205]

[0206] Among them, β0 represents the attraction force at zero distance.

[0207] (3) Position update strategy

[0208] When there are N firefly individuals in the space forming a population, the fireflies will move according to the attraction force, and their position update method is as follows:

[0209]

[0210] In the formula, α is the step size, rand is a random number between [0, 1], and β i,j is the attraction force of firefly i to firefly j. After each iteration, the brightest firefly moves randomly, as shown in the formula:

[0211]

[0212] Perform CFA-MVMD decomposition on the vibration signals of the multi-channel circuit breaker to obtain the fused modal components and their spectral characteristics. The modal component diagrams obtained by multivariate variational mode decomposition are as Figure 4 .

[0213] Modal decomposition separates the components in the frequency domain. These characteristics can reflect the characteristic differences between different modal components, that is, they represent the characteristics of the components with different frequencies after frequency domain decomposition. However, the vibration sound signal of the circuit breaker has strong non-periodicity, and the characteristics such as the amplitude and energy density of a single modal component in the time domain will also change greatly. Therefore, it is necessary to process the characteristic changes in the time domain inside the component. Perform segmentation on the obtained fused modal components so that the segmentation range covers the main vibration stage. The segmented diagram of the fused modal components is as Figure 5 . Extract four types of eigenvalue features, namely energy entropy, kurtosis, envelope entropy, and spectral kurtosis, for each segment. The statistical diagrams of each feature are as Figures 6 - 9 . Construct a segmented time-frequency feature matrix. The flow chart of segmented time-frequency feature extraction is as Figure 10 .

[0214] (1) Segmentation of fused modal components

[0215] It has been verified that 0.3 s after the signal mutation point on the dataset can cover the main vibration stage of one action. The segmentation range is selected within 0.3 s after the mutation point and divided into 6 segments.

[0216] (2) Extraction of segmented time-frequency features

[0217] (a) Energy entropy

[0218] The segmented energy entropy feature can characterize the energy complexity of each IMF in the time domain after the circuit breaker action starts, and thus the change of the complexity of each IMF over time can be analyzed. The segmented energy entropy E energy_entropy is calculated by the formula:

[0219]

[0220] where x i represents the vibration amplitude, p i represents the ratio of the current position energy to the total energy, and beg and end represent the starting sequence numbers of the segmented signals.

[0221] (b) Kurtosis

[0222] Kurtosis is a measure of the fourth-order standard moment, used to describe the sharpness of the signal peaks, and can be used to capture non-periodic mutations in the signal. It is used to detect instantaneous faults in the signal in fault diagnosis. The calculation formula for kurtosis K is as follows:

[0223]

[0224] where N represents the number of data sequence samples in the segmented signal.

[0225] (c) Envelope entropy

[0226] Calculate the envelope of the segmented IMF obtained, and calculate the segmented envelope entropy for the segmented envelopes. The calculation formula for the segmented envelope entropy is as follows:

[0227]

[0228] where p i is the probability distribution of the envelope signal.

[0229] (d) Spectral kurtosis

[0230] Spectral kurtosis can represent the sharpness of the frequency spectrum and can be approximately implemented as the kurtosis of the power spectral density. The calculation formula for spectral kurtosis is as follows:

[0231]

[0232] where N represents the number of data sequence samples in the segmented signal.

[0233] In the formula, f i represents the frequency, psd(f i ) represents the power spectral density of the current frequency, represents the average power spectral density.

[0234] (3) Construct the segmented time-frequency feature matrix

[0235] The method of adjacent single-category features is used. In the arrangement of single-category features, starting from the features extracted from the first time-domain segment to the features extracted from the last time-domain segment. Specifically, define the energy entropy of a single modal component as e i , i = 1, 2,..., 6, and define the kurtosis as k i, where \(i = 1, 2, \ldots, 6\), the envelope entropy is defined as \(v\) i , where \(i = 1, 2, \ldots, 6\), the spectral kurtosis is defined as \(sk\) i , where \(i = 1, 2, \ldots, 6\), the vector form of the feature extracted from the \(t\)-th modal component can be expressed as:

[0236] v IMFt =[e t1 , \ldots, e t6 , k t1 , \ldots, k t6 , v t1 , \ldots v t6 , sk t1 , \ldots sk t6 T (33)

[0237] Assume there are \(N\) vibration sampling channels, the matrix form of the segmented time-frequency features is:

[0238] V = [v IMF1 , v IMF2 , \ldots, v IMFN T (34)

[0239] The dimension of the segmented time-frequency feature matrix is \([N\times24]\).

[0240] The global features include two categories: time-domain features and frequency-domain features. Extract time-domain features such as standard deviation, root mean square, skewness, impulse factor, shape factor, and frequency-domain features such as center frequency, root mean square frequency, frequency variance, frequency standard deviation, and construct a global time-frequency feature matrix.

[0241] (1) Global time-frequency feature extraction

[0242] The global features include two categories: time-domain features and frequency-domain features. Extract 9-dimensional global features including 5-dimensional time-domain features and 4-dimensional frequency-domain features for the signals of each channel. The extracted time-domain features are: standard deviation, root mean square, skewness, impulse factor, shape factor, and the frequency-domain features are: center frequency, root mean square frequency, frequency variance, frequency standard deviation. The specific calculation formulas are as follows.

[0243]

[0244]

[0245] (2) Construct the global time-frequency feature matrix

[0246] Extract 9-dimensional global features including 5-dimensional time-domain features and 4-dimensional frequency-domain features for the signals of each channel. Assume there are \(N\) vibration sampling channels, then the global time-frequency feature \(M\) of the \(k\)-th channel​​k = [T k1 , T k2 ,... T k5 , F k1 ,..., F k4 T , the global time-frequency feature matrix is

[0247] M = [M1, M2,..., M N T (35)

[0248] The dimension of the global time-frequency feature matrix is [N × 9].

[0249] Construct a multi-channel fusion composite feature space of breaker vibration signals by combining segmental and global time-frequency feature extraction.

[0250] Input the multi-channel vibration signal composite features of the breaker into the CFA-KELM fault diagnosis model to obtain the early fault state of the breaker.

[0251] Let the connection weights of the hidden layer of ELM with the input layer and the output layer be ω and β respectively, set the bias of the hidden layer as b, and the activation function as f(x). The ELM model can be expressed as:[[]]

[0252]

[0253] In the formula, ω i = [ω i1 , ω i2 ,..., ω in , x j = [x 1j , x 2j ,..., x n,j T , t j = [t 1j , t 2j ,..., t m,j T , N represents the number of output nodes, L represents the number of neurons in the hidden layer, m represents the dimension of the output t j and n represents the dimension of the input sample x j . Let H be the output of the hidden layer, and the output of the network is T = [t1, t2,..., t n . The matrix form of ELM can be expressed as:[[]]

[0254] Hβ = T T (37)

[0255] The output weight β can be obtained through the least squares solution, and its solution can be expressed as:[[]]

[0256] β = H​​​​T (HH T ) -1 T(38)

[0257] To avoid the matrix inversion operation, the generalized inverse matrix H + of H is introduced, and the calculation formula is simplified to:

[0258]

[0259] The kernel extreme learning machine (KELM) realizes the non-linear mapping from input to output by introducing the kernel function as the activation function in ELM. The selection of the kernel function needs to satisfy Mercer's theorem. Representing K(x i , x j ) as the kernel function, the Gaussian kernel function is used in this paper, and its formula is as follows:

[0260]

[0261] where σ is the bandwidth parameter of the Gaussian kernel.

[0262] Furthermore, taking the F1 score index F_score as the fitness function for CFA optimization, it is used to optimize the parameter selection of KELM. According to the true category of the samples and the predicted result categories of the fault diagnosis model, it is divided into four cases: true positive (TP), false positive (FP), true negative (TN), and false negative (FN). The F1 score index F_score is the harmonic mean of the precision rate and the recall rate, and the calculation formula is as follows:

[0263]

[0264] where P is the precision rate and R is the recall rate;

[0265] The precision rate is the ratio of the number of correctly classified results to all classification results of a certain category in the prediction results of a certain category, and the calculation formula is as follows:

[0266]

[0267] The recall rate is the ratio of the number of correctly classified samples to all classification results of a certain actual category, and the calculation formula is as follows:

[0268]

[0269] Finally, the fault diagnosis confusion matrix based on the composite feature space is as Figure 11 , and the fault diagnosis flow chart based on the composite feature space of the breaker vibration signal is as Figure 12 .

Claims

1. A method for constructing a composite feature space of breaker vibration signals based on multi-channel fusion, characterized in that, The method includes the following steps: 1) Collect vibration signals of the multi-channel breaker operating mechanism, and construct a noise-free reference signal and a standard signal-to-noise ratio test signal; 2) Decompose the multi-channel breaker vibration signals by CFA-MVMD to obtain the fused modal components and their spectral characteristics; 3) Perform a segmentation operation on the obtained fused modal components, extract the segmented time-frequency characteristics, and construct a segmented time-frequency characteristic matrix; 4) Extract the global time-frequency characteristics and construct a global time-frequency characteristic matrix; Construct a composite feature space of the multi-channel fused breaker vibration signals by combining the segmented and global time-frequency feature extraction; 5) Input the composite features of the multi-channel breaker vibration signals into the CFA-KELM fault diagnosis model to obtain the early fault state of the breaker.

2. The method for constructing a composite feature space of breaker vibration signals based on multi-channel fusion according to claim 1, characterized in that In step 1), the dual-tree complex wavelet transform algorithm is used to denoise the vibration signals. The steps of the dual-tree complex wavelet decomposition denoising algorithm include: ψ(t) = ψ h (t) + iψ g (t) where ψ(t) is a complex wavelet, ψ h (t) and ψ g (t) are two real wavelets, and i is the imaginary unit. The wavelet coefficients and scale coefficients of the real part tree transform above the dotted line are calculated by equations (2) and (3) The wavelet coefficients and scale coefficients of the imaginary part tree transform below the dotted line are calculated by equations (4) and (5) Therefore, the wavelet coefficients and scale coefficients of the dual-tree complex wavelet transform are The wavelet coefficients and scale coefficients of the dual-tree complex wavelet transform are reconstructed by equations (8) and (9). The reconstructed signal of the dual-tree complex wavelet transform can be expressed as 3. The method for constructing a composite feature space of breaker vibration signals based on multi-channel fusion according to claim 1, wherein The specific method of step 2) is as follows: (1) Let the multi-channel signal vector be X(t), expressed as X(t) = [x1(t), x2(t),..., x c (t), ], and extract the predefined K multi-variable modulated oscillations from the multi-channel signal as follows: where u k is the k-th multi-channel joint modal component, representing the component with specific frequency characteristics extracted from the multi-channel signal. Each component is strongly coupled with its center frequency, and the total bandwidth of each mode is kept to a minimum; (2) Construct the analytical signal of the modal component denotes the analytic signal of the k-th modal component after Hilbert transform processing, denotes the original signal u k (t) of the Hilbert transform; (3) The square of the second norm after frequency shifting is the estimated value of the bandwidth. Add the bandwidth estimated values of all channels and modes to obtain the total cost function g. The formula is The cost function \(g\) represents the sum of the bandwidths of all modal components in all channels and is the object to be minimized in the MVMD algorithm; \(k\) is the modal component index, \(c\) is the channel index, and \(\omega\) k represents the center frequency of the \(k\)-th modal component, which is used to shift the spectrum of the analytic signal to the baseband; (4) Represent the constrained optimization problem as the following constrained variational model: Construct an augmented function to improve the reconstruction accuracy: λ is the Lagrangian operator; β is the penalty parameter; The update methods of the modal component, center frequency, and Lagrangian operator are as follows: is the updated value of the modal component of the k-th mode and the c-th channel at the (n + 1)-th iteration in the frequency domain; is the frequency-domain representation of the original signal of the c-th channel; is the frequency-domain representation of the Lagrange operator λ in the c-th channel; α is the penalty parameter that controls the frequency-domain filtering characteristics; τ is the relaxation parameter that controls the update step size; The iteration stop condition in MVMD becomes: ε is the preset convergence threshold to determine whether the iteration stops; Select the minimum envelope entropy as the index for optimizing the MVMD model parameters. Introduce the chaotic firefly algorithm CFA to adaptively select the model mode number K and the penalty factor α to obtain the optimal values of each parameter; (1) Generate a chaotic firefly population Adopt the Logistic chaotic mapping strategy. The expression of the Logistic mapping is: x n+1 = l c x n (1 - x n )(22) where x n ∈[0,1] represents the numerical solution of the n-th iteration of the chaotic map, x0 is the initial iteration value, and l c is the control parameter of the chaotic map; Assume that the required firefly population size is m, and the space dimension of the firefly variable is d. The chaotic mapping steps are as follows: (a) First, randomly generate a firefly individual \(y_1=(y 1,1 ,y 1,2 ,\cdots,y 1,d )\) within the interval \([0,1]\) with dimension \(d\); (b) Let the firefly numbers be \(i = 1, 2,\cdots,m\), and the dimension of the firefly variables be \(j = 1, 2,\cdots,n\). Generate the remaining \(m - 1\) firefly individuals by iterating the Logistic map to obtain the firefly chaotic sequence \(Y=(y_1,y_2,\cdots,y m )\), and the iterative calculation formula is: y i+1,j = 4y i,j (1 - y i,j ) (23) (c) Map the generated chaotic sequence Y of firefly population to the firefly search range [x min , x max according to the following formula, and obtain the chaotic variable X=(x1, x2,..., x m ) of the firefly population. (2) Population brightness mechanism and attraction mechanism (a) Brightness mechanism The absolute brightness I of a firefly is directly proportional to the objective function f(x). The higher the brightness, the better the objective function. The absolute brightness of the i-th firefly individual is I i = f(x i ); Relative Brightness refers to the brightness felt on other individuals. The relative brightness is inversely proportional to the distance between two individuals. The farther the distance, the lower the brightness. The relative brightness formula is defined as: $I_0$ is the initial light intensity of the firefly, $\gamma$ is the light absorption coefficient, and $l$ i,j is the Euclidean distance between two firefly individuals, and $l$ i,j $ = ||\mathbf{x}$ i $ - \mathbf{x}$ j $||$; (b) Gravitational mechanism Any two individuals among the randomly generated firefly populations will attract each other. The attraction between individuals is proportional to the relative brightness and inversely proportional to the distance. The definition of attraction is: where β0 represents the attraction without distance; (3) Position update strategy When there are N firefly individuals in the space forming a population, the fireflies will move according to the attraction. The position update method is: where α is the step size, rand is a random number between [0, 1], and β i,j is the attractiveness of firefly i to firefly j. After each iteration, the brightest firefly moves randomly, as shown in the equation:

4. The method for constructing a composite feature space of breaker vibration signals based on multi-channel fusion according to claim 1, wherein The method for extracting the segmented time-frequency characteristics and constructing the segmented time-frequency characteristic matrix in step 3) is as follows: (1) Segment the fused modal components: The segmentation range is selected within the time t after the mutation point and is divided into 6 segments; (2) Extraction of segmented time-frequency features (a) Energy entropy The segmented energy entropy feature can characterize the energy complexity of each IMF in the time domain after the breaker action starts. From this, the variation of the complexity of each IMF with time can be analyzed. The segmented energy entropy E energy_entropy is calculated by the formula: where x i represents the vibration amplitude, and p i represents the ratio of the current position energy to the total energy; beg and end represent the starting sequence numbers of the segmented signals. (b) Kurtosis Kurtosis is a measure of the fourth-order standard moment, used to describe the sharpness of a signal, and can be used to capture non-periodic mutations in a signal. It is used to detect instantaneous faults in a signal during fault diagnosis. The formula for calculating kurtosis K is as follows: where N represents the number of data sequence samples in the segmented signal; (c) Envelope entropy Calculate the envelope of the segmented IMF obtained, and calculate the segmented envelope entropy for the segmented envelope. The formula for the segmented envelope entropy is as follows: where p i is the probability distribution of the envelope signal; (d) Spectral kurtosis Spectral kurtosis represents the sharpness of the frequency spectrum. The formula for calculating spectral kurtosis is as follows: where N represents the number of data sequence samples in the segmented signal; where f i represents the frequency, and psd(f i ) represents the power spectral density at the current frequency, represents the average power spectral density; (3)Construct a segmented time-frequency feature matrix: Define the energy entropy of a single modal component as e i , i = 1, 2, ..., 6, and the kurtosis is defined as k i , i = 1, 2, ..., 6, and the envelope entropy is defined as v i , i = 1, 2, ..., 6, and the spectral kurtosis is defined as sk i , i = 1, 2, ..., 6. Thus, the vector form of the features extracted from the t-th modal component can be expressed as: v IMFt = [e t1 ,..., e t6 , k t1 ,..., k t6 , v t1 ,... v t6 , sk t1 ,... sk t6 T (33)​ Assume there are N vibration sampling channels. The matrix form of the segmented time-frequency features is: V = [v IMF1 , v IMF2 ,..., v IMFN T (34)​ The dimension of the segmented time-frequency feature matrix is [N×24].

5. The method for constructing a composite feature space of breaker vibration signals based on multi-channel fusion according to claim 1, wherein Step 4) Specifically includes: (1) Extraction of global time-frequency features Extract global time-domain features: standard deviation, root mean square, skewness, impulse factor, shape factor, and frequency-domain features: center frequency, root mean square frequency, frequency variance, frequency standard deviation; (2) Construction of the global time-frequency feature matrix Suppose there are N vibration sampling channels, then the global time-frequency feature M of the k-th channel k = [T k1 , T k2 ,... T k5 , F k1 ,..., F k4 T , and the global time-frequency feature matrix is​ M = [M1, M2,..., M N T (35)​ The dimension of the global time-frequency feature matrix is [N×9].

6. The method for constructing a composite feature space of breaker vibration signals based on multi-channel fusion according to claim 1, characterized in that Step 5) Specifically includes: Let the connection weights between the hidden layer and the input layer and output layer of ELM be ω and β respectively. Set the bias of the hidden layer as b and the activation function as f(x). The ELM model is expressed as: where ω i = [ω i1 , ω i2 ,..., ω in , x j = [x 1j , x 2j ,..., x n,j T , t j = [t 1j , t 2j ,..., t m,j T , N represents the number of output nodes, L represents the number of neurons in the hidden layer, m represents the dimension of the output t j , n represents the dimension of the input sample x j ; Let H be the output of the hidden layer, the output of the network is T = [t1, t2,..., t n , and the matrix form of ELM is expressed as:​​ Hβ = T T (37) Obtain the output weight β through the least squares solution, and its solution is expressed as: β = H T (HH T ) -1 T (38) To avoid matrix inversion operations, the generalized inverse matrix H of H is introduced + , and the calculation formula is simplified to: The kernel extreme learning machine realizes the non-linear mapping from input to output by introducing a kernel function as the activation function in ELM; the selection of the kernel function needs to satisfy Mercer's theorem. Denote K(x i , x j ) as the kernel function, and use the Gaussian kernel function, whose formula is as follows: where σ is the bandwidth parameter of the Gaussian kernel; Use the F1 score index F_score as the fitness function for CFA optimization to optimize the parameter selection of KELM; According to the true category of the sample and the predicted result category of the fault diagnosis model, it is divided into four cases: true positive TP, false positive FP, true negative TN, and false negative FN. The F1 score index F_score is the harmonic mean of precision and recall. The calculation formula is as follows: where P is precision and R is recall; The formula for calculating precision is as follows: The formula for calculating recall is as follows:

Images