Power system broadband measurement method and system based on wavelet packet and SVMD

Abstract

The present invention relates to a method for measuring broadband power system based on the combination of wavelet packets and SVMD, comprising the following steps: obtaining frequency information corresponding to circuit signals at various measurement points in the circuit system; performing wavelet packet decomposition on the frequency information to determine the optimal wavelet basis for broadband power system measurement; performing multi-layer decomposition on the frequency information based on the optimal wavelet basis to obtain multiple wavelet frequency bands; performing component extraction on the wavelet frequency bands to obtain harmonic components corresponding to the wavelet frequency bands; performing denoising on the harmonic components using a wavelet coefficient correlation algorithm to obtain de-noised harmonics; decomposing the de-noised harmonics according to an EMD algorithm to obtain signal residuals that meet a preset stop condition; and obtaining measurement data corresponding to each harmonic waveform based on the signal residuals. In the technical solution of the present invention, the time-frequency characteristics of the signal can be captured more accurately, and better performance is achieved for non-stationary signals and signals with complex frequency structures.

Description

Technical Field

[0001] The present invention relates to the technical field of electronic measurement of power systems, and in particular to a method and system for broadband measurement of power systems based on a combination of wavelet packets and SVMD. Background Art

[0002] For an ideal three-phase power system, the basic parameters that measurement devices need to measure are the amplitude and phase angle of the power frequency voltage. However, the voltage and current waveforms obtained at various measurement points in the power grid are typically not three-phase symmetrical sinusoidal waves. This is related to the connected equipment. If the power supply is connected by power electronic devices and the load is nonlinear, a large number of harmonics will be injected into the system, causing distortion of the voltage and current waveforms. With the widespread application of nonlinear loads and power electronic devices, monitoring high-order harmonic information in broadband signals has become an important research direction in new power systems. Compared with traditional wavelet transforms, wavelet packet transforms can more comprehensively describe the time-frequency characteristics of signals. However, while wavelet decomposition algorithms can effectively divide the frequency band of broadband signals, they cannot extract and detect harmonics of specific frequencies, which affects the accuracy and reliability of measurement results.

[0003] The present invention introduces the SVMD (Successive Variational Mode Decomposition, SVMD) algorithm to further process the wavelet decomposition band, and combines 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 effective measurement of broadband 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 art.

[0005] To this end, the purpose of the present invention is to provide a power system broadband measurement method and system based on the combination of wavelet packets 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 objectives, the technical solution of the first aspect of the present invention provides a power system broadband measurement method based on a combination of wavelet packets and SVMD, comprising the following steps:

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

[0008] Step 2: performing wavelet packet decomposition on the frequency information to determine the optimal wavelet basis for broadband measurement of the power system;

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

[0010] 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;

[0011] Step 5: De-noising the harmonic components using a wavelet coefficient correlation algorithm to obtain de-noised harmonics;

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

[0013] Step 7: Calculate the measured data corresponding to each harmonic waveform based on the signal residual.

[0014] The component extraction of the wavelet frequency band in the above step 4 includes:

[0015] First, the SVMD algorithm is used to perform successive modal decomposition on the wavelet frequency band to obtain the component signals u k (t), the component signal u k The expression of (t) is:

[0016]

[0017] Where 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;

[0018] The SVMD algorithm satisfies the following constrained variational problem:

[0019]

[0020] 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, that is, the target signal to be decomposed, and the equality condition ensures that the sum of all components is the original function;

[0021] 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 α:

[0022]

[0023] (22) is defined as the augmented Lagrangian function L({u k},{ω k},λ), and solve it by iteration. The iteration formula is expressed as:

[0024]

[0025] Each time an iteration is completed, the current modal component is calculated to see if it meets the following convergence conditions:

[0026]

[0027] 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 k-th mode obtained in the current iteration, is the representation of Lagrange multiplier in 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;

[0028] When the decomposed modal components meet the given accuracy, K components are extracted and the signal f(t) can be decomposed successively by the SVMD algorithm in the following way:

[0029]

[0030] Where 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;

[0031] 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:

[0032]

[0033] Where Q1 is the value of the filter β when the filter β is applied to the residual signal. 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 degree of 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, and the better the 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 adjacent extracted components u K (t) and u K-1 (t) filter impulse response;

[0034] The above filter is the key constraint to ensure that the SVMD algorithm can accurately extract each component. The expression of the filter is:

[0035]

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

[0037]

[0038] 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 K Diffusion within the frequency range; is a filter 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. n hour, A larger value indicates that the mode is in ω nIt has a larger weight near ω, keeping the main frequency components. n hour, Decreases rapidly so that the mode does 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;

[0039] Then the penalty factor α is optimized using the NRBO optimization algorithm, and the minimum envelope entropy is used as the fitness function:

[0040]

[0041] 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, i is 0, 1, 2...;

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

[0043] In the above technical solution, preferably, determining the optimal wavelet basis for broadband measurement of the power system includes:

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

[0045]

[0046] (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 layers, k is the number of points in each layer, Z is an integer, m is an integer in Z, ∑ is the accumulation symbol, and 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} are called low-pass filters and high-pass filters.

[0047] According to the preset conditions, the decomposition relationship between the standard orthogonal bases is decomposed into a transformation space to obtain a standard orthogonal base of a closed subspace on the finite function space;

[0048]

[0049] In formulas (3)-(4), j is the number of decomposition layers, 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, 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.

[0050] According to the nature of the filter itself: s 2i+1 (t)=ψ j，k (t), it is deduced that the standard orthogonal bases s2i(t) and s2i+1(t) take integer translations to form two mutually orthogonal standard orthogonal systems and the subspace relationship of the corresponding standard orthogonal systems;

[0051]

[0052] 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, 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.

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

[0054] The decomposition result is expressed as:

[0055]

[0056] The optimal wavelet basis is expressed as:

[0057]

[0058] 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, 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.

[0059] In the above technical solution, preferably, performing multi-layer decomposition on the frequency information according to the optimal wavelet basis includes:

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

[0061]

[0062] 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 .

[0063] The function representation of the entire harmonic waveform is obtained based on the single period function representation:

[0064]

[0065] 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 analysis time range, 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.

[0066] According to the dynamic window characteristics, the function representation of the entire harmonic waveform is transformed by continuous wavelet transform to obtain the continuous wavelet dependent on the scale parameter and the function representation of the corresponding continuous wavelet;

[0067] The dynamic window characteristic is expressed as: ψ(t)∈L 2 (R)

[0068]

[0069] The continuous wavelet expression that depends on the scale parameter is:

[0070]

[0071] 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, corresponding to the continuous wavelet function representation:

[0072]

[0073] The wavelet transform window is centered at the coordinate (b, ±ω0 / a), with a time domain window width of aΔψ and a frequency domain window width of Δψ / a. In the formula, a is the scale parameter and b is the displacement parameter. is the complex conjugate of ψ(t).

[0074] Perform inverse transformation on the continuous wavelet and its corresponding function representation to obtain the inverse function representation of the corresponding continuous wavelet:

[0075]

[0076] In the formula To invert the wavelet packet energy, the wavelet transform shows the change of the signal f(t) around the time point b as the scale parameter a changes continuously; the inverse function is discretized to obtain multiple wavelet frequency bands.

[0077] In the above technical solution, preferably, discretizing the inverse function includes:

[0078] Define the preset function representation that constitutes the standard orthogonal basis of the corresponding wavelet frequency band:

[0079]

[0080] According to the preset function representation, the linear representation of the corresponding wavelet frequency band and the inverse function after discretization is obtained:

[0081]

[0082] Among them, ψ(t)∈L 2 (R), then ψ(t) is an orthogonal wavelet;

[0083] According to the coefficient sequence of the preset standard orthogonal basis, the wavelet frequency band is obtained:

[0084]

[0085] 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 continuous wavelet function.

[0086] The technical solution of the second aspect of the present invention provides a power system broadband measurement system based on the combination of wavelet packets and SVMD, comprising:

[0087] an acquisition module configured to acquire frequency information corresponding to circuit signals at various measurement points in the circuit system;

[0088] 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;

[0089] a multi-layer decomposition module configured to perform multi-layer decomposition on the frequency information according to the optimal wavelet basis to obtain a plurality of wavelet frequency bands;

[0090] An SVMD principal component extraction module is configured to perform modal decomposition based on the wavelet frequency band using an SVMD algorithm and extract each harmonic component in sequence;

[0091] An NRBO optimization module is configured to optimize a penalty parameter α of the SVMD algorithm and use envelope entropy minimization as an objective function to reduce harmonic aliasing;

[0092] a noise reduction module configured to perform noise reduction processing on the harmonic components using a wavelet coefficient correlation algorithm to obtain noise-reduced harmonics;

[0093] a harmonic decomposition module configured to decompose the noise reduction harmonics according to an EMD algorithm to obtain a signal residual that meets a preset stop condition;

[0094] The calculation module is configured to calculate the measurement data corresponding to each harmonic waveform based on the signal residual.

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

[0096] 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;

[0097]

[0098] (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 layers, k is the number of points in each layer, Z is an integer, m is an integer in Z, ∑ is the accumulation symbol, and 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}, {gm ; m∈Z} is called a low-pass filter and a high-pass filter

[0099] 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 the finite function space;

[0100]

[0101] In formulas (3)-(4), j is the number of decomposition layers, 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, 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.

[0102] The subspace decomposition unit is set to be used according to the properties of the filter itself: s 2i+1 (t)=ψ j，k (t), it is deduced that the standard orthogonal bases s2i(t) and s2i+1(t) take integer translations to form two mutually orthogonal standard orthogonal systems and the subspace relationship of the corresponding standard orthogonal systems;

[0103]

[0104] 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, 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.

[0105] The optimal wavelet basis generating unit is configured to obtain the decomposition result and the optimal wavelet basis of the corresponding wavelet packet according to the standard orthogonal system and its subspace relationship;

[0106] The decomposition result is expressed as:

[0107]

[0108] The optimal wavelet basis is expressed as:

[0109]

[0110] 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, 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.

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

[0112] 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:

[0113]

[0114] 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 .

[0115] The waveform generation unit is configured to obtain a functional representation of the entire harmonic waveform according to a single periodic functional representation:

[0116]

[0117] 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 analysis time range, 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.

[0118] a continuous wavelet generation unit 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;

[0119] The dynamic window characteristic is expressed as: ψ(t)∈L 2 (R)

[0120]

[0121] Scale-parameter-dependent continuous wavelet representation:

[0122]

[0123] 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:

[0124]

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

[0126] Where a is the scale parameter and b is the displacement parameter. is the complex conjugate of ψ(t).

[0127] The inverse transform unit is configured to perform an inverse transform on the continuous wavelet and its corresponding function representation to obtain the inverse function representation of the corresponding continuous wavelet:

[0128]

[0129] 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.

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

[0131] The wavelet band generation subunit is configured to define a preset function representation of a standard orthogonal basis constituting a corresponding wavelet band:

[0132]

[0133] The frequency band linear function generating subunit is configured to obtain a linear representation of the corresponding wavelet frequency band and the discretized inverse function according to a preset function representation:

[0134]

[0135] Among them, ψ(t)∈L 2 (R), then ψ(t) is an orthogonal wavelet;

[0136] The orthogonal coefficient generating subunit is configured to obtain a wavelet frequency band according to a coefficient sequence of a preset standard orthogonal basis:

[0137]

[0138] 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 continuous wavelet function.

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

[0140] Variational constraint unit: used to construct the constrained variational problem of SVMD based on the input wavelet decomposition band signal;

[0141] 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;

[0142] Parameter optimization unit: used to use the NRBO optimization algorithm to adaptively optimize the penalty factor α, and use the minimum envelope entropy as the objective function;

[0143] The harmonic extraction unit is used to extract the harmonic components corresponding to the frequency band according to the modal components obtained by SVMD decomposition, and output them to the noise reduction module.

[0144] Compared with existing technologies, the advantages of the power system broadband measurement method and system based on the combination of wavelet packets and SVMD provided by the present invention are: the wavelet packet transform has higher resolution and richer frequency information, can more accurately capture the time-frequency characteristics of the signal, and has better performance for non-stationary signals and signals with complex frequency structures. It also has greater flexibility and adjustability, allowing for deeper or coarser decomposition of the signal as needed, thereby obtaining time-frequency representations at different scales.

[0145] The basic wavelet packet decomposition algorithm can only complete the frequency band division and cannot decompose the subcomponents in detail. By combining the frequency band decomposition function of the wavelet packet algorithm and the adaptive extraction function of SVMD, the problem of selecting the wavelet basis and the number of decomposition layers can be effectively solved.

[0146] The proposed method, based on combining these two algorithms, not only reduces the number of wavelet decomposition layers but also significantly improves the detection efficiency of broadband signals, adaptively extracting the principal components of each frequency band. The improved broadband measurement algorithm effectively extracts each subharmonic component. Analysis and calculations show that the maximum relative error in amplitude detection is 1.45%, and the maximum error in phase angle detection is approximately 0.1°, validating the algorithm's rationality and accuracy in processing broadband, dynamic signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0147] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the following description of the embodiments with reference to the accompanying drawings, in which:

[0148] Figure 1 A flowchart of a method according to an embodiment of the present invention is shown;

[0149] Figure 2 shows a flowchart of step S2 involved in an embodiment of the present invention;

[0150] Figure 3 shows a flowchart of step S3 involved in an embodiment of the present invention;

[0151] Figure 4 shows a flowchart of step S35 involved in an embodiment of the present invention;

[0152] Figure 5 shows a flowchart of step S4 involved in an embodiment of the present invention;

[0153] Figure 6 shows a structural block diagram of a system involved in an embodiment of the present invention;

[0154] Figure 7 shows a structural block diagram of a wavelet decomposition module involved in an embodiment of the present invention;

[0155] Figure 8 It shows a structural block diagram of a multi-layer decomposition module involved in an embodiment of the present invention;

[0156] Figure 9 shows a structural block diagram of a discrete processing unit involved in an embodiment of the present invention;

[0157] Figure 10 shows a structural block diagram of a component extraction module involved in an embodiment of the present invention;

[0158] Figure 11 A schematic diagram of the decomposition of a single operation involved in an embodiment of the present invention is shown;

[0159] Figure 12 A schematic diagram of wavelet decomposition according to an embodiment of the present invention is shown;

[0160] Figure 13 A schematic diagram of wavelet packet tree decomposition according to an embodiment of the present invention is shown;

[0161] Figure 14 shows a wavelet packet decomposition diagram of a signal model involved in an embodiment of the present invention;

[0162] Figure 15 The figure shows the effect of extracting the main components of the decomposed frequency bands involved in the embodiment of the present invention;

[0163] Figure 16 A flowchart of correlation noise reduction processing according to an embodiment of the present invention is shown;

[0164] Figure 17 The figure shows the effect diagram of adaptive EMD decomposition involved in the embodiment of the present invention;

[0165] Figure 18 shows a curve diagram of the permutation entropy values ​​of various components involved in an embodiment of the present invention;

[0166] Figure 19 A comparison diagram of the noise reduction effects involved in the embodiments of the present invention is shown. DETAILED DESCRIPTION

[0167] In order to more clearly understand the above-mentioned objects, features and advantages of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be noted that, in the absence of conflict, the embodiments of the present application and the features therein can be combined with each other.

[0168] In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention may also be implemented in other ways different from those described herein. Therefore, the scope of protection of the present invention is not limited to the specific embodiments disclosed below.

[0169] like Figure 1 As shown, according to an embodiment of the present invention, a power system broadband measurement method based on a combination of wavelet packets and SVMD includes the following steps:

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

[0171] S2, performs wavelet packet decomposition on the frequency information to determine the optimal wavelet basis for broadband measurement of the power system;

[0172] S3, performing multi-layer decomposition of the frequency information according to the optimal wavelet basis to obtain multiple wavelet frequency bands;

[0173] S4, performing component extraction on the wavelet frequency band, adaptively extracting harmonic components of specific frequencies using an SVMD algorithm, and obtaining harmonic components corresponding to the wavelet frequency band;

[0174] S5, use the wavelet coefficient correlation algorithm to denoise the harmonic components and obtain the denoised harmonics:

[0175]

[0176] Where c j (k) is the wavelet coefficient, j is the scale, k is the number of points, ∑ is the accumulation sign, and N is a natural number. Noise reduction is achieved by comparing the correlation of wavelet coefficients with a set threshold to determine whether to retain the current wavelet coefficient. Wavelet coefficient correlation denoising is more effective for non-stationary signals and signals containing non-local features. Compared with wavelet modulus maximum denoising, it is more effective in preserving signal details. However, it may cause signal distortion in cases where the correlation between noise and signal is weak.

[0177] However, this algorithm can easily lead to the loss of high-frequency parameters in practice. To address this issue, we introduced an improved EMD algorithm to perform a preliminary decomposition of the signal and calculate its permutation entropy. For the noisy portions of the signal, we used a correlation denoising method. The resulting reconstructed denoised signal effectively preserves broadband information, thereby improving the accuracy and reliability of signal processing.

[0178] First, traditional EMD decomposition is a step-by-step decomposition process in the time domain that targets different fluctuation trends. It decomposes the overall non-stationary signal into stationary components. These decomposition results are called intrinsic mode functions (IMFs). The detailed decomposition process is as follows: First, the local extreme values ​​of f(t) are determined through interpolation and form an envelope. The overall trend component c(t) of the waveform is averaged. After removing c(t) from the original signal, a convergence test is performed. The above decomposition process is repeated until the final signal residual r(t) is left after the screening stop condition is met, completing the decomposition process.

[0179] Through this process, it can be clearly seen that the advantage of EMD is that it is an adaptive signal decomposition method that does not require the prior assumption of the signal model or basis function. Therefore, it is applicable to various types of signals, especially nonlinear and non-stationary signals. In addition, EMD is completely based on the time domain decomposition algorithm, which makes its decomposition results complete. All output components can be directly restored by direct addition without losing the original signal. As shown in formula (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 process will not have any impact on the correlation noise reduction algorithm. In addition, this iterative decomposition and direct addition synthesis process is simple and easy to implement, and will not cause too much computational and time pressure.

[0180]

[0181] Where c i (t) is the overall trend component of the waveform, r K (t) is the signal residual, k is 0, 1, 2...

[0182] Secondly, similar to most decomposition algorithms, Fourier spectrum analysis of the IMF reveals that the frequency distribution of its components varies. This decomposition method is adaptive to the characteristics of the waveform data, and therefore does not follow the strict frequency band divisions required by wavelet decomposition. This adaptive decomposition process is precisely one of the characteristics required in this paper. However, the screening and stopping criteria for this adaptive decomposition process are difficult to determine. Typically, they are determined based on a given fixed value, which often loses some of the decomposition algorithm's adaptability. This can lead to under- or over-decomposition when decomposing signals with different characteristics, introducing errors. Therefore, an adaptive stopping mechanism is key to ensuring the stability of each component.

[0183] In EMD, the number of screening iterations is directly determined by the screening stop criterion SSC. In both "underscreening" and "overscreening" cases, the number of iterations will affect the decomposition effect. If the number of iterations is insufficient, the decomposition will be incomplete, and multiple single-component signals will be decomposed into one IMF, resulting in too many unrelated components. Too many iterations will cause the single-component signal to be decomposed into multiple IMFs, resulting in a waste of computing resources and reduced decomposition accuracy, and will also affect the efficiency of the subsequent correlation noise reduction process. Therefore, it is crucial to reduce the amount of computation to improve the noise reduction efficiency without sacrificing the decomposition accuracy. Soft SSC can monitor the screening process of EMD separation, and more importantly, it can select the optimal number of iterations. Based on 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 will always be smaller than that of the component extracted the previous time, thereby defining a signal-related function:

[0184]

[0185] Where n is the point value, k is the number of iterations, and N k Total number of iterations, E k is the kth iteration EMD decomposition component, c k (n) is the value of the n-point component of the current k-th iteration. k When the iteration is smaller than before, it can be identified as a stopping condition, and the screening process is automatically stopped according to the value of the objective function, thereby obtaining better signal decomposition performance.

[0186] The concept of soft SSC achieves adaptability in the screening process. Traditional EMD, due to fixed screening conditions, is often susceptible to noise and error in the decomposition process. Adaptive stopping conditions significantly improve the robustness of signal decomposition. Each screening process determines an optimal number of iterations to avoid generating meaningless IMFs. CPU time is proportional to both the number of iterations and the number of IMFs. This eliminates errors introduced by under-screening and wastes computing resources and time due to over-screening.

[0187] In summary, the decomposition process of adaptive stopping EMD can make up for the defects of correlation denoising and optimize its denoising effect. According to the principle of correlation denoising, signal denoising is performed in each decomposition domain of the signal, and the noise of each segment can be deeply denoised without waveform distortion. In order to reduce unnecessary calculation amount and loss of signal feature information, the components of the signal after adaptive EMD decomposition are subjected to permutation entropy detection, and correlation denoising is performed only on the noise-dominated components. The specific algorithm flow chart is as follows: Figure 16 shown.

[0188] First, the noisy signal needs to be adaptively decomposed using the EMD component method to obtain components arranged from low frequency to high frequency. For components containing a high concentration of Gaussian white noise, correlation denoising is required. Therefore, the entropy of each component is calculated and arranged according to the entropy value. Components exceeding the noise limit are set to undergo correlation denoising. Finally, the denoised signal is obtained by summing the components. Since each component has already undergone an adaptive decomposition, the number of wavelet decomposition layers can be reduced when performing correlation processing on the noisy signal, thereby saving computational effort and processing time while better preserving the true information.

[0189] In order to verify the effect of the modified correlation noise reduction algorithm, the noise reduction test was carried out on the signal containing 10dB Gaussian white noise. After experimental comparison, the db4 wavelet basis function was selected, and the high noise signal was subjected to 3-layer wavelet decomposition correlation noise reduction processing. First, the adaptive preliminary decomposition of the signal is performed as follows Figure 17 shown.

[0190] The high-frequency components in the original waveform decomposition diagram usually contain more noise. After calculating the permutation entropy value of each component, a comparison diagram of all components can be obtained, such as the permutation entropy value curve of each component. Figure 18 As shown, it can be seen that wavelet correlation denoising processing is required for the IMF1-IMF2 component.

[0191] After the improved correlation algorithm is used for noise reduction, the signal-to-noise ratio is 23.485dB and the root mean square value is 0.43. Compared with the correlation noise reduction effect of the traditional wavelet noise reduction algorithm mentioned above, the improved noise reduction algorithm retains more effective information and ensures the noise reduction effect. After the noise reduction process, each component is reconstructed to obtain the low-noise signal and the original signal. The comparison of the improved correlation noise reduction effect is shown in the figure below. Figure 19 As shown in Figure 3, the effect of the improved noise reduction algorithm is verified.

[0192] S6, decomposing the noise reduction harmonics according to the EMD algorithm to obtain the signal residual that meets the preset stopping condition;

[0193] S7, obtaining the measured data corresponding to each harmonic waveform according to the signal residual calculation.

[0194] In the above technical solution, preferably, Figure 2 As shown in Figure 2, S2 determines the optimal wavelet basis for broadband measurement of power systems, including:

[0195] S21, using the Mallat decomposition algorithm, determining the decomposition relationship between the standard orthogonal bases and the decomposition and composition relationship of the subspace coordinates corresponding to the decomposition relationship between the standard orthogonal bases;

[0196]

[0197] 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 layers, k is the number of points in each layer, Z is an integer, m is an integer in Z, ∑ is the accumulation symbol, and 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} are called low-pass filters and high-pass filters.

[0198] From this, the decomposition relationship between its standard orthogonal bases is derived, and then the decomposition and composition relationship between coordinates can be obtained, that is, the Mallat decomposition algorithm and composition algorithm:

[0199] S22, performing 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 the finite function space;

[0200]

[0201] In formulas (3)-(4), j is the number of decomposition layers, 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, c j,kis 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.

[0202] The Mallat algorithm can realize the decomposition operation of wavelet coefficients. This process is repeated until the preset decomposition level is reached, thereby constructing a wavelet transform decomposition tree, as shown in the wavelet decomposition diagram. Figure 12 As shown in Figure 2, 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 characteristics and time-frequency information of the signal, and is suitable for broadband and dynamic signal processing tasks.

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

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

[0205] The root node of the tree represents the original signal, and each node represents the output of a filter. Different decomposition layers or stopping conditions can be selected as needed to obtain decomposition results suitable for specific applications.

[0206] S23, according to the nature of the filter itself: s 2i+1 (t)=ψ j，k (t), it is deduced that the standard orthogonal bases s2i(t) and s2i+1(t) take integer translations to form two mutually orthogonal standard orthogonal systems and the subspace relationship of the corresponding standard orthogonal systems;

[0207]

[0208] Correspondingly, from the perspective of space, 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 the new symbol Sij if:

[0209]

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

[0211]

[0212] 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, 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.

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

[0214] The decomposition result is expressed as:

[0215]

[0216] The optimal wavelet basis is expressed as:

[0217]

[0218] 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, 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.

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

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

[0221] like Figure 3 As shown, in the above technical solution, preferably, S3, performing multi-layer decomposition of the frequency information according to the optimal wavelet basis, includes:

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

[0223]

[0224] 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 .

[0225] S32, based on the single period function representation, obtain the function representation of the entire harmonic waveform:

[0226]

[0227] 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 analysis time range, 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.

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

[0229] The dynamic window characteristic is expressed as: ψ(t)∈L 2 (R)

[0230]

[0231] The continuous wavelet expression that depends on the scale parameter is:

[0232]

[0233] 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:

[0234]

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

[0236] Where a is the scale parameter and b is the displacement parameter. is the complex conjugate of ψ(t).

[0237] S34, performing an inverse transformation on the continuous wavelet and its corresponding function representation to obtain the inverse function representation of the corresponding continuous wavelet:

[0238]

[0239] In the formula The wavelet transform is used to invert the wavelet packet energy. The wavelet transform shows how the signal f(t) changes around the time point b as the scale parameter a changes. S35 discretizes the inverse function to obtain multiple wavelet frequency bands.

[0240] like Figure 4 As shown, in the above technical solution, preferably, S35, discretizing the inverse function includes:

[0241] S351, defining a preset function representation of a standard orthogonal basis constituting a corresponding wavelet frequency band:

[0242]

[0243] S352: Obtain a linear representation of the corresponding wavelet frequency band and the inverse function after discretization according to the preset function representation:

[0244]

[0245] Among them, ψ(t)∈L 2 (R), then ψ(t) is an orthogonal wavelet;

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

[0247]

[0248] In formulas (17)-(19), the coefficient sequence α j，k Indicates 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 fis a continuous wavelet function.

[0249] like Figure 5 As shown, in the above technical solution, preferably, S4, performing component extraction on the wavelet frequency band, using the SVMD algorithm to adaptively extract the harmonic components of a specific frequency, and obtaining the harmonic components corresponding to the wavelet frequency band includes:

[0250] S41, using the SVMD algorithm 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:

[0251]

[0252] Where 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:

[0253]

[0254] 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, that is, the target signal to be decomposed, and the equality condition ensures that the sum of all components is the original function;

[0255] 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 α:

[0256]

[0257] (22) is defined as the augmented Lagrangian function L({u k},{ω k},λ), and solve it by iteration. The iteration formula is expressed as:

[0258]

[0259] S42, after each iteration, calculate whether the current modal component meets the following convergence conditions:

[0260]

[0261] 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 k-th mode obtained in the current iteration, is the representation of Lagrange multiplier in 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;

[0262] S43, when the decomposed modal components meet the given accuracy, K components are extracted, and the signal f(t) can be decomposed successively by the SVMD algorithm in the following manner:

[0263]

[0264] Where 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;

[0265] S44, extracting the Kth component u K (t), construct the filter impulse response β K (t) and β n (t), forming a new constraint criterion for SVMD:

[0266]

[0267] Where Q1 is the value of the filter β when the filter β is applied to the residual signal. 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 degree of overlap between the new component and the previously extracted component. The smaller Q2 is, the new component u KThe lower the overlap or aliasing between (t) and previously extracted components, the better the distinction between components, and the better the 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 adjacent extracted components u K (t) and u K-1 (t) filter impulse response;

[0268] The above filter is the key constraint to ensure that the SVMD algorithm can accurately extract each component. The expression of the filter is:

[0269]

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

[0271]

[0272] 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 K Diffusion within the frequency range; is a filter 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. 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 does 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;

[0273] S45, using the NRBO optimization algorithm to optimize the penalty factor α, and taking the minimum envelope entropy as the fitness function:

[0274]

[0275] 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, i is 0, 1, 2...;

[0276] S46, extracting accurate harmonic components from the wavelet frequency band through the above-mentioned SVMD and NRBO optimization steps, and performing component extraction on the wavelet frequency band according to the harmonic components.

[0277] In order to verify the detection accuracy and efficiency of the improved algorithm for broadband signals, a simulation analysis is performed through the constructed signal model. First, the waveform sampling frequency is set to 12.8kHz, and extraction is performed every 10 fundamental cycles. Therefore, 2560 sampling points need to be processed. After wavelet packet decomposition, it can be obtained that the harmonics within 64 times have a total of 8 frequency bands. The waveform and spectrum diagram are as shown in the signal model wavelet packet decomposition. Figure 14 shown.

[0278] 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. By extracting the main components of the four frequency bands respectively, all harmonic components can be obtained, such as Figure 15 (a), 15(b), 15(c), and 15(d).

[0279] The optimization results of the penalty parameter α based on NRBO are: 6949, 12866, 18764, and 20480, respectively. After decomposition, the harmonic waveforms of the 10 fundamental wave periods can be obtained. However, due to the endpoint effects at both ends, distortion occurs. The extraction results are corrected by removing the first and last two fundamental wave periods, and the amplitude and phase angle of each waveform are analyzed and calculated. The test results statistics are shown in the table.

[0280] Broadband algorithm test results statistics table

[0281]

[0282]

[0283] Simulation results show that the improved broadband measurement algorithm can effectively extract each harmonic component. After 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.

[0284] Simulations verified that the algorithm demonstrated excellent dynamic extraction capabilities in broadband measurements. Programming simulations based on the algorithm enabled processing of measured data, verifying the algorithm's practical processing capabilities and laying a solid foundation for further application of subsequent detection algorithms. Subsequently, three wavelet denoising algorithms were explored and, after in-depth analysis and comparison, the correlation denoising algorithm was ultimately selected.

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

[0286] An acquisition module 10 is configured to acquire frequency information corresponding to circuit signals at various measurement points in the circuit system;

[0287] The wavelet decomposition module 20 is configured to perform wavelet packet decomposition on the frequency information to determine the optimal wavelet basis for broadband measurement of the power system;

[0288] The multi-layer decomposition module 30 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;

[0289] a component extraction module 40 configured to perform modal decomposition based on the wavelet frequency band using the SVMD algorithm and sequentially extract each harmonic component, and optimize the penalty parameter α of the SVMD algorithm while taking envelope entropy minimization as the objective function to reduce harmonic aliasing;

[0290] The noise reduction module 50 is configured to perform noise reduction processing on the harmonic components using a wavelet coefficient correlation algorithm to obtain noise-reduced harmonics;

[0291] The harmonic decomposition module 60 is configured to decompose the noise reduction harmonics according to the EMD algorithm to obtain a signal residual that meets a preset stop condition;

[0292] The calculation module 70 is configured to calculate the measurement data corresponding to each harmonic waveform based on the signal residual.

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

[0294] The orthogonal basis decomposition unit 21 is configured 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 the standard orthogonal bases using the Mallat decomposition algorithm;

[0295]

[0296] (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 layers, k is the number of points in each layer, Z is an integer, m is an integer in Z, ∑ is the accumulation symbol, and 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} are called low-pass filters and high-pass filters.

[0297] The spatial transformation unit 22 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 the finite function space;

[0298]

[0299] In formulas (3)-(4), j is the number of decomposition layers, 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, 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.

[0300] The subspace decomposition unit 23 is configured to: s 2i+1 (t)=ψ j，k (t), it is deduced that the standard orthogonal bases s2i(t) and s2i+1(t) take integer translations to form two mutually orthogonal standard orthogonal systems and the subspace relationship of the corresponding standard orthogonal systems;

[0301]

[0302] 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, 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.

[0303] The optimal wavelet basis generating unit 24 is configured to obtain the decomposition result and the optimal wavelet basis of the corresponding wavelet packet according to the standard orthogonal system and its subspace relationship;

[0304] The decomposition result is expressed as:

[0305]

[0306] The optimal wavelet basis is expressed as:

[0307]

[0308] 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, 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.

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

[0310] The wavelet transform unit 31 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:

[0311]

[0312] 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 .

[0313] The waveform generating unit 32 is configured to obtain a functional representation of the entire harmonic waveform according to a single period functional representation:

[0314]

[0315] 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 analysis time range, 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.

[0316] The continuous wavelet generating unit 33 is configured to perform a continuous wavelet transform on the function representation of the entire harmonic waveform according to the dynamic window characteristics to obtain a scale-parameter-dependent continuous wavelet and a function representation corresponding to the continuous wavelet;

[0317] The dynamic window characteristic is expressed as: ψ(t)∈L 2 (R)

[0318]

[0319] The continuous wavelet expression that depends on the scale parameter is:

[0320]

[0321] 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.

[0322] The function representation corresponding to the continuous wavelet is:

[0323]

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

[0325] Where a is the scale parameter and b is the displacement parameter. is the complex conjugate of ψ(t).

[0326] The inverse transform unit 34 is configured to perform an inverse transform on the continuous wavelet and its corresponding function representation to obtain an inverse function representation of the corresponding continuous wavelet:

[0327]

[0328] In the formula To invert the wavelet packet energy, the wavelet transform shows how the signal f(t) changes around time point b as the scale parameter a changes. The wavelet transform shows how the signal f(t) changes around time point b as the scale parameter a changes.

[0329] The discrete processing unit 35 is configured to discretize the inverse function to obtain a plurality of wavelet frequency bands.

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

[0331] The wavelet frequency band generating subunit 351 is configured to define a preset function representation of a standard orthogonal basis constituting a corresponding wavelet frequency band:

[0332]

[0333] The frequency band linear function generating subunit 352 is configured to obtain a linear representation of the corresponding wavelet frequency band and the discretized inverse function according to a preset function representation:

[0334] Among them, ψ(t)∈L 2 (R), then ψ(t) is an orthogonal wavelet;

[0335] The orthogonal coefficient generating subunit 353 is configured to obtain the wavelet frequency band according to the coefficient sequence of the preset standard orthogonal basis:

[0336]

[0337] 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 continuous wavelet function.

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

[0339] The SVMD principal component extraction unit 41 is configured to perform modal decomposition based on the wavelet frequency band using the SVMD algorithm and extract each harmonic component in sequence:

[0340] First, the SVMD algorithm is used to perform successive modal decomposition on the wavelet frequency band to obtain the component signals u k (t), the component signal u k The expression of (t) is:

[0341]

[0342] Where 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:

[0343]

[0344] 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, that is, the target signal to be decomposed, and the equality condition ensures that the sum of all components is the original function;

[0345] 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 α:

[0346]

[0347] (22) is defined as the augmented Lagrangian function L({u k},{ω k},λ), and solve it by iteration. The iteration formula is expressed as:

[0348]

[0349] Each time an iteration is completed, the current modal component is calculated to see if it meets the following convergence conditions:

[0350] 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 k-th mode obtained in the current iteration, is the representation of Lagrange multiplier in 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;

[0351] When the decomposed modal components meet the given accuracy, K components are extracted and the signal f(t) can be decomposed successively by the SVMD algorithm in the following way:

[0352]

[0353] Where 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;

[0354] 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:

[0355]

[0356] Where Q1 is the value of the filter β when the filter β is applied to the residual signal. 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 degree of 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, and the better the 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 adjacent extracted components u K (t) and u K-1 (t) filter impulse response;

[0357] The above filter is the key constraint to ensure that the SVMD algorithm can accurately extract each component. The expression of the filter is:

[0358]

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

[0360]

[0361] 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 K Diffusion within the frequency range; is a filter 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. 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 does 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;

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

[0363]

[0364] 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, i is 0, 1, 2...;

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

[0366] Based on the above Figure 5 and Figure 6 The method shown, accordingly, an embodiment of the present application also provides a computer-readable storage medium, on which a computer program is stored, which, when executed by a processor, implements the steps of the power system broadband measurement method based on the combination of wavelet packets and SVMD in any of the above embodiments.

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

[0368] Based on the above Figure 5 and Figure 6 The method shown, and Figure 7 In the virtual device embodiment shown, in order to achieve the above-mentioned purpose, 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 computer programs; the processor is used to execute the computer program to implement the steps of the power system wideband measurement method based on the combination of wavelet packets and SVMD in any of the above-mentioned embodiments.

[0369] Optionally, the computer device may further include a user interface, a network interface, a camera, a radio frequency (RF) circuit, a sensor, an audio circuit, a Wi-Fi module, etc. The user interface may include a display, an input unit such as a keyboard, etc., and the optional user interface may also 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.

[0370] Those skilled in the art will understand that the computer device structure provided in this embodiment does not constitute a limitation on the computer device, and may include more or fewer components, or a combination of certain components, or different component arrangements.

[0371] The storage medium may also include an operating system and a network communication module. An operating system is a program that manages and stores the hardware and software resources of a computer device, supporting the execution of information processing programs and other software and / or programs. The network communication module facilitates communication between components within the storage medium, as well as with other hardware and software within the physical device.

[0372] Throughout this specification, terms such as "one embodiment," "some embodiments," and "specific embodiments" mean that the specific features, structures, materials, or characteristics described in conjunction with that embodiment or example are included in at least one embodiment or example of the present invention. In this specification, schematic representations of these terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in any one or more embodiments or examples.

[0373] The above description is only the preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations.

[0374] Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the invention shall be

[0375] It is included in the protection scope of the present invention.

Claims

1. A power system broadband measurement method based on the combination of wavelet packets 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 optimal wavelet basis for broadband measurement of the power system; Step 3: performing 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 the preset stop condition; Step 7: Calculate the measured data corresponding to each harmonic waveform based on the signal residual; The component extraction of the wavelet frequency band in the above step 4 includes: First, the SVMD algorithm is used to perform successive modal decomposition on the wavelet frequency band to obtain the component signals u k (t), the component signal u k The expression of (t) is: Where 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, that is, 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 k-th mode obtained in the current iteration, is the representation of Lagrange multiplier in 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 by the SVMD algorithm in the following way: Where 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 value of the filter β when the filter β is applied to the residual signal. 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 degree of 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, and the better the 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 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 K Diffusion within the frequency range; is a filter 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. 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 does 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: 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, i is 0, 1, 2...; Through the above SVMD and NRBO optimization steps, accurate harmonic components are extracted from the wavelet frequency band, and component extraction is performed on the wavelet frequency band according to the harmonic components.

2. The power system broadband measurement method 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 composition 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 layers, k is the number of points in each layer, Z is an integer, m is an integer in Z, ∑ is the accumulation symbol, and 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 low-pass filter and 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 layers, 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, 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), derive the orthonormal bases s2i(t) and s2i+1(t) which, after integer translation, form two mutually orthogonal orthonormal systems and the subspace relationship corresponding to the orthonormal systems; 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, 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, 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.

3. The power system broadband measurement method 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 analysis time range, 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 functional representation of the entire harmonic waveform according to a dynamic window characteristic to obtain a scale-parameter-dependent continuous wavelet and a functional representation corresponding to the continuous wavelet; The dynamic window characteristic is expressed as: ψ(t)∈L 2 (R) The continuous wavelet expression 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 wavelet transform window is a rectangular window centered at the 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 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 change 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 power system broadband measurement method according to claim 3, characterized in that: Discretizing the inverse function includes: The preset function expression constituting the standard orthogonal basis corresponding to the wavelet frequency band is defined as follows: 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, ψ(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 continuous wavelet function.

5. A power system broadband measurement system for the power system broadband measurement method based on the combination of wavelet packets and SVMD as claimed in claim 1, characterized in that: include: an acquisition module 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 configured to perform multi-layer decomposition on the frequency information according to the optimal wavelet basis to obtain a plurality of wavelet frequency bands; An SVMD principal component extraction module is configured to perform modal decomposition based on the wavelet frequency band using an SVMD algorithm and extract each harmonic component in sequence; An NRBO optimization module is configured to optimize a penalty parameter α of the SVMD algorithm and use envelope entropy minimization as an objective function to reduce harmonic aliasing; a noise reduction module 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 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 based on the signal residual.

6. The power system broadband measurement system according to claim 5, 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 layers, k is the number of points in each layer, Z is an integer, m is an integer in Z, ∑ is the accumulation symbol, and 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 low-pass filter and 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 the finite function space; In formulas (3)-(4), j is the number of decomposition layers, 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, 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), derive the orthonormal bases s2i(t) and s2i+1(t) which, after integer translation, form two mutually orthogonal orthonormal systems and the subspace relationship corresponding to the orthonormal systems; 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, 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 of a corresponding wavelet packet and an optimal wavelet basis 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, 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.

7. 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 functional representation of the entire harmonic waveform according to the single period functional 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 analysis time range, 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 configured to perform a continuous wavelet transform on the functional representation of the entire harmonic waveform according to a dynamic window characteristic to obtain a scale-parameter-dependent continuous wavelet and a functional representation corresponding to the continuous wavelet; The dynamic window characteristic is expressed as: ψ(t)∈L 2 (R) The continuous wavelet expression 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 wavelet transform window is a rectangular window centered at the 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.

8. The power system broadband measurement system according to claim 7, characterized in that: The discrete processing unit includes: The wavelet frequency band generating subunit is configured to define a preset function representation of the standard orthogonal basis corresponding to the 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, ψ(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 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 continuous wavelet function.

9. The power system broadband measurement system according to claim 5, characterized in that: The SVMD principal component extraction module includes: Variational constraint unit: used to construct the constrained variational problem of SVMD based on 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 use the minimum envelope entropy as the objective function; The harmonic extraction unit is used to extract the harmonic components corresponding to the frequency band according to the modal components obtained by SVMD decomposition, and output them to the noise reduction module.

Images