Wideband measurement method and system for power system based on wavelet packet and SVMD

Abstract

The invention relates to a power system broadband measurement method based on combination of a wavelet packet and SVMD. The method comprises the following steps: acquiring frequency information corresponding to circuit signals of all measurement points in a circuit system; performing wavelet packet decomposition on the frequency information, and determining an optimal wavelet basis for broadband measurement of the power system; performing multilayer decomposition on the frequency information according to the optimal wavelet basis to obtain a plurality of wavelet frequency bands; performing component extraction on the wavelet frequency band to obtain a harmonic component corresponding to the wavelet frequency band; de-noising the harmonic components by using a wavelet coefficient correlation algorithm to obtain de-noised harmonics; decomposing the noise reduction harmonic wave according to an EMD algorithm to obtain a signal residual error meeting a preset stop condition; and calculating to obtain measurement and calculation data corresponding to each harmonic waveform according to the signal residual error. According to the technical scheme, the time-frequency characteristics of the signals can be captured more accurately, and the method has better performance on non-stationary signals and signals with complex frequency structures.

Description

Technical Field

[0001] The present invention relates to the technical field of electronic measurement in power systems, and particularly to a wide-band measurement method and system for power systems based on the combination of wavelet packet and SVMD. Background Art

[0002] For an ideal three-phase power system, the basic parameters that a measuring device needs to measure are the amplitude and phase angle of the power frequency voltage. However, usually, the voltage and current waveforms obtained at each measurement point in the power grid are not three-phase symmetric sine waves, which is related to the devices connected thereto. If the power supply is connected by a power electronic device and the load is non-linear, a large amount of harmonics will be injected into the system, causing distortion of the voltage and current waveforms. With the wide application of non-linear loads and power electronic devices, the monitoring of high-order harmonic information of wide-band signals has become an important direction in the research of current new power systems. Compared with the traditional wavelet transform, the wavelet packet transform can more comprehensively describe the time-frequency characteristics of signals. However, although the wavelet decomposition algorithm can effectively divide the frequency band of wide-band signals, the algorithm cannot complete the extraction and detection of specific frequency harmonics, thus affecting the accuracy and reliability of the measurement results.

[0003] The present invention will introduce the SVMD (Successive Variational Mode Decomposition, SVMD) algorithm to further process the wavelet decomposition frequency band, and combine the adaptive extraction characteristics of the introduced algorithm with the frequency band decomposition characteristics of the wavelet algorithm to fully extract each frequency component, thereby realizing the effective measurement of wide-band signals. Summary of the Invention

[0004] The present invention aims to solve at least one of the technical problems existing in the prior art or related technologies.

[0005] To this end, the object of the present invention is to provide a wide-band measurement method and system for power systems based on the combination of wavelet packet and SVMD, which can more accurately capture the time-frequency characteristics of signals and have better performance for non-stationary signals and signals with complex frequency structures.

[0006] To achieve the above object, the technical solution of the first aspect of the present invention provides a wide-band measurement method for power systems based on the combination of wavelet packet and SVMD, including the following steps:

[0007] Step 1: Obtain the frequency information corresponding to the circuit signals at each measurement point in the circuit system;

[0008] Step 2: Perform wavelet packet decomposition on the frequency information to determine the optimal wavelet basis for wide-band measurement of the power system;

[0009] Step 3: Perform multi-level decomposition on the frequency information according to the optimal wavelet basis to obtain multiple wavelet frequency bands;

[0010] Step 4: Extract components from the wavelet frequency bands, and adaptively extract harmonic components of specific frequencies using the SVMD algorithm to obtain the harmonic components corresponding to the wavelet frequency bands;

[0011] Step 6: Denoise the harmonic components using the wavelet coefficient correlation algorithm to obtain denoised harmonics;

[0012] Step 9: Decompose the denoised harmonics according to the EMD algorithm to obtain a signal residual that meets the preset stop condition;

[0013] Step 12: Calculate the measurement data corresponding to each harmonic waveform according to the signal residual.

[0014] In the above technical solution, preferably, determining the optimal wavelet basis for wide-band measurement of the power system includes:

[0015] Using the Mallat decomposition algorithm, determine the decomposition relationship between the standard orthogonal bases and the decomposition and synthesis relationship of the subspace coordinates corresponding to the decomposition relationship between the standard orthogonal bases;

[0016]

[0017] In equations (1)-(2), t is time, φ(t) is the wavelet function, ψ(t) is the wavelet function of the subspace coordinates, i is the number of nodes, j is the decomposition level, k is the number of points in each layer, Z is an integer, m is an integer in Z, ∑ is the summation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter, and the coordinate sequences {h m ; m∈Z}, {g m ; m∈Z} are called the low-pass filter and the high-pass filter.

[0018] Perform transform space decomposition on the decomposition relationship between the standard orthogonal bases according to the preset conditions to obtain the standard orthogonal basis of a closed subspace on the finite function space;

[0019]

[0020]

[0021] In equations (3)-(4), j is the decomposition level, k is the number of points in each layer, Z is an integer, m is an integer in Z, α jk is the coefficient sequence, c j,k is the scaling coefficient, d j,k is the wavelet coefficient, hm-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

[0022] According to the properties of the filter itself: s 2i+1 (t) = ψ j,k (t), the orthonormal bases s2i(t) and s2i+1(t) are deduced to form two mutually orthogonal orthonormal systems and the subspace relationship of the corresponding orthonormal systems after integer translation;

[0023]

[0024] (5)-(7) In the formula, S is the original signal of the root node representing the output of a filter, S i j represents a closed subspace on L2(R), i is the number of nodes, j is the decomposition layer number, k is the number of points per layer, Z is an integer, m is an integer in Z, ∑ is the summation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter.

[0025] According to the orthonormal system and its subspace relationship, the decomposition result of the corresponding wavelet packet and the optimal wavelet basis are obtained;

[0026] The decomposition result is expressed as:

[0027]

[0028] The optimal wavelet basis is expressed as:

[0029]

[0030] (8)-(9) In the formula, α j，k is the coefficient sequence, d j,k is the wavelet coefficient, i is the number of nodes, j is the decomposition layer number, k is the number of points per layer, Z is an integer, m is an integer in Z, ∑ is the summation symbol, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

[0031] In the above technical solution, preferably, the frequency information is decomposed into multiple layers according to the optimal wavelet basis, including:

[0032] Perform wavelet transform on the optimal wavelet basis to obtain the representation of a single - period function for one period of the harmonic waveform:

[0033]

[0034] where \(f(t)\) is a periodic function, \(e\) ikt is the basis, \(k\) is the number of points in each layer, \(\alpha\) k can be expressed as the coordinate sequence of the periodic function \(f(t)\) in the basis \(e\) ikt .

[0035] Obtain the function representation of the entire harmonic waveform from the single - period function representation:

[0036]

[0037] where \(G(t\) o , \(\omega)\) represents the short - time Fourier transform, \(t\) is time, \(\omega\) is the angular frequency, \(f(t)\) is the waveform signal to be analyzed; \(g(t)\) is the window function that restricts the analysis time range, \(t_0\) is the specified time point indicating the magnitude of the component with frequency \(\omega\) of the signal \(f(t)\) near the time point \(t_0\) and within the range determined by the function \(g(t)\), \(e\) -jωt is the basis function.

[0038] Perform continuous wavelet transform on the function representation of the entire harmonic waveform according to the dynamic window characteristics to obtain the continuous wavelet dependent on the scale parameter and the function representation of the corresponding continuous wavelet;

[0039] The dynamic window characteristics are expressed as: \(\psi(t)\in L\) 2 (\mathbb{R})

[0040]

[0041] The continuous wavelet dependent on the scale parameter is expressed as:

[0042]

[0043] In equations (12)-(14), \(C\) ψ is the wavelet packet energy, \(\psi(t)\) is the wavelet basis function, is the Fourier transform of the wavelet basis function, \(\omega\) is the angular frequency, \(\mathbb{R}^*\) represents the set of all non - zero real numbers, \(\mathbb{R}\) represents the set of all real numbers, and at the same time, scale parameter \(a\) and displacement parameter \(b\) are set, \(L\) is the space of all signals, \(e\) -jωt is the basis function, and the function representation of the corresponding continuous wavelet:

[0044]

[0045] The window of the wavelet transform is a rectangular window centered at the coordinates (b, ±ω0 / a), with a time-domain window width of aΔψ and a frequency-domain window width of Δψ / a, where a is the scale parameter and b is the displacement parameter. is the complex conjugate of ψ(t).

[0046] Performing an inverse transform on the continuous wavelet and its corresponding function representation yields the inverse function representation corresponding to the continuous wavelet:

[0047]

[0048] where is the inverse of the wavelet packet energy. The wavelet transform shows how the signal f(t) changes as the scale parameter a varies near the time point b; discretizing the inverse function gives multiple wavelet frequency bands.

[0049] In the above technical solution, preferably, the discretization process of the inverse function includes:

[0050] Defining the preset function representation that constitutes the orthonormal basis of the corresponding wavelet frequency band:

[0051]

[0052] Obtaining the linear representation of the corresponding wavelet frequency band and the discretized inverse function according to the preset function representation:

[0053]

[0054] where ψ(t) ∈ L 2 (R) forms an orthonormal basis, so ψ(t) is an orthogonal wavelet;

[0055] Obtaining the wavelet frequency band according to the coefficient sequence of the preset orthonormal basis:

[0056]

[0057] In equations (17)-(19), the coefficient sequence α j，k represents the continuous wavelet transform of f(t) when the scale parameter a = 2 -j , and the displacement parameter b = 2 -j k, is the complex conjugate of ψ(t), ψ(t) is the wavelet basis function, i is the number of nodes, j is the decomposition level, k is the number of points per layer, and W f is the function of the continuous wavelet.

[0058] In the above technical solution, preferably, for component extraction of the wavelet frequency band, the SVMD algorithm is used to adaptively extract harmonic components of a specific frequency, and obtaining the harmonic components corresponding to the wavelet frequency band includes:

[0059] First, the SVMD algorithm is used to perform successive mode decomposition on the wavelet frequency band to obtain each component signal u k (t). The expression of the component signal u k (t) is as follows:

[0060]

[0061] In the formula, u k (t) is the k-th mode component, U k (t) is the amplitude function of the k-th component signal, φ k (t) is the phase function, then the frequency ω k (t) = dφ k (t) / dt, and all the above functions have non-negative characteristics;

[0062] Among them, the SVMD algorithm satisfies the following constrained variational problem:

[0063]

[0064] In the formula, {ω k} is the column vector of the center frequencies of each component; δ(t) is the unit impulse signal, is the derivative operator with respect to the time variable t, j is the imaginary unit, f is the original input signal, that is, the target signal to be decomposed, and the equality condition ensures that the sum of all components is the original function;

[0065] To solve the constrained variational problem, an unconstrained variational problem with a penalty function is constructed by introducing the Lagrange multiplier λ and the penalty factor α:

[0066]

[0067] (Equation (22)) is defined as the augmented Lagrangian function L({u k},{ω k},λ), and it is solved by an iterative method. The iterative formula is expressed as:

[0068]

[0069] Each time an iteration is completed, it is calculated whether the current mode component satisfies the following convergence condition:

[0070]

[0071] In equations (22)-(24), f(t) is the original input signal, λ k (t) is the Lagrange multiplier of the k-th mode constraint, is the updated result of the k-th mode component in the frequency domain in the (n + 1)-th iteration, is the Fourier transform of the original signal f(t). is the sum of the spectra of the remaining modes except the k-th mode obtained in the current iteration. is the representation of the Lagrange multiplier in the frequency domain, ω k is the central frequency of the k-th mode component, ω is the angular frequency. is the updated value of the Lagrange multiplier at the (n + 1)-th iteration, τ is the step size coefficient, K is the total number of modes to be decomposed, n is the number of iterations, k is the index of the mode component. is the updated value of the central frequency corresponding to the k-th mode component at the (n + 1)-th iteration. is the representation of the k-th mode in the frequency domain at the n-th iteration. is the representation of the k-th mode in the frequency domain at the (n + 1)-th iteration, ε is the convergence threshold, a small positive number. When the updates of all modes are small enough, the algorithm satisfies the convergence condition and the iteration can be stopped.

[0072] When the decomposed mode components meet the given accuracy, extract K components. The signal f(t) can be successively decomposed by the SVMD algorithm in the following way:

[0073]

[0074] In the formula, u K (t) is the K-th extracted mode component, x r (t) is the residual component, that is, the remaining part of the original signal after extracting the K-th component, x u (t) is the noise term or the residual term not yet extracted.

[0075] When extracting the K-th component u K (t), construct the filter impulse responses β K (t) and β n (t) to form a new constraint criterion for SVMD:

[0076]

[0077] In the formula, Q1 is the magnitude of the energy obtained after applying the filter β K (t) to the remaining signal. The smaller Q1 is, the less ineffective extraction of the remaining signal by the new filter, and it can be better distinguished from the remaining part. Q2 is the degree of overlap between the new component and the previously extracted components. The smaller Q2 is, the lower the overlap or aliasing degree between the new component u K (t) and the previously extracted components, and the better the distinguishability between components. β K (t) is to avoid the frequency of the extracted component u K (t) and the residual component xr (t) Overlapped filter impulse response, β n (t) To effectively distinguish adjacent extracted components u K (t) and u K-1 (t)'s filter impulse response;

[0078] The above filter is the key constraint condition to ensure the accurate extraction of each component by the SVMD algorithm. The expression of the filter is:

[0079]

[0080] Therefore, formula (21) is transformed into:

[0081]

[0082] In equations (27)-(28), is the weight function of the k-th modal component, which is used to adjust the weight of the signal in the frequency domain so that the mode is mainly concentrated around its center frequency ω K Nearby, when ω is close to ω K When,[[]] Obtains a large value, indicating that the signal energy is mainly concentrated around ω K Nearby, when ω is far from ω K When,[[]] Decays rapidly, so that the mode will not spread in the frequency range far from ω K ; Is the filtering weight function used to adjust the modal components in the frequency domain, which is used to control the frequency distribution of the extracted modal components to reduce the frequency aliasing between adjacent modes. When ω is close to ω n When,[[]] Obtains a large value, indicating that the mode has a large weight near ω n Nearby, maintaining the main frequency components. When ω is far from ω n When,[[]] Decreases rapidly, so that the mode will not expand to other frequency regions, thereby reducing the spectral overlap between adjacent modes; ω K Is the center frequency of the k-th mode, ω n Is the center frequency of the n-th mode, and α is the penalty factor, which is used to adjust the influence of the filter bandwidth;

[0083] Then, the NRBO optimization algorithm is used to optimize the penalty factor α, and the minimum envelope entropy is used as the fitness function:

[0084]

[0085] In the formula, p iFor the probability of calculating the envelope amplitude distribution of discrete components, S(i) is the envelope entropy of each component, N is a natural number, and i is 0, 1, 2...;

[0086] Through the above SVMD and NRBO optimization steps, accurate harmonic components are extracted from the wavelet frequency band, and the wavelet frequency band is subjected to component extraction according to the harmonic components.

[0087] The technical solution of the second aspect of the present invention provides a power system wide-frequency measurement system based on the combination of wavelet packet and SVMD, including:

[0088] An acquisition module, configured to acquire the frequency information corresponding to the circuit signals at each measurement point in the circuit system;

[0089] A wavelet decomposition module, configured to perform wavelet packet decomposition on the frequency information to determine the optimal wavelet basis for power system wide-frequency measurement;

[0090] A multi-layer decomposition module, configured to perform multi-layer decomposition on the frequency information according to the optimal wavelet basis to obtain multiple wavelet frequency bands;

[0091] An SVMD principal component extraction module, configured to perform modal decomposition on the basis of the wavelet frequency band by using the SVMD algorithm and sequentially extract harmonic components;

[0092] An NRBO optimization module, configured to optimize the penalty parameter α of the SVMD algorithm and use the minimum envelope entropy as the objective function to reduce harmonic aliasing;

[0093] A noise reduction module, configured to perform noise reduction processing on the harmonic components by using the wavelet coefficient correlation algorithm to obtain noise-reduced harmonics;

[0094] A harmonic decomposition module, configured to decompose the noise-reduced harmonics according to the EMD algorithm to obtain a signal residual that meets the preset stop condition;

[0095] A measurement module, configured to calculate measurement data corresponding to each harmonic waveform according to the signal residual.

[0096] In the above technical solution, preferably, the wavelet decomposition module includes:

[0097] An orthogonal basis decomposition unit, configured to use the Mallat decomposition algorithm to determine the decomposition relationship between standard orthogonal bases and the decomposition and synthesis relationship of the subspace coordinates corresponding to the decomposition relationship between standard orthogonal bases;

[0098]

[0099] In equations (1)-(2), t is time, φ(t) is a wavelet function, ψ(t) is a wavelet function of subspace coordinates, i is the number of nodes, j is the decomposition level, k is the number of points per layer, Z is an integer, and m is an integer in Z. ∑ is the summation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter, and the coordinate sequences {h m ; m ∈ Z} and {g m ; m ∈ Z} are called the low-pass filter and the high-pass filter

[0100] space transformation unit, which is set to perform space decomposition on the decomposition relationship between orthonormal bases according to preset conditions to obtain the orthonormal basis of a closed subspace on the finite function space;

[0101]

[0102] In equations (3)-(4), j is the decomposition level, k is the number of points per layer, Z is an integer, m is an integer in Z, α jk is the coefficient sequence, c j,k is the scaling coefficient, d j,k is the wavelet coefficient, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

[0103] Subspace decomposition unit, which is set to, according to the properties of the filter itself: s 2i+1 (t) = ψ j，k (t), deduce that the orthonormal bases s2i(t) and s2i + 1(t) form two mutually orthogonal orthonormal systems after integer translation and the subspace relationship of the corresponding orthonormal systems;

[0104]

[0105] In equations (5)-(7), S is the original signal of the root node representing the output of a filter, and S i j represents a closed subspace on L2(R), i is the number of nodes, j is the decomposition level, k is the number of points per layer, Z is an integer, and m is an integer in Z. ∑ is the summation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter.

[0106] The optimal wavelet basis generation unit is configured to obtain the decomposition result of the corresponding wavelet packet and the optimal wavelet basis according to the orthonormal system and its subspace relationship;

[0107] The decomposition result is expressed as:

[0108]

[0109] The optimal wavelet basis is expressed as:

[0110]

[0111] In equations (8)-(9), α j，k is the coefficient sequence, d j,k is the wavelet coefficient, i is the number of nodes, j is the decomposition level, k is the number of points per layer, Z is an integer, m is an integer in Z, ∑ is the summation symbol, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

[0112] In the above technical solution, preferably, the multi-layer decomposition module includes:

[0113] The wavelet transform unit is configured to perform wavelet transform on the optimal wavelet basis to obtain a single-period function representation of one period of the harmonic waveform:

[0114]

[0115] In the formula, f(t) is a periodic function, e ikt is the basis, k is the number of points per layer, α k can be expressed as the coordinate sequence of the periodic function f(t) in the basis e ikt ;

[0116] The waveform generation unit is configured to obtain a function representation of the entire harmonic waveform according to the single-period function representation:

[0117]

[0118] In the formula, G(t o , w) represents the short-time Fourier transform, t is time, w is the angular frequency, f(t) is the waveform signal to be analyzed; g(t) is the window function that restricts the analysis time range, t0 is the specified time point indicating the magnitude of the component with frequency ω of the signal f(t) near the time point t0 and within the range determined by the function g(t), e -jωt is the basis function.

[0119] A continuous wavelet generation unit, which is configured to perform continuous wavelet transform on the functional representation of the entire harmonic waveform according to dynamic window characteristics, to obtain continuous wavelets dependent on a scale parameter and the functional representation corresponding to the continuous wavelets;

[0120] The dynamic window characteristics are expressed as: ψ(t) ∈ L 2 (R)

[0121]

[0122] The continuous wavelets dependent on the scale parameter are expressed as:

[0123]

[0124] In equations (12)-(14), C ψ is the wavelet packet energy, ψ(t) is the wavelet basis function, is the Fourier transform of the wavelet basis function, w is the angular frequency, R* represents all non-zero real numbers, R represents all real numbers, and at the same time the scale parameter a and the displacement parameter b are set, L is the entire signal space, e -jωt is the basis function. The functional representation corresponding to the continuous wavelets:

[0125]

[0126] The window of the wavelet transform is a rectangular window centered at the coordinates (b, ±ω0 / a), with a time-domain window width of aΔψ and a frequency-domain window width of Δψ / a.

[0127] In the formula, a is the scale parameter and b is the displacement parameter, is the complex conjugate of ψ(t).

[0128] An inverse transform unit, which is configured to perform an inverse transform on the continuous wavelets and their corresponding functional representations to obtain the inverse functional representation corresponding to the continuous wavelets:

[0129]

[0130] In the formula is the inverse of the wavelet packet energy. The wavelet transform shows the variation of the signal f(t) near the time point b as the scale parameter a continuously changes; A discrete processing unit, which is configured to perform discretization processing on the inverse function to obtain multiple wavelet frequency bands.

[0131] In the above technical solution, preferably, the discrete processing unit includes:

[0132] A wavelet frequency band generation sub-unit, which is configured to define a preset functional representation that constitutes an orthonormal basis corresponding to the wavelet frequency band:

[0133]

[0134] A band linear function generation subunit, which is configured to obtain a linear representation of the corresponding wavelet band and the discretized inverse function according to a preset function representation:

[0135]

[0136] Among them, the orthonormal basis that constitutes ψ(t) ∈ L 2 (R), then ψ(t) is an orthogonal wavelet;

[0137] An orthogonal coefficient generation subunit, which is configured to obtain a wavelet band according to a coefficient sequence of a preset orthonormal basis:

[0138]

[0139] In equations (17)-(19), the coefficient sequence α j，k represents the continuous wavelet transform of f(t) when the scale parameter a = 2 -j , and the displacement parameter b = 2 -j k, is the complex conjugate of ψ(t), ψ(t) is the wavelet basis function, i is the number of nodes, j is the decomposition level, k is the number of points per layer, and W f is the function of the continuous wavelet.

[0140] In the above technical solution, preferably, the SVMD principal component extraction module includes:

[0141] A variational constraint unit: used to construct a constrained variational problem of SVMD for the input wavelet decomposition band signal;

[0142] An alternating multiplier iteration unit: used to iteratively solve each modal component by using the alternating direction multiplier algorithm based on the introduced Lagrange multiplier λ and penalty factor α;

[0143] A parameter optimization unit: used to adaptively optimize the penalty factor α by using the NRBO optimization algorithm, and use the minimum envelope entropy as the objective function;

[0144] A harmonic extraction unit: used to extract the harmonic components corresponding to the frequency band according to each modal component decomposed by SVMD and output them to the noise reduction module.

[0145] Compared with the prior art, the advantages of the power system wide-band measurement method and system based on the combination of wavelet packet and SVMD provided by the present invention are as follows: The wavelet packet transform has higher resolution and richer frequency information, can capture the time-frequency characteristics of signals more accurately, and has better performance for non-stationary signals and signals with complex frequency structures. At the same time, it has greater flexibility and adjustability, and can decompose signals more deeply or more roughly according to needs, so as to obtain time-frequency representations at different scales.

[0146] The basic wavelet packet decomposition algorithm can only complete frequency band division and cannot decompose each component in detail. By combining the frequency band decomposition function of the wavelet packet algorithm and the adaptive extraction function of SVMD, the difficult problems of wavelet basis and decomposition layer selection can be effectively solved.

[0147] The method proposed by the present invention based on the combination of these two algorithms not only reduces the decomposition layer of wavelets, but also significantly improves the detection efficiency of wide-band signals and adaptively completes the extraction of the main components in each frequency band. The improved wide-band measurement algorithm can effectively extract each harmonic component. Through analysis and calculation, the maximum relative error of amplitude detection is 1.45%, and the maximum error of phase angle detection is about 0.1°, verifying the rationality and accuracy of the algorithm for processing wide-band and dynamic signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0148] The above and / or additional aspects and advantages of the present invention will become obvious and easy to understand from the description of the embodiments in conjunction with the following drawings, where:

[0149] Figure 1 Shows the flowchart of the method involved in the embodiment of the present invention;

[0150] Figure 2 Shows the flowchart of step S2 involved in the embodiment of the present invention;

[0151] Figure 3 Shows the flowchart of step S3 involved in the embodiment of the present invention;

[0152] Figure 4 Shows the flowchart of step S35 involved in the embodiment of the present invention;

[0153] Figure 5 Shows the flowchart of step S4 involved in the embodiment of the present invention;

[0154] Figure 6 Shows the structural block diagram of the system involved in the embodiment of the present invention;

[0155] Figure 7 Shows the structural block diagram of the wavelet decomposition module involved in the embodiment of the present invention;

[0156] Figure 8 The structural block diagram of the multi-layer decomposition module involved in the embodiments of the present invention is shown;

[0157] Figure 9 The structural block diagram of the discrete processing unit involved in the embodiments of the present invention is shown;

[0158] Figure 10 The structural block diagram of the component extraction module involved in the embodiments of the present invention is shown;

[0159] Figure 11 The schematic diagram of single-operation decomposition involved in the embodiments of the present invention is shown;

[0160] Figure 12 The schematic diagram of wavelet decomposition involved in the embodiments of the present invention is shown;

[0161] Figure 13 The schematic diagram of wavelet packet tree decomposition involved in the embodiments of the present invention is shown;

[0162] Figure 14 The wavelet packet decomposition diagram of the signal model involved in the embodiments of the present invention is shown;

[0163] Figure 15 The effect diagram of main component extraction of the decomposition frequency band involved in the embodiments of the present invention is shown;

[0164] Figure 16 The flow chart of correlation noise reduction processing involved in the embodiments of the present invention is shown;

[0165] Figure 17 The effect diagram of adaptive EMD decomposition involved in the embodiments of the present invention is shown;

[0166] Figure 18 The curve diagram of the permutation entropy value of each component involved in the embodiments of the present invention is shown;

[0167] Figure 19 The comparison diagram of noise reduction effects involved in the embodiments of the present invention is shown. Detailed implementation manners

[0168] In order to be able to more clearly understand the above-mentioned objects, features and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific implementation manners. It should be noted that, without conflict, the embodiments of the present application and the features in the embodiments can be combined with each other.

[0169] Many specific details are set forth in the following description in order to provide a thorough understanding of the present invention. However, the present invention can also be implemented in other ways different from those described herein. Therefore, the protection scope of the present invention is not limited to the limitations of the specific embodiments disclosed below.

[0170] Such as Figure 1As shown, a wide - band measurement method for power systems based on the combination of wavelet packet and SVMD according to an embodiment of the present invention includes the following steps:

[0171] S1, obtain the frequency information corresponding to the circuit signals at each measurement point in the circuit system;

[0172] S2, perform wavelet packet decomposition on the frequency information to determine the optimal wavelet basis for wide - band measurement of the power system;

[0173] S3, perform multi - layer decomposition on the frequency information according to the optimal wavelet basis to obtain multiple wavelet frequency bands;

[0174] S4, perform component extraction on the wavelet frequency bands, and adaptively extract harmonic components of specific frequencies using the SVMD algorithm to obtain harmonic components corresponding to the wavelet frequency bands;

[0175] S5, perform denoising processing on the harmonic components using the wavelet coefficient correlation algorithm to obtain denoised harmonics:

[0176]

[0177] In the formula, c j (k) is the wavelet coefficient, j is the scale, k is the number of points, ∑ is the summation symbol, N is a natural number. By comparing the wavelet coefficient correlation with a set threshold, it is judged whether to retain the current wavelet coefficient, so as to achieve the denoising effect. Wavelet coefficient correlation denoising has better effects on non - stationary signals and signals containing non - local features. Compared with wavelet modulus maximum denoising, it can better retain the detailed features of the signal. For the case where the correlation between noise and signal is weak, it may cause signal distortion.

[0178] However, this algorithm is prone to loss of high - frequency parameters in application. To solve this problem, an improved EMD algorithm is introduced to perform preliminary decomposition on the signal and calculate its permutation entropy value. For the part with a higher noise - to - signal ratio, a correlation denoising method is used for processing. The denoised signal obtained after reconstruction better retains the wide - band information, thereby improving the accuracy and reliability of signal processing.

[0179] First, traditional EMD decomposition is a step - by - step decomposition process in the time domain for different fluctuation trends, which decomposes the overall non - stationary signal into stationary components, and these decomposition results are called intrinsic mode functions (IMFs). The detailed decomposition process is as follows: First, determine the local extrema of f(t) by interpolation to form the envelope line, take the average to obtain the overall trend component c(t) of the waveform, remove c(t) from the original signal and perform a convergence test, and continue to repeat the above decomposition process until the final signal residue r(t) remains after reaching the screening stop condition, completing the decomposition process.

[0180] Through this process, it can be clearly found that the advantage of EMD is that it is an adaptive signal decomposition method, which does not require prior assumptions about the signal model or basis functions. Therefore, it is applicable to various types of signals, especially non-linear and non-stationary signals. Moreover, EMD is a decomposition algorithm completely based on the time domain, making its decomposition results complete. All the output components can be directly restored without loss of the original signal by direct summation, as shown in Equation (31). Compared with the uncertainty of the inverse process of wavelet transform, it avoids the algorithm adaptation problem caused by the double decomposition process of the signal. Therefore, this processing process will not have any impact on the correlation noise reduction algorithm. In addition, this iterative decomposition and direct summation synthesis process is simple and easy to implement, without causing too much computational load and time pressure.

[0181]

[0182] In the formula, c i (t) is the overall trend component of the waveform, r K (t) is the signal residue, and k is 0, 1, 2...

[0183] Secondly, similar to most decomposition algorithms, performing Fourier spectrum analysis on the IMF reveals that the frequency distributions of its components are different. This decomposition method is mutually adaptable to the waveform data characteristics. Therefore, it does not have a strict frequency band division like wavelet decomposition. This adaptive decomposition process is one of the characteristics required in this paper. However, it is relatively difficult to determine the screening stop criterion for this adaptive decomposition process. Usually, it is judged based on a given fixed value, but this often loses part of the self-adaptability of the decomposition algorithm. As a result, during the decomposition process of different characteristic signals, it often causes under-decomposition or over-decomposition, both of which will introduce errors. Therefore, the adaptive stop mechanism is the key to ensuring the stability of each component.

[0184] In EMD, the number of screening iterations is directly determined by the screening stop criterion SSC. In the cases of "under-screening" and "over-screening", the number of iterations will affect the decomposition effect. If the number of iterations is insufficient and the decomposition is incomplete, decomposing multiple single-component signals into 1 IMF will result in too many uncorrelated components. While too many iterations will cause single-component signals to be decomposed into multiple IMFs, leading to waste of computing resources and reduction of decomposition accuracy, and will also affect the efficiency of the subsequent correlation noise reduction process. Therefore, it is crucial to reduce the computational load without sacrificing the decomposition accuracy to improve the noise reduction efficiency. Soft SSC can monitor the screening process of EMD separation. More importantly, it can select the optimal number of iterations. Starting from the principle of successive decomposition of EMD, the sum of the root mean square value and the absolute value of the kurtosis of the component extracted each time is always smaller than that of the component extracted in the previous time. Thus, a function related to the signal is defined:

[0185]

[0186] In the formula, n is the point value, k is the number of iterations, N k is the total number of iterations, and E k is the EMD decomposition component at the k-th iteration, and c k (n) is the n-point value of the current component at the k-th iteration. When α generated in this iteration k is relatively small compared to the previous iteration, it can be regarded as the stopping condition, and the screening process is automatically stopped according to the value of the objective function, so as to obtain better signal decomposition performance.

[0187] The idea of soft SSC realizes the self-adaptability of the screening process. Due to the fixed screening conditions of traditional EMD, the decomposition process is often easily affected by noise and error values. The adaptive stopping condition greatly improves the robustness of signal decomposition. Each screening will determine an optimal number of iterations, and no meaningless IMFs will be generated. The CPU time is proportional to both the number of iterations and the number of IMFs. It will not introduce errors due to under-screening, nor will it cause waste of computing resources and time due to over-screening.

[0188] To sum up, the decomposition process of adaptive stopping EMD can make up for the defects of correlation noise reduction and optimize its noise reduction effect. According to the principle of correlation noise reduction, signal noise reduction is carried out in each decomposition domain of the signal, and the noise of each segment can be deeply reduced on the premise of avoiding waveform distortion. In order to reduce unnecessary computational complexity and loss of signal feature information, permutation entropy detection is performed on the components after the signal is decomposed by adaptive EMD, and only the components dominated by noise are processed by correlation noise reduction. The specific algorithm flow chart is as Figure 16 shown.

[0189] First, the noisy signal needs to be decomposed into adaptive EMD components to obtain each component arranged from low frequency to high frequency. For the part with more Gaussian white noise, it needs to be removed by correlation noise reduction. Therefore, calculate the entropy value of each component, and arrange according to the entropy value. Set the components exceeding the noise limit value to be processed by correlation noise reduction. Finally, the denoised signal is obtained by adding each component. Since each component has been adaptively decomposed once, the number of wavelet decomposition layers can be reduced when the noisy signal is subjected to correlation processing, thus saving computational complexity and processing time, and better retaining the true information.

[0190] To verify the effect of the improved correlation noise reduction algorithm, a noise reduction test was carried out on the signal containing 10 dB Gaussian white noise. After experimental comparison, the db4 wavelet basis function was selected, and the correlation noise reduction processing of 3-layer wavelet decomposition was carried out on the high-noise signal. First, the adaptive preliminary decomposition of the signal is as Figure 17 shown.

[0191] In the original waveform decomposition diagram, the high-frequency components usually contain more noise. By calculating the permutation entropy values of each component, a comparison diagram of all components, such as the permutation entropy value curves of each component, can be obtained. Figure 18 As shown, it can be seen that for the IMF1-IMF2 components, wavelet correlation noise reduction processing is required.

[0192] After noise reduction by the improved correlation algorithm, the signal-to-noise ratio is 23.485 dB and the root mean square value is 0.43. Comparing with the record of the correlation noise reduction effect of the traditional wavelet noise reduction algorithm in the previous text, the improved noise reduction algorithm retains more effective information and ensures the noise reduction effect. After noise reduction processing, each component is reconstructed to obtain a comparison diagram of the improved correlation noise reduction effect between the low-noise signal and the original signal, as shown in Figure 19 As shown, the effect of the improved noise reduction algorithm is verified.

[0193] S6. Decompose the noise-reduced harmonic according to the EMD algorithm to obtain the signal residual that meets the preset stop condition;

[0194] S7. Calculate the measurement data corresponding to the waveforms of each harmonic according to the signal residual.

[0195] In the above technical solution, preferably, as shown in Figure 2 S2. Determine the optimal wavelet basis for wide-band measurement of the power system, including:

[0196] S21. Use the Mallat decomposition algorithm to determine the decomposition relationship between the standard orthogonal bases and the decomposition and synthesis relationship of the subspace coordinates corresponding to the decomposition relationship between the standard orthogonal bases;

[0197]

[0198] In equations (1)-(2), t is time, φ(t) is the wavelet function, ψ(t) is the wavelet function of the subspace coordinates, i is the number of nodes, j is the decomposition layer, k is the number of points in each layer, Z is an integer, m is an integer in Z, ∑ is the summation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter, and the coordinate sequences {h m ; m ∈ Z}, {g m ; m ∈ Z} are called the low-pass filter and the high-pass filter.

[0199] Thus, the decomposition relationship between its standard orthogonal bases is derived, and then the decomposition and synthesis relationship between the coordinates can be obtained, that is, the Mallat decomposition algorithm and the synthesis algorithm:

[0200] S22. Perform transform space decomposition on the decomposition relationship between the standard orthogonal bases according to the preset conditions to obtain the standard orthogonal bases of a closed subspace on the finite function space;

[0201]

[0202] (3)-(4) In the equations, j is the decomposition level, k is the number of points in each layer, Z is an integer, m is an integer in Z, and α jk is the coefficient sequence, and c j,k is the scaling coefficient, and d j,k is the wavelet coefficient, and h m-2k is the impulse response function of the low-pass filter, and g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

[0203] The Mallat algorithm can implement the decomposition operation of wavelet coefficients. This process is repeated continuously until the preset decomposition level is reached, thereby constructing a decomposition tree of wavelet transform, as shown in the wavelet decomposition schematic Figure 12 shown, where each node represents the approximate or detailed component of the signal at different scales. The algorithm has good stability and locality, can accurately capture the local features and time-frequency information of the signal, and is suitable for broadband and dynamic signal processing tasks.

[0204] Based on the idea of the above algorithm, the orthonormal basis and coordinate sequence of the low-frequency part of this level are always used to decompose the high-frequency and low-frequency information of the next level. The information shown by the decomposition tree obtained by repeating this process is that the low-frequency part is continuously refined and the high-frequency part is always maintained. If you want to obtain the high-frequency part information and improve the frequency resolution in the same way, you need to use the idea of orthogonal wavelet packets.

[0205] Compared with the traditional wavelet transform, the wavelet packet transform can more comprehensively describe the time-frequency characteristics of the signal. Its basic principle is to decompose and reconstruct the high-frequency and low-frequency signals simultaneously through a series of filter banks. When the preset decomposition level is reached or a certain stop condition is met, a wavelet packet decomposition tree can be constructed, as shown in the wavelet packet tree decomposition schematic Figure 13 shown.

[0206] The root node of the tree represents the original signal, and each node represents the output of a filter. You can select different decomposition levels or stop conditions according to needs to obtain the decomposition result suitable for specific applications.

[0207] S23. According to the properties of the filter itself: s 2i+1 (t) = ψ j，k (t), the orthonormal basis s2i(t) and s2i + 1(t) are deduced to form two mutually orthogonal orthonormal systems and the subspace relationship of the corresponding orthonormal systems after integer translation;

[0208]

[0209] Correspondingly, from a spatial perspective, in order to facilitate the analysis of the spatial decomposition relationship of wavelet packet transform, the energy finite function space L 2 A closed subspace on (R) is represented by a new symbol Sij if:

[0210]

[0211] The standard orthogonal basis of Sij can be deduced based on the properties of the high-pass filter and the low-pass filter: s2i(t) and s2i+1(t) can form two mutually orthogonal standard orthogonal systems after integer translation, and the corresponding subspace relationship is:

[0212]

[0213] In formulas (5)-(7), S is the root node original signal representing the output of a filter, S i j represents a closed subspace on L2(R), i is the number of nodes, j is the number of decomposition layers, k is the number of points in each layer, Z is an integer, and m is an integer in Z. ∑ is the accumulation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter.

[0214] S24, obtaining the decomposition result and the best wavelet basis of the corresponding wavelet packet according to the standard orthogonal system and its subspace relationship;

[0215] The decomposition result is expressed as:

[0216]

[0217] The optimal wavelet basis is expressed as:

[0218]

[0219] In formulas (8)-(9), α j，k is the coefficient sequence, d j,k is the wavelet coefficient, i is the number of nodes, j is the number of decomposition layers, k is the number of points in each layer, Z is an integer, and m is an integer in Z. ∑ is the accumulation symbol, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

[0220] The wavelet packet transform has higher resolution and richer frequency information, can capture the time-frequency characteristics of signals more accurately, and has better performance for non-stationary signals and signals with complex frequency structures.

[0221] From the perspective of qualitative analysis, a more applicable wavelet basis function was selected, and the application of wavelet analysis in signal analysis was implemented using MATLAB programming, verifying the feasibility of signal decomposition.

[0222] As Figure 3 shown, in the above technical solution, preferably, S3, perform multi-layer decomposition on the frequency information according to the optimal wavelet basis, including:

[0223] S31, perform wavelet transform on the optimal wavelet basis to obtain a single-period function representation of one period of the harmonic waveform:

[0224]

[0225] In the formula, f(t) is a periodic function, e ikt is the base, k is the number of points in each layer, and α k can be expressed as the coordinate sequence of the periodic function f(t) in the base e ikt .

[0226] S32, obtain the function representation of the entire harmonic waveform according to the single-period function representation:

[0227]

[0228] In the formula, G(t o , w) represents the short-time Fourier transform, t is time, w is the angular frequency, and f(t) is the waveform signal to be analyzed; g(t) is a window function that restricts the analysis time range, t0 is the specified time point indicating the magnitude of the component with frequency ω of the signal f(t) near the time point t0 and within the range determined by the function g(t), and e -jωt is the basis function.

[0229] S33, perform continuous wavelet transform on the function representation of the entire harmonic waveform according to the dynamic window characteristics to obtain the continuous wavelet dependent on the scale parameter and the function representation of the corresponding continuous wavelet;

[0230] The dynamic window characteristics are expressed as: ψ(t) ∈ L 2 (R)

[0231]

[0232] The continuous wavelet dependent on the scale parameter is expressed as:

[0233]

[0234] (12)-(14) In the formulas, C ψ is the wavelet packet energy, ψ(t) is the wavelet basis function, is the Fourier transform of the wavelet basis function, w is the angular frequency, R* represents all non-zero real numbers, R represents all real numbers. At the same time, the scale parameter a and the displacement parameter b are set. L is the entire signal space, e -jωt is the basis function. The function representation corresponding to the continuous wavelet is:

[0235]

[0236] The window of the wavelet transform is a rectangular window centered at the coordinates (b, ±ω0 / a), with a time-domain window width of aΔψ and a frequency-domain window width of Δψ / a.

[0237] In the formula, a is the scale parameter and b is the displacement parameter, is the complex conjugate of ψ(t).

[0238] S34, perform the inverse transform on the continuous wavelet and its corresponding function representation to obtain the inverse function representation corresponding to the continuous wavelet:

[0239]

[0240] In the formula is the inverse of the wavelet packet energy. The wavelet transform shows the change of the signal f(t) near the time point b as the scale parameter a continuously changes. S35, discretize the inverse function to obtain multiple wavelet frequency bands.

[0241] As Figure 4 shown, in the above technical solution, preferably, S35, the discretization process of the inverse function includes:

[0242] S351, define the preset function representation of the standard orthogonal basis that constitutes the corresponding wavelet frequency band:

[0243]

[0244] S352, obtain the linear representation of the corresponding wavelet frequency band and the discretized inverse function according to the preset function representation:

[0245]

[0246] Among them, the standard orthogonal basis that constitutes ψ(t) ∈ L 2 (R), then ψ(t) is an orthogonal wavelet;

[0247] S353, obtain the wavelet frequency band according to the coefficient sequence of the preset standard orthogonal basis:

[0248]

[0249] In equations (17)-(19), the coefficient sequence α j，k represents the continuous wavelet transform of f(t) when the scale parameter a = 2 -j , the displacement parameter b = 2 -j and k, where is the complex conjugate of ψ(t), ψ(t) is the wavelet basis function, i is the number of nodes, j is the decomposition level, k is the number of points per layer, and W f is the function of the continuous wavelet.

[0250] As Figure 5 shown, in the above technical solution, preferably, in S4, component extraction is performed on the wavelet frequency band, and the SVMD algorithm is used to adaptively extract harmonic components of a specific frequency, and the harmonic components corresponding to the wavelet frequency band include:

[0251] S41, the SVMD algorithm is used to perform successive modal decomposition on the wavelet frequency band to obtain each component signal u k (t), and the expression of the component signal u k (t) is:

[0252]

[0253] wherein, u k (t) is the k-th modal component, U k (t) is the amplitude function of the k-th component signal, φ k (t) is the phase function, then the frequency ω k (t) = dφ k (t) / dt, and all the above functions have non-negative characteristics;

[0254] Among them, the SVMD algorithm satisfies the following constrained variational problem:

[0255]

[0256] wherein, {ω k} is the column vector of the center frequencies of each component; δ(t) is the unit impulse signal, is the derivative operator with respect to the time variable t, j is the imaginary unit, f is the original input signal, that is, the target signal to be decomposed, and the equality condition ensures that the sum of all components is the original function;

[0257] For the solution of the constrained variational problem, an unconstrained variational problem with a penalty function is constructed by introducing the Lagrange multiplier λ and the penalty factor α:

[0258]

[0259] (22) is defined as the augmented Lagrangian function L({uk},{ω k}, λ), it is solved by an iterative method, and the iterative formula is expressed as:

[0260]

[0261] S42. Each time an iteration is completed, calculate whether the current modal component meets the following convergence conditions:

[0262]

[0263] In equations (22)-(24), f(t) is the original input signal, and λ k (t) is the Lagrange multiplier of the k-th modal constraint, is the updated result of the k-th modal component in the frequency domain in the (n + 1)-th iteration, is the Fourier transform of the original signal f(t), is the sum of the spectra obtained for the remaining modes except the k-th mode in the current iteration, is the representation of the Lagrange multiplier in the frequency domain, ω k is the center frequency of the k-th modal component, ω is the angular frequency, is the updated value of the Lagrange multiplier at the (n + 1)-th iteration, τ is the step size coefficient, K is the total number of modes to be decomposed, n is the number of iterations, k is the index of the modal component, is the updated value of the center frequency corresponding to the k-th modal component at the (n + 1)-th iteration, is the representation of the k-th mode in the frequency domain at the n-th iteration, is the representation of the k-th mode in the frequency domain at the (n + 1)-th iteration, ε is the convergence threshold, which is a small positive number. When holds, the updates of all modes are small enough, and the algorithm meets the convergence conditions and can stop the iteration;

[0264] S43. When the decomposed modal components meet the given accuracy, extract K components. The signal f(t) can be successively decomposed by the SVMD algorithm in the following way:

[0265]

[0266] In the formula, u K (t) is the K-th extracted modal component, x r (t) is the residual component, that is, the remaining part of the original signal after extracting the K-th component, x u (t) is the noise term or the residual term that has not been extracted yet;

[0267] S44. When extracting the K-th component u KAt time (t), construct the filter impulse response β K (t) and β n (t), to form a new constraint criterion for SVMD:

[0268]

[0269] In the formula, Q1 is the magnitude of the energy obtained after applying the filter β K (t) to the remaining signal. The smaller Q1 is, the less ineffective extraction of the remaining signal by the new filter, and it can be better distinguished from the remaining part. Q2 is the degree of overlap between the new component and the previously extracted components. The smaller Q2 is, the lower the degree of overlap or aliasing between the new component u K (t) and the previously extracted components, and the better the distinguishability between components. β K (t) is the filter impulse response to avoid the overlap of the frequency of the extracted component u K (t) and the residual component x r (t). β n (t) is the filter impulse response to effectively distinguish adjacent extracted components u K (t) and u K-1 (t);

[0270] The above filter is the key constraint condition to ensure the accurate extraction of each component by the SVMD algorithm. The expression of the filter is:

[0271]

[0272] Therefore, formula (21) is transformed into:

[0273]

[0274] In equations (27)-(28), is the weight function of the k-th modal component, which is used to adjust the weight of the signal in the frequency domain, so that the mode is mainly concentrated around its center frequency ω K Nearby, when ω is close to ω K , Obtains a larger value, indicating that the signal energy is mainly concentrated around ω K Nearby, when ω is far from ω K , Decays rapidly, so that the mode will not spread in the frequency range far from ω K ; Is the filtering weight function used to adjust the modal components in the frequency domain, which is used to control the frequency distribution of the extracted modal components to reduce the frequency aliasing between adjacent modes. When ω is close to ω n , Obtains a larger value, indicating that the mode is at ω nhas a larger weight nearby, maintaining the main frequency components. When ω is far from ω n , it rapidly decreases, so that this mode will not spread to other frequency regions, thereby reducing the spectral overlap between adjacent modes; ω K is the central frequency of the k-th mode, ω n is the central frequency of the n-th mode, and α is a penalty factor used to adjust the bandwidth effect of the filter;

[0275] S45. Use the NRBO optimization algorithm to optimize the penalty factor α, and take the minimum envelope entropy as the fitness function:

[0276]

[0277] where p i is the probability calculated for the envelope amplitude distribution of discrete components, S(i) is the envelope entropy of each component, N is a natural number, and i is 0, 1, 2...;

[0278] S46. Through the above SVMD and NRBO optimization steps, extract accurate harmonic components from the wavelet frequency band, and perform component extraction on the wavelet frequency band according to the harmonic components.

[0279] To verify the detection accuracy and efficiency of the improved algorithm for broadband signals, perform simulation analysis through the constructed signal model. First, set the waveform sampling frequency to 12.8 kHz and extract every 10 fundamental wave periods. Therefore, 2560 sampling points need to be processed. After wavelet packet decomposition, there are 8 frequency bands for harmonics within 64 times. The waveform diagram and spectrum diagram are as shown in the wavelet packet decomposition of the signal model Figure 14 .

[0280] The harmonics of the experimental signal model are evenly distributed in the 1st, 2nd, 4th, and 7th frequency bands, and there is a slight frequency band aliasing phenomenon between the frequency bands. Perform principal component extraction on the 4 frequency bands respectively, and all harmonic components can be obtained, as shown in Figure 15 (a), 15(b), 15(c), 15(d).

[0281] The optimization results of the penalty parameter α based on NRBO are: 6949, 12866, 18764, 20480 respectively. After decomposition, the harmonic waveforms of 10 fundamental wave periods can be obtained, but due to the endpoint effect at both ends, there is a distortion situation. Correct the extraction results by removing the first and last two fundamental wave periods, analyze and calculate the amplitude and phase angle of each waveform, and the test result statistical table is as shown in the table

[0282] Statistical Table of Test Results of Broadband Algorithm

[0283]

[0284]

[0285] The simulation results show that the improved broadband measurement algorithm can effectively extract each harmonic component. Through analysis and calculation, the maximum relative error of amplitude detection is 1.45%, and the maximum error of phase angle detection is about 0.1°, which verifies the rationality and accuracy of the algorithm in processing broadband and dynamic signals.

[0286] Through simulation verification, the algorithm demonstrates good dynamic extraction capabilities in broadband measurement. Based on the programming simulation of the algorithm, the processing of measured data is realized, verifying the actual processing ability of the algorithm and laying a solid foundation for the further application of subsequent detection algorithms. Subsequently, three algorithms for wavelet denoising technology were understood. After in-depth analysis and comparison, the correlation denoising algorithm was finally selected.

[0287] As Figure 6 shown, a power system broadband measurement system 100 based on the combination of wavelet packet and SVMD according to another embodiment of the present invention includes:

[0288] An acquisition module 10, configured to acquire the frequency information corresponding to the circuit signals at each measurement point in the circuit system;

[0289] A wavelet decomposition module 20, configured to perform wavelet packet decomposition on the frequency information to determine the optimal wavelet basis for power system broadband measurement;

[0290] A multi-layer decomposition module 30, configured to perform multi-layer decomposition on the frequency information according to the optimal wavelet basis to obtain multiple wavelet frequency bands;

[0291] A component extraction module 40, configured to perform modal decomposition using the SVMD algorithm on the basis of the wavelet frequency bands and sequentially extract each harmonic component, optimize the penalty parameter α of the SVMD algorithm, and minimize the envelope entropy as the objective function to reduce harmonic aliasing;

[0292] A denoising module 50, configured to perform denoising processing on the harmonic components using the wavelet coefficient correlation algorithm to obtain denoised harmonics;

[0293] A harmonic decomposition module 60, configured to decompose the denoised harmonics according to the EMD algorithm to obtain a signal residual that meets the preset stop condition;

[0294] A measurement module 70, configured to calculate measurement data corresponding to each harmonic waveform according to the signal residual.

[0295] As Figure 7 shown, in the above technical solution, preferably, the wavelet decomposition module 20 includes:

[0296] The orthogonal basis decomposition unit 21 is configured to determine the decomposition relationship between orthonormal bases and the decomposition and synthesis relationship of subspace coordinates corresponding to the decomposition relationship between orthonormal bases by using the Mallat decomposition algorithm;

[0297]

[0298] In equations (1)-(2), t is time, φ(t) is the wavelet function, ψ(t) is the subspace coordinate wavelet function, i is the number of nodes, j is the decomposition level, k is the number of points per layer, Z is an integer, m is an integer in Z, ∑ is the summation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter, and the coordinate sequences {h m ; m ∈ Z}, {g m ; m ∈ Z} are called the low-pass filter and the high-pass filter.

[0299] The space transformation unit 22 is configured to perform a transformed space decomposition on the decomposition relationship between orthonormal bases according to preset conditions to obtain an orthonormal basis of a closed subspace in the finite function space;

[0300]

[0301] In equations (3)-(4), j is the decomposition level, k is the number of points per layer, Z is an integer, m is an integer in Z, α jk is the coefficient sequence, c j,k is the scale coefficient, d j,k is the wavelet coefficient, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

[0302] The subspace decomposition unit 23 is configured to, according to the properties of the filter itself: s 2i+1 (t) = ψ j，k (t), deduce that the orthonormal bases s2i(t) and s2i + 1(t) form two mutually orthogonal orthonormal systems after integer translation and the subspace relationship of the corresponding orthonormal systems;

[0303]

[0304] In equations (5)-(7), S is the original signal of the root node representing the output of a filter, S ij represents a closed subspace on L2(R), where i is the number of nodes, j is the decomposition level, k is the number of points per layer, Z is an integer, and m is an integer in Z. ∑ is the summation symbol, and h m is the impulse response function of the low-pass filter, and g m is the impulse response function of the high-pass filter.

[0305] The optimal wavelet basis generation unit 24 is configured to obtain the decomposition result of the corresponding wavelet packet and the optimal wavelet basis according to the orthonormal system and its subspace relationship.

[0306] The decomposition result is expressed as:

[0307]

[0308] The optimal wavelet basis is expressed as:

[0309]

[0310] In equations (8)-(9), α j，k is the coefficient sequence, d j,k is the wavelet coefficient, i is the number of nodes, j is the decomposition level, k is the number of points per layer, Z is an integer, and m is an integer in Z. ∑ is the summation symbol, and h m-2k is the impulse response function of the low-pass filter, and g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

[0311] As Figure 8 shown, in the above technical solution, preferably, the multi-layer decomposition module 30 includes:

[0312] The wavelet transform unit 31 is configured to perform wavelet transform on the optimal wavelet basis to obtain the representation of a single-period function of one period of the harmonic waveform:

[0313]

[0314] In the formula, f(t) is a periodic function, e ikt is the basis, k is the number of points per layer, and α k can be expressed as the coordinate sequence of the periodic function f(t) in the basis e ikt In.

[0315] The waveform generation unit 32 is configured to obtain the function representation of the entire harmonic waveform according to the representation of the single-period function:

[0316]

[0317] In the formula, G(t o , w) represents the short-time Fourier transform, t is time, w is the angular frequency, and f(t) is the waveform signal to be analyzed; g(t) is the window function that restricts the analysis time range, t0 is the specified time point indicating the magnitude of the component with frequency ω whose signal f(t) is near the time point t0 and the range is determined by the function g(t), and e -jωt is the basis function.

[0318] The continuous wavelet generation unit 33 is configured to perform a continuous wavelet transform on the functional representation of the entire harmonic waveform according to the dynamic window characteristics, obtaining a continuous wavelet dependent on the scale parameter and the functional representation corresponding to the continuous wavelet;

[0319] The dynamic window characteristics are expressed as: ψ(t) ∈ L 2 (R)

[0320]

[0321] The continuous wavelet dependent on the scale parameter is expressed as:

[0322]

[0323] In equations (12)-(14), C ψ is the wavelet packet energy, ψ(t) is the wavelet basis function, is the Fourier transform of the wavelet basis function, w is the angular frequency, R* represents all non-zero real numbers, R represents all real numbers, and at the same time, the scale parameter a and the displacement parameter b are set, L is the entire signal space, and e -jωt is the basis function.

[0324] The functional representation corresponding to the continuous wavelet:

[0325]

[0326] The window of the wavelet transform is a rectangular window centered at the coordinates (b, ±ω0 / a), with a time-domain window width of aΔψ and a frequency-domain window width of Δψ / a.

[0327] In the formula, a is the scale parameter and b is the displacement parameter, is the complex conjugate of ψ(t).

[0328] The inverse transform unit 34 is configured to perform an inverse transform on the continuous wavelet and its corresponding functional representation, obtaining the inverse functional representation corresponding to the continuous wavelet:

[0329]

[0330] In the formula For the inverse of wavelet packet energy, the wavelet transform shows the variation of the signal f(t) near the time point b as the scale parameter a continuously changes. The wavelet transform shows the variation of the signal f(t) near the time point b as the scale parameter a continuously changes;

[0331] The discrete processing unit 35 is configured to discretize the inverse function to obtain multiple wavelet frequency bands.

[0332] As Figure 9 shown, in the above technical solution, preferably, the discrete processing unit 35 includes:

[0333] The wavelet frequency band generation sub-unit 351 is configured to define the preset function representation of the orthonormal basis constituting the corresponding wavelet frequency band:

[0334]

[0335] The frequency band linear function generation sub-unit 352 is configured to obtain the linear representation of the corresponding wavelet frequency band and the discretized inverse function according to the preset function representation:

[0336]

[0337] Among them, the orthonormal basis constituting ψ(t) ∈ L 2 (R), then ψ(t) is an orthogonal wavelet;

[0338] The orthogonal coefficient generation sub-unit 353 is configured to obtain the wavelet frequency band according to the coefficient sequence of the preset orthonormal basis:

[0339]

[0340] (17)-(19), the coefficient sequence α j,k represents the continuous wavelet transform of f(t) when the scale parameter a = 2 -j , the displacement parameter b = 2 -j k, is the complex conjugate of ψ(t), ψ(t) is the wavelet basis function, i is the number of nodes, j is the decomposition level, k is the number of points per layer, and W f is the function of the continuous wavelet.

[0341] As Figure 10 shown, in the above technical solution, preferably, the component extraction module 40 includes:

[0342] The SVMD principal component extraction unit 41 is configured to perform modal decomposition on the basis of the wavelet frequency band by using the SVMD algorithm and sequentially extract each harmonic component:

[0343] First, the wavelet frequency band is subjected to successive mode decomposition using the SVMD algorithm to obtain each component signal u k (t). The expression of the component signal u k (t) is as follows:

[0344]

[0345] In the formula, u k (t) is the k-th mode component, U k (t) is the amplitude function of the k-th component signal, φ k (t) is the phase function, then the frequency ω k (t) = dφ k (t) / dt, and all the above functions have non-negative characteristics;

[0346] Among them, the SVMD algorithm satisfies the following constrained variational problem:

[0347]

[0348] In the formula, {ω k} is the column vector of the center frequencies of each component; δ(t) is the unit impulse signal, is the derivative operator with respect to the time variable t, j is the imaginary unit, f is the original input signal, that is, the target signal to be decomposed, and the equality condition ensures that the sum of all components is the original function;

[0349] To solve the constrained variational problem, an unconstrained variational problem with a penalty function is constructed by introducing the Lagrange multiplier λ and the penalty factor α:

[0350]

[0351] (Equation (22)) is defined as the augmented Lagrangian function L({u k},{ω k},λ), and it is solved by an iterative method. The iterative formula is expressed as:

[0352]

[0353] Each time an iteration is completed, calculate whether the current mode component satisfies the following convergence condition:

[0354]

[0355] In equations (22)-(24), f(t) is the original input signal, λ k (t) is the Lagrange multiplier of the k-th mode constraint, is the updated result of the k-th mode component in the frequency domain in the (n + 1)-th iteration, is the Fourier transform of the original signal f(t), is the sum of the spectra of the remaining modes except the k-th mode obtained in the current iteration, is the representation of the Lagrange multiplier in the frequency domain, ω k is the center frequency of the k-th mode component, ω is the angular frequency, is the updated value of the Lagrange multiplier at the (n + 1)-th iteration, τ is the step size coefficient, K is the total number of modes to be decomposed, n is the number of iterations, k is the index of the mode component, is the updated value of the center frequency corresponding to the k-th mode component at the (n + 1)-th iteration, is the representation of the k-th mode in the frequency domain at the n-th iteration, is the representation of the k-th mode in the frequency domain at the (n + 1)-th iteration, ε is the convergence threshold, which is a small positive number. When the updates of all modes are small enough, the algorithm meets the convergence condition and can stop iterating;

[0356] When the decomposed mode components meet the given accuracy, extract K components. The signal f(t) can be successively decomposed by the SVMD algorithm in the following way:

[0357]

[0358] where u K (t) is the K-th extracted mode component, x r (t) is the residual component, that is, the part remaining of the original signal after extracting the K-th component, x u (t) is the noise term or the residual term not yet extracted;

[0359] When extracting the K-th component u K (t), construct the filter impulse responses β K (t) and β n (t) to form a new constraint criterion for SVMD:

[0360]

[0361] where Q1 is the magnitude of the energy obtained after applying the filter β K (t) to the remaining signal. The smaller Q1 is, the less ineffective extraction of the remaining signal by the new filter, and it can be better distinguished from the remaining part. Q2 is the degree of overlap between the new component and the previously extracted components. The smaller Q2 is, the lower the overlap or aliasing degree between the new component u K (t) and the previously extracted components, and the better the distinguishability between components. β K (t) is to avoid the frequency of the extracted component u K (t) from being the same as that of the residual component xr (t) Overlapping filter impulse responses, β n (t) To effectively distinguish adjacent extracted components u K (t) and u K-1 (t) of the filter impulse response;

[0362] The above filter is the key constraint condition to ensure the accurate extraction of each component by the SVMD algorithm. The expression of the filter is:

[0363]

[0364] Therefore, formula (21) is transformed into:

[0365]

[0366] In equations (27)-(28), is the weight function of the k-th modal component, which is used to adjust the weight of the signal in the frequency domain so that the mode is mainly concentrated around its center frequency ω K Nearby, when ω is close to ω K , obtains a large value, indicating that the signal energy is mainly concentrated around ω K Nearby, when ω is far from ω K , rapidly decays, so that the mode will not spread in the frequency range far from ω K ; is the filtering weight function used to adjust the modal components in the frequency domain, which is used to control the frequency distribution of the extracted modal components to reduce the frequency aliasing between adjacent modes. When ω is close to ω n , obtains a large value, indicating that the mode has a large weight near ω n , maintaining the main frequency components. When ω is far from ω n , rapidly decreases, so that the mode will not expand to other frequency regions, thereby reducing the spectral overlap between adjacent modes; ω K is the center frequency of the k-th mode, ω n is the center frequency of the n-th mode, and α is the penalty factor, which is used to adjust the influence of the filter bandwidth;

[0367] The NRBO optimization unit 42 is set to optimize the penalty parameter α of the SVMD algorithm and use the minimum envelope entropy as the objective function to reduce harmonic aliasing:

[0368]

[0369] In the formula, p iFor the probability of calculating the envelope amplitude distribution of discrete components, S(i) is the envelope entropy of each component, N is a natural number, and i is 0, 1, 2...;

[0370] Through the above SVMD and NRBO optimization steps, accurate harmonic components are extracted from the wavelet frequency band, and the wavelet frequency band is component-extracted according to the harmonic components.

[0371] Based on the above as Figure 5 and Figure 6 shown in the method, correspondingly, the embodiment of the present application also provides a computer-readable storage medium, on which a computer program is stored, and when the program is executed by a processor, the steps of the power system wide-band measurement method based on the combination of wavelet packet and SVMD in any of the above embodiments are implemented.

[0372] Based on such an understanding, the technical solution of the present application can be embodied in the form of a software product, and the software product can be stored in a non-volatile storage medium (which can be a CD-ROM, a USB flash drive, a mobile hard disk, etc.), including several instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods in various implementation scenarios of the present application.

[0373] Based on the above as Figure 5 and Figure 6 shown in the method, and Figure 7 shown in the virtual device embodiment, for the purpose of achieving the above object, the embodiment of the present application also provides a computer device, including a storage medium and a processor; the storage medium is used to store a computer program; the processor is used to execute the computer program to implement the steps of the power system wide-band measurement method based on the combination of wavelet packet and SVMD in any of the above embodiments.

[0374] Optionally, the computer device may further include a user interface, a network interface, a camera, a radio frequency (RF) circuit, sensors, an audio circuit, a WI-FI module, etc. The user interface may include a display screen (Display), an input unit such as a keyboard (Keyboard), etc., and the optional user interface may further include a USB interface, a card reader interface, etc. The network interface may optionally include a standard wired interface, a wireless interface (such as a Bluetooth interface, a WI-FI interface), etc.

[0375] Those skilled in the art can understand that the structure of a computer device provided in this embodiment does not limit the computer device, and it may include more or fewer components, or combine some components, or arrange different components.

[0376] The storage medium may further include an operating system and a network communication module. The operating system is a program for managing and storing the hardware and software resources of a computer device, and supports the operation of information processing programs and other software and / or programs. The network communication module is used to implement communication between components within the storage medium, as well as communication with other hardware and software in the entity device.

[0377] In the description of this specification, the descriptions of terms such as "one embodiment", "some embodiments", "specific embodiments", etc. mean that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or instance. Moreover, the specific features, structures, materials, or characteristics described can be combined in a suitable manner in any one or more embodiments or examples.

[0378] The foregoing is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention can have various changes and modifications. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A broadband measurement method for power system based on the combination of wavelet packet and SVMD, characterized in that: The following steps are involved: Step 1: Obtain frequency information corresponding to the circuit signal at each measurement point in the circuit system; Step 2: performing wavelet packet decomposition on the frequency information to determine the best wavelet basis for broadband measurement of the power system; Step 3: Perform multi-layer decomposition on the frequency information according to the optimal wavelet basis to obtain multiple wavelet frequency bands; Step 4: extracting components from the wavelet frequency band, and using the SVMD algorithm to adaptively extract harmonic components of specific frequencies to obtain harmonic components corresponding to the wavelet frequency band; Step 5: De-noising the harmonic components using a wavelet coefficient correlation algorithm to obtain de-noised harmonics; Step 6: Decompose the noise reduction harmonics according to the EMD algorithm to obtain a signal residual that meets a preset stop condition; Step 7: Calculate the measured data corresponding to each harmonic waveform based on the signal residual.

2. The method for broadband measurement of electric power system according to claim 1, characterized in that: Determine the optimal wavelet basis for broadband measurements in power systems, including: Using the Mallat decomposition algorithm, determining the decomposition relationship between standard orthogonal bases and the decomposition and synthesis relationship of the subspace coordinates corresponding to the decomposition relationship between the standard orthogonal bases; (1)-(2) Where t is time, φ(t) is the wavelet function, ψ(t) is the subspace coordinate wavelet function, i is the number of nodes, j is the number of decomposition levels, k is the number of points in each level, Z is an integer, and m is an integer in Z. ∑ is the accumulation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter, the coordinate sequence {h m ; m∈Z}, {g m ;m∈Z} is called a low-pass filter and a high-pass filter; Performing transformation space decomposition on the decomposition relationship between the standard orthogonal bases according to preset conditions to obtain a standard orthogonal base of a closed subspace on the finite function space; In formulas (3)-(4), j is the number of decomposition levels, k is the number of points in each level, Z is an integer, m is an integer in Z, and α jk is the coefficient sequence, c j,k is the scale factor, d j,k is the wavelet coefficient, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions; According to the nature of the filter itself: s 2i+1 (t) = ψ j，k (t), deriving that the standard orthogonal bases s2i(t) and s2i+1(t) are integer-shifted to form two mutually orthogonal standard orthogonal systems and the subspace relationship corresponding to the standard orthogonal systems; In formulas (5)-(7), S is the root node original signal representing the output of a filter, represents a closed subspace on L2(R), i is the number of nodes, j is the number of decomposition layers, k is the number of points in each layer, Z is an integer, and m is an integer in Z. ∑ is the accumulation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter; Finally, the decomposition result of the corresponding wavelet packet and the optimal wavelet basis are obtained according to the standard orthogonal system and its subspace relationship; The decomposition result is expressed as: The optimal wavelet basis is expressed as: In formulas (8)-(9), α j，k is the coefficient sequence, d j,k is the wavelet coefficient, i is the number of nodes, j is the number of decomposition layers, k is the number of points in each layer, Z is an integer, and m is an integer in Z. ∑ is the accumulation symbol, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

3. The method for broadband measurement of electric power system according to claim 2, characterized in that: The optimal wavelet basis performs multi-layer decomposition on the frequency information, including: The optimal wavelet basis is subjected to wavelet transform to obtain a single periodic function representation of one period of the harmonic waveform: Where f(t) is a periodic function, e ikt is the basis, k is the number of points in each layer, α k It can be expressed as the periodic function f(t) in the basis e ikt The coordinate sequence in ; The function representation of the entire harmonic waveform is obtained according to the single period function representation: In the formula, G(t o , w) represents short-time Fourier transform, t is time, w is angular frequency, f(t) is the waveform signal to be analyzed; g(t) is a window function that limits the time range of analysis, t0 is the specified time point, indicating that the signal f(t) is near the time point t0 and the range is determined by the function g(t) The size of the component with a frequency of ω, e -jωt is the basis function; Performing a continuous wavelet transform on the function representation of the entire harmonic waveform according to the dynamic window characteristics to obtain a continuous wavelet that depends on a scale parameter and a function representation corresponding to the continuous wavelet; The dynamic window characteristic is expressed as: ψ(t)∈L 2 (R) The continuous wavelet representation that depends on the scale parameter is: In formulas (12)-(14), C ψ is the wavelet packet energy, ψ(t) is the wavelet basis function, is the Fourier transform of the wavelet basis function, w is the angular frequency, R* represents all non-zero real numbers, R represents all real numbers, and at the same time, the scale parameter a and the displacement parameter b are set, L is the entire signal space, e -jωt is the basis function; The function corresponding to the continuous wavelet is expressed as: The window of wavelet transform is a rectangular window centered at coordinate (b, ±ω0 / a), with a time domain window width of aΔψ and a frequency domain window width of Δψ / a; Where a is the scale parameter and b is the displacement parameter. is the complex conjugate of ψ(t); Perform an inverse transformation on the continuous wavelet and its corresponding function representation to obtain the inverse function representation corresponding to the continuous wavelet: In the formula To invert the wavelet packet energy, the wavelet transform shows the changes of the signal f(t) around the time point b as the scale parameter a changes continuously; The inverse function is discretized to obtain a plurality of wavelet frequency bands.

4. The method for broadband measurement of electric power system according to claim 3, characterized in that: Discretizing the inverse function includes: The preset function that constitutes the standard orthogonal basis corresponding to the wavelet frequency band is defined as: According to the preset function representation, a linear representation of the inverse function corresponding to the wavelet frequency band and the discretization processing is obtained: Among them, the composition of ψ(t)∈L 2 (R), then ψ(t) is an orthogonal wavelet; According to the preset coefficient sequence of the standard orthogonal basis, the wavelet frequency band is obtained: In formulas (17)-(19), the coefficient sequence α j，k Indicates that when the scale parameter a=2 -j , displacement parameter b = 2 -j The continuous wavelet transform of f(t) at time k, is the complex conjugate of ψ(t), ψ(t) is the wavelet basis function, i is the number of nodes, j is the number of decomposition layers, k is the number of points in each layer, W f is a function of continuous wavelet.

5. The method for broadband measurement of electric power system according to claim 4, characterized in that: Extracting components from the wavelet frequency band includes: First, the SVMD algorithm is used to perform successive modal decomposition on the wavelet frequency band to obtain each component signal u k (t), the component signal u k The expression of (t) is: In the formula, u k (t) is the kth modal component, U k (t) is the amplitude function of the kth component signal, is a phase function, then the frequency And the above functions all have non-negative characteristics; The SVMD algorithm satisfies the following constrained variational problem: In the formula, {ω k } is the column vector of the center frequency of each component; δ(t) is the unit pulse signal, is the derivative operator with respect to the time variable t, j is the imaginary unit, f is the original input signal, i.e. the target signal to be decomposed, and the equality condition ensures that the sum of all components is the original function; The solution to the constrained variational problem is to construct an unconstrained variational problem with a penalty function by introducing the Lagrange multiplier λ and the penalty factor α: (22) is defined as the augmented Lagrangian function L({u k },{ω k },λ), and solve it by iteration. The iteration formula is expressed as: Each time an iteration is completed, the current modal component is calculated to see if it meets the following convergence conditions: In formulas (22)-(24), f(t) is the original input signal, λ k (t) is the Lagrange multiplier of the kth modal constraint, is the updated result of the kth modal component in the frequency domain in the n+1th iteration, is the Fourier transform of the original signal f(t), is the sum of the spectra of the modes except the kth mode obtained in the current iteration, is the representation of Lagrange multipliers in the frequency domain, ω k is the center frequency of the kth modal component, ω is the angular frequency, is the updated value of the Lagrange multiplier at the n+1th iteration, τ is the step coefficient, K is the total number of modes to be decomposed, n is the number of iterations, k is the index of the modal component, is the updated value of the center frequency corresponding to the kth modal component at the n+1th iteration, is the representation of the kth mode in the frequency domain at the nth iteration, is the representation of the kth mode in the frequency domain at the n+1th iteration, ε is the convergence threshold, which is a small positive number. When , the updates of all modes are small enough, the algorithm meets the convergence condition and the iteration can be stopped; When the decomposed modal components meet the given accuracy, K components are extracted, and the signal f(t) can be decomposed successively in the following way through the SVMD algorithm: In the formula, u K (t) is the Kth extracted modal component, x r (t) is the residual component, that is, the remaining part of the original signal after extracting the Kth component, x u (t) is the noise term or the residual term that has not been extracted; When extracting the Kth component u K (t), construct the filter impulse response β K (t) and β n (t), forming a new constraint criterion for SVMD: Where Q1 is the residual signal when filter β is applied to it. K (t), the smaller Q1 is, the less invalid the new filter extracts from the remaining signal, and the better it can be distinguished from the remaining part. Q2 is the overlap between the new component and the previously extracted component. The smaller Q2 is, the new component u K The lower the overlap or aliasing between (t) and previously extracted components, the better the distinction between components, β K (t) To avoid extracting component u K (t) frequency and residual component x r (t) Overlapped filter impulse response, β n (t) is to effectively distinguish the adjacent extracted components u K (t) and u K-1 (t) filter impulse response; The above filter is the key constraint to ensure that the SVMD algorithm can accurately extract each component. The expression of the filter is: Therefore, formula (21) is transformed into: In formulas (27)-(28), is the weight function of the kth modal component, which is used to adjust the weight of the signal in the frequency domain so that the mode is mainly concentrated at its center frequency ω K Nearby, when ω is close to ω K hour, A larger value indicates that the signal energy is mainly concentrated in ω K Nearby, when ω is far away from ω K hour, Rapidly decays, so that the mode will not move away from ω K Diffusion within the frequency range; is the filter weight function used to adjust the modal component in the frequency domain, which is used to control the frequency distribution of the extracted modal component to reduce the frequency aliasing between adjacent modes. n hour, A larger value indicates that the mode is in ω n It has a larger weight near ω, keeping the main frequency components. n hour, Decreases rapidly, so that the mode will not extend to other frequency regions, thereby reducing the spectrum overlap between adjacent modes; ω K is the center frequency of the kth mode, ω n is the center frequency of the nth mode, and α is the penalty factor used to adjust the bandwidth effect of the filter; Then, the penalty factor α is optimized using the NRBO optimization algorithm, and the minimum envelope entropy is used as the fitness function: In the formula, p i is the probability of calculating the envelope amplitude distribution of discrete components, S(i) is the envelope entropy of each component, N is a natural number, i is 0, 1, 2...; Through the above SVMD and NRBO optimization steps, accurate harmonic components are extracted from the wavelet frequency band, and components of the wavelet frequency band are extracted according to the harmonic components.

6. A power system broadband measurement system based on the combination of wavelet packets and SVMD, characterized in that: include: An acquisition module is configured to acquire frequency information corresponding to circuit signals at various measurement points in the circuit system; A wavelet decomposition module is configured to perform wavelet packet decomposition on the frequency information to determine an optimal wavelet basis for broadband measurement of the power system; A multi-layer decomposition module is configured to perform multi-layer decomposition on the frequency information according to the optimal wavelet basis to obtain a plurality of wavelet frequency bands; The SVMD principal component extraction module is configured to perform modal decomposition based on the wavelet frequency band using the SVMD algorithm and extract each harmonic component in turn; An NRBO optimization module is configured to optimize a penalty parameter α of the SVMD algorithm and to reduce harmonic aliasing by minimizing envelope entropy as an objective function; A noise reduction module is configured to perform noise reduction processing on the harmonic components using a wavelet coefficient correlation algorithm to obtain noise-reduced harmonics; A harmonic decomposition module is configured to decompose the noise reduction harmonics according to an EMD algorithm to obtain a signal residual that meets a preset stop condition; The calculation module is configured to calculate the measurement data corresponding to each harmonic waveform according to the signal residual.

7. The power system broadband measurement system according to claim 6, characterized in that: The wavelet decomposition module includes: An orthogonal basis decomposition unit is configured to determine a decomposition relationship between standard orthogonal bases and a decomposition and synthesis relationship of subspace coordinates corresponding to the decomposition relationship between the standard orthogonal bases using a Mallat decomposition algorithm; (1)-(2) Where t is time, φ(t) is the wavelet function, ψ(t) is the subspace coordinate wavelet function, i is the number of nodes, j is the number of decomposition levels, k is the number of points in each level, Z is an integer, and m is an integer in Z. ∑ is the accumulation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter, the coordinate sequence {h m ; m∈Z}, {g m ;m∈Z} is called a low-pass filter and a high-pass filter; A space transformation unit is configured to perform a transformation space decomposition on the decomposition relationship between the standard orthogonal bases according to a preset condition to obtain a standard orthogonal base of a closed subspace on a finite function space; In formulas (3)-(4), j is the number of decomposition levels, k is the number of points in each level, Z is an integer, m is an integer in Z, and α jk is the coefficient sequence, c j,k is the scale factor, d j,k is the wavelet coefficient, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions; The subspace decomposition unit is set to be used according to the properties of the filter itself: s 2i+1 (t) = ψ j，k (t), deriving that the standard orthogonal bases s2i(t) and s2i+1(t) are integer-shifted to form two mutually orthogonal standard orthogonal systems and the subspace relationship corresponding to the standard orthogonal systems; In formulas (5)-(7), S is the root node original signal representing the output of a filter, represents a closed subspace on L2(R), i is the number of nodes, j is the number of decomposition layers, k is the number of points in each layer, Z is an integer, and m is an integer in Z. ∑ is the accumulation symbol, h m is the impulse response function of the low-pass filter, g m is the impulse response function of the high-pass filter; An optimal wavelet basis generating unit is configured to obtain a decomposition result and an optimal wavelet basis of a corresponding wavelet packet according to the standard orthogonal system and its subspace relationship; The decomposition result is expressed as: The optimal wavelet basis is expressed as: In formulas (8)-(9), α j，k is the coefficient sequence, d j,k is the wavelet coefficient, i is the number of nodes, j is the number of decomposition layers, k is the number of points in each layer, Z is an integer, and m is an integer in Z. ∑ is the accumulation symbol, h m-2k is the impulse response function of the low-pass filter, g m-2k is the impulse response function of the high-pass filter, and are the conjugate functions of their corresponding impulse response functions respectively.

8. The power system broadband measurement system according to claim 6, characterized in that: The multi-layer decomposition module includes: The wavelet transform unit is configured to perform wavelet transform on the optimal wavelet basis to obtain a single periodic function representation of one period of the harmonic waveform: Where f(t) is a periodic function, e ikt is the basis, k is the number of points in each layer, α k It can be expressed as the periodic function f(t) in the basis e ikt The coordinate sequence in ; A waveform generating unit is configured to obtain a function representation of the entire harmonic waveform according to the single period function representation: In the formula, G(t o , w) represents short-time Fourier transform, t is time, w is angular frequency, f(t) is the waveform signal to be analyzed; g(t) is a window function that limits the time range of analysis, t0 is the specified time point, indicating that the signal f(t) is near the time point t0 and the range is determined by the function g(t) The size of the component with a frequency of ω, e -jωt is the basis function; A continuous wavelet generation unit is configured to perform a continuous wavelet transform on the function representation of the entire harmonic waveform according to a dynamic window characteristic to obtain a scale-parameter-dependent continuous wavelet and a function representation corresponding to the continuous wavelet; The dynamic window characteristic is expressed as: ψ(t)∈L 2 (R) The continuous wavelet representation that depends on the scale parameter is: In formulas (12)-(14), C ψ is the wavelet packet energy, ψ(t) is the wavelet basis function, is the Fourier transform of the wavelet basis function, w is the angular frequency, R* represents all non-zero real numbers, R represents all real numbers, and at the same time, the scale parameter a and the displacement parameter b are set, L is the entire signal space, e -jωt is the basis function; the function corresponding to the continuous wavelet is expressed as: The window of wavelet transform is a rectangular window centered at coordinate (b, ±ω0 / a), with a time domain window width of aΔψ and a frequency domain window width of Δψ / a; Where a is the scale parameter and b is the displacement parameter. is the complex conjugate of ψ(t); The inverse transform unit is configured to perform an inverse transform on the continuous wavelet and its corresponding function representation to obtain an inverse function representation corresponding to the continuous wavelet: In the formula To invert the wavelet packet energy, the wavelet transform represents the change of the signal f(t) around the time point b as the scale parameter a changes continuously; the discrete processing unit is configured to discretize the inverse function to obtain multiple wavelet frequency bands.

9. The power system broadband measurement system according to claim 8, characterized in that: The discrete processing unit comprises: The wavelet frequency band generating subunit is configured to define a preset function representation of the standard orthogonal basis constituting the corresponding wavelet frequency band: The frequency band linear function generating subunit is configured to obtain a linear representation of the inverse function corresponding to the wavelet frequency band and the discretization processing according to the preset function representation: Among them, the composition of ψ(t)∈L 2 (R), then ψ(t) is an orthogonal wavelet; The orthogonal coefficient generating subunit is configured to obtain the wavelet frequency band according to a preset coefficient sequence of the standard orthogonal basis: In formulas (17)-(19), the coefficient sequence α j，k It means that when the scale parameter α=2 -j , displacement parameter b = 2 -j The continuous wavelet transform of f(t) at time k, is the complex conjugate of ψ(t), ψ(t) is the wavelet basis function, i is the number of nodes, j is the number of decomposition layers, k is the number of points in each layer, W f is a function of continuous wavelet.

10. The power system broadband measurement system according to claim 6, characterized in that: The SVMD principal component extraction module includes: Variational constraint unit: used to construct the constrained variational problem of SVMD for the input wavelet decomposition band signal; Alternating multiplier iteration unit: used to iteratively solve each modal component based on the introduced Lagrange multiplier λ and penalty factor α using the alternating direction multiplier algorithm; Parameter optimization unit: used to use the NRBO optimization algorithm to adaptively optimize the penalty factor α, and take the minimum envelope entropy as the objective function; Harmonic extraction unit: used for extracting the harmonic components corresponding to the frequency band according to each modal component obtained by SVMD decomposition, and outputting them to the noise reduction module.

Images