A Broadband Signal Harmonic Analysis Method Based on VMD-SG Filtering and Improved HHT

Abstract

This invention discloses a broadband signal harmonic analysis method based on VMD-SG filtering-improved HHT. The broadband signal harmonic analysis method includes the following steps: acquiring a discrete sampling sequence of the broadband signal to be analyzed, and performing sliding frame processing on the discrete sampling sequence to obtain multiple sets of continuous single-frame signal sequences; performing variational mode decomposition on each single-frame signal sequence to obtain k IMF components; performing SG filtering on each IMF component to obtain a filtered mode; performing improved HHT processing on the filtered mode to identify reconstruction parameters; and reconstructing the target broadband signal based on the identified reconstruction parameters to complete the harmonic analysis of the broadband signal.

Description

Technical Field

[0001] This invention relates to the field of power system harmonic analysis technology, and in particular to a broadband signal harmonic analysis method based on VMD-SG filtering-improved HHT. Background Technology

[0002] Driven by the "dual carbon" goal, with the grid connection of new energy sources and the integration of numerous power electronic devices (such as frequency converters, inverters, charging piles, flexible DC transmission, etc.) and nonlinear loads, actual signals often contain fundamental, harmonic, interharmonic, and transient attenuation components, and are frequently contaminated by noise, exhibiting characteristics of broadband, non-stationarity, and complexity. Frequency components of signals such as harmonics and interharmonics may experience amplitude changes, frequency shifts, or transient abrupt changes in a short period, posing significant challenges to power quality monitoring and harmonic parameter assessment. However, for broadband signals from new energy sources, traditional frequency domain analysis methods such as Fast Fourier Transform (FFT) are limited by spectral leakage and the picket fence effect, making it difficult to handle non-stationary signals. While time-frequency analysis methods such as wavelet transform possess time-frequency analysis capabilities, their time-frequency resolution is constrained by fixed basis functions, resulting in insufficient adaptability. Empirical Mode Decomposition (EMD), although possessing adaptive capabilities, suffers from three inherent defects: mode aliasing, noise sensitivity, and endpoint effects, leading to problems of mode aliasing, noise, and harmonic mixing. Summary of the Invention

[0003] In view of the above-mentioned prior art, the present invention provides a broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT, which mainly solves the technical problems existing in the background art.

[0004] To achieve the above objectives, the technical solution of this invention is implemented as follows: This invention discloses a broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT. The broadband signal harmonic analysis method includes the following steps: The discrete sampling sequence of the broadband signal to be analyzed is obtained, and the discrete sampling sequence is subjected to sliding frame division processing to obtain multiple sets of continuous single-frame signal sequences. Variational mode decomposition is performed on each frame of the single-frame signal sequence to obtain k IMF components; SG filtering is performed on each IMF component to obtain the filtered mode; The filter mode is subjected to improved HHT processing to identify the reconstruction parameters. The target broadband signal is then reconstructed based on the identified reconstruction parameters to complete the harmonic analysis of the broadband signal.

[0005] Optionally, a discrete sampling sequence of the broadband signal to be analyzed is obtained, and the discrete sampling sequence is subjected to sliding frame division processing to obtain multiple sets of continuous single-frame signal sequences, specifically including: Based on the characteristics of broadband signals, construct the expression for broadband signals:

[0006] With preset sampling frequency The above-mentioned broadband signal to be analyzed is sampled to obtain a discrete sampling sequence of the broadband signal:

[0007] The discrete sampling sequence y[n] is subjected to sliding frame processing. Based on the preset window length N and frame step size H, the m-th frame signal sequence is generated, as shown in the following expression: +

[0008] In the formula, M is the total number of signal components. Let be the amplitude of the i-th signal component. Let i be the frequency of the i-th signal component. Let be the initial phase of the i-th signal component, n(t) be Gaussian white noise, n[n] be a discrete-domain Gaussian white noise sequence, and n represent the n-th sampling point. Let be the attenuation factor of the i-th signal component. The sampling period.

[0009] Optionally, variational mode decomposition is performed on each frame of the single-frame signal sequence to obtain k IMF components, as shown in the following equation:

[0010] In the formula, Let r[n] be the k-th IMF component obtained from the decomposition, and r[n] be the residual.

[0011] Optionally, SG filtering is performed on each IMF component to obtain the filtered mode, specifically including: Perform a fast Fourier transform on each IMF component to obtain the spectrum of the corresponding IMF component, and calculate the center frequency and bandwidth of the IMF component based on the spectrum. Based on the center frequency and bandwidth of the IMF component, determine the filtering parameters of the corresponding IMF component; Based on the determined filtering parameters, in length Within the window, perform p-order polynomial least squares fitting to obtain the filtered IMF components; Amplitude compensation is performed by comparing the energy of the filtered modes to obtain the final pure IMF component.

[0012] Optionally, the SG filtering results can be improved by HHT processing, specifically including: Constructing the mirrored signal sequence based on the endpoint mirroring extension method. ; For the extended signal sequence Perform improved EMD sieving process and determine sieving stop based on sieving convergence index. Output the current IMF component when the index requirements are met; otherwise, update the residual signal. The output IMF component is truncated back to the original signal window length to obtain an IMF component of the same length as the original single-frame signal sequence; The HHT transform is performed on the IMF components that are of the same length as the original single-frame signal sequence to obtain the analytical signals of each IMF component.

[0013] Optionally, for the extended signal sequence Implement improved EMD screening processes, specifically including: The residual signal is initialized by assigning the signal sequence after mirror extension to the initial residual signal of the outer layer. An outer-layer sieving loop is performed on the initial residual signal. During the outer-layer sieving loop, all local maxima and local minima of the current residual signal are extracted, the total number of local extrema is counted, and the total number of local extrema is compared with a threshold. If the threshold is exceeded, an inner-layer sieving loop is performed. The final residual signal obtained during the outer layer screening cycle is assigned as the initial candidate signal of the inner layer. The inner layer initial candidate signal is subjected to an inner layer screening loop, and all local maxima and local minima of the current candidate signal are extracted to form a set of local maxima and a set of local minima; Interpolation fitting is performed on the set of local maxima and the set of local minima respectively to obtain upper and lower envelopes of the same length as the candidate signal, and the envelope mean of the upper and lower envelopes at each discrete sampling point is calculated respectively. Candidate signals are updated based on the envelope mean, and the screening convergence index is calculated based on the updated candidate signals.

[0014] Optionally, screening stop determination can be performed based on screening convergence indices, specifically including: The calculated screening convergence index is compared with the preset convergence threshold. If the screening convergence index is less than the preset convergence threshold, the current inner screening cycle is terminated, and the candidate signal obtained by iteration is used as the current IMF component. If the screening convergence index is greater than or equal to the preset convergence threshold, the residual signal is updated and the improved EMD screening process is repeated.

[0015] Optionally, the reconstruction parameters include amplitude parameters, phase parameters, and frequency parameters, which are obtained as follows: Based on the demodulation of the analytical signal, the instantaneous amplitude sequence, instantaneous phase sequence, and instantaneous frequency sequence of the corresponding signal components are obtained; The median of the instantaneous amplitude sequence within a single frame signal window is taken as the amplitude parameter corresponding to the IMF component; the instantaneous phase sequence is differentially processed to obtain the instantaneous frequency sequence, and the median of the instantaneous frequency sequence within a single frame signal window is taken as the frequency parameter corresponding to the IMF component; the instantaneous phase at the center sampling point of the single frame signal window is taken as the phase parameter corresponding to the IMF component.

[0016] Optionally, based on the amplitude parameters, frequency parameters, and phase parameters of each identified IMF component, and combined with the attenuation factor of the corresponding IMF component, a reconstruction sequence of the target broadband signal is constructed to complete the reconstruction of the target broadband signal.

[0017] Optionally, the attenuation factor can be solved using the least squares method based on the amplitude parameters of the IMF components.

[0018] The beneficial effects of this invention are as follows: by using sliding frame processing, long-sequence non-stationary broadband signals are transformed into multiple sets of short-sequence approximately stationary signals. Combined with the adaptive non-recursive decomposition characteristics of variational mode decomposition, complex broadband signals are decomposed into multiple narrowband intrinsic mode function components with independent center frequencies. This fundamentally suppresses the inherent mode aliasing defects of traditional empirical mode decomposition. At the same time, it overcomes the technical bottleneck of traditional fast Fourier transform being constrained by spectral leakage and picket fence effect and unable to effectively process non-stationary signals. It also solves the problem that wavelet transform time-frequency resolution is limited by fixed basis functions and is not adaptable to complex broadband signal scenarios. It achieves accurate separation of different frequency band components in broadband signals. Based on the adaptive matching of SG filter parameters according to the spectral characteristics of each intrinsic mode function component, and in conjunction with the amplitude compensation operation after filtering, differentiated and accurate noise reduction of components in different frequency bands is achieved. This solves the problems of noise and harmonic component coupling and the easy occurrence of filtered waves or insufficient filtering by fixed filter parameters in traditional methods. While effectively filtering out Gaussian white noise interference in the field, it completely preserves the amplitude, phase and transient change characteristics of the signal, greatly improves the signal-to-noise ratio of broadband signals in complex noise environments, and effectively improves the defect of traditional empirical mode decomposition being highly sensitive to noise. By using the endpoint mirror extension method associated with the SG filter window length, combined with the improved EMD sieving process of double-layer loop, the inherent endpoint effect of the traditional HHT transform is specifically suppressed, the pollution of the instantaneous feature calculation of the entire signal range by the envelope fitting distortion at both ends of the signal is eliminated, and the full-dimensional parameter accurate identification of steady-state harmonics and transient attenuation components is achieved, solving the problem of insufficient accuracy of the traditional method in identifying transient broadband component parameters. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only preferred embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0020] Figure 1 This is a schematic diagram of the overall process of the broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT in the embodiments of this application; Figure 2 This is the signal graph after adding 20dB Gaussian white noise; Figure 3 This is a fitting graph of a broadband signal after noise reduction; Figure 4 This is a detailed flowchart illustrating an embodiment of this application. Detailed Implementation

[0021] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used in this specification is for the purpose of describing particular embodiments only and is not intended to limit the invention. In the following description, the expression "some embodiments" refers to a subset of all possible embodiments; however, it should be understood that "some embodiments" can be the same subset or different subsets of all possible embodiments and can be combined with each other without conflict.

[0022] In the following description, numerous specific details are set forth in order to provide a more thorough understanding of the invention. However, it will be apparent to those skilled in the art that the invention can be practiced without one or more of these details. In other instances, certain technical features well-known in the art have not been described in order to avoid obscuring the invention.

[0023] It should be understood that the present invention can be embodied in various forms and should not be construed as being limited to the embodiments set forth herein. Rather, providing these embodiments will make the disclosure thorough and complete, and will fully convey the scope of the invention to those skilled in the art. Furthermore, the terminology used herein is intended only to describe particular embodiments and is not intended to limit the invention. When used herein, the singular forms “a,” “an,” and “the” are also intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the terms “compose” and / or “comprising,” when used in this specification, identify the presence of the stated features, integers, steps, operations, elements, and / or components, but do not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and / or groups. When used herein, the term “and / or” includes any and all combinations of the associated listed items.

[0024] It should also be noted that when an element is referred to as being "fixed to" another element, it can be directly attached to the other element or there may be an intervening element. When an element is referred to as being "connected to" another element, it can be directly connected to the other element or there may be an intervening element. The terms "vertical," "horizontal," "inner," "outer," "left," "right," and similar expressions used herein are for illustrative purposes only and do not represent the only possible implementation.

[0025] To fully understand this invention, a detailed structure will be presented in the following description to illustrate the technical solution proposed by this invention. Optional embodiments of the invention are described in detail below; however, in addition to these detailed descriptions, the invention may have other embodiments.

[0026] Please refer to the attached document. Figure 1 and Figure 4 This application provides a broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT. The broadband signal harmonic analysis method includes the following steps: S1. Obtain the discrete sampling sequence of the broadband signal to be analyzed, and perform sliding frame division processing on the discrete sampling sequence to obtain multiple sets of continuous single-frame signal sequences. Specifically, the discrete sampling sequence of the broadband signal to be analyzed is obtained, and the discrete sampling sequence is subjected to sliding frame division processing to obtain multiple sets of continuous single-frame signal sequences, specifically including: Based on the characteristics of broadband signals, construct the expression for broadband signals:

[0027] With preset sampling frequency The above-mentioned broadband signal to be analyzed is sampled to obtain a discrete sampling sequence of the broadband signal:

[0028] The discrete sampling sequence y[n] is subjected to sliding frame processing. Based on the preset window length N and frame step size H, the m-th frame signal sequence is generated, as shown in the following expression: +

[0029] In the formula, M is the total number of signal components. Let be the amplitude of the i-th signal component. Let i be the frequency of the i-th signal component. Let be the initial phase of the i-th signal component, n(t) be Gaussian white noise, n[n] be the discrete domain Gaussian white noise sequence, and n represent the discrete sampling number. Let be the attenuation factor of the i-th signal component. The sampling period.

[0030] It should be noted that the attenuation factor for the steady-state fundamental and steady-state harmonic components in a power system... The value is 0, which represents the attenuation factor for transient attenuation components generated during the grid connection of new energy sources and the switching of power electronic devices. It is a non-zero real number, thereby achieving a unified representation of all types of broadband signals in the power system.

[0031] In some implementations, the window length N is typically a power of 2, such as 256, 512, or 1024, to facilitate efficient execution of subsequent fast Fourier transform operations. The frame step size H is less than the window length N. Preferably, the frame step size H is half of the window length N, i.e., a sliding frame method with a 50% overlap rate is used to reduce the truncation error of inter-frame signals and improve the continuity of time-frequency analysis of non-stationary broadband signals.

[0032] S2. Perform variational mode decomposition (VMD) on each frame of the single-frame signal sequence to obtain k mode functions (IMFs). Specifically, variational mode decomposition is performed on each frame of the single-frame signal sequence to obtain k IMF components and residuals, as shown in the following equation:

[0033] In the formula, Let r[n] be the k-th IMF component obtained by decomposition, r[n] be the residual signal, and K be the number of modes (IMFs) obtained by VMD decomposition.

[0034] It should be noted that variational mode decomposition, as an adaptive non-recursive signal decomposition method, can effectively suppress the inherent defects of mode aliasing compared with traditional empirical mode decomposition methods, achieve accurate separation of different frequency band components in broadband signals, and adapt to the decomposition requirements of complex broadband signals containing harmonics, interharmonics and transient attenuation components in new energy grid-connected scenarios.

[0035] Preferably, before performing variational mode decomposition (VMD) processing, the core decomposition parameters need to be configured in advance, including key parameters such as the number of mode decompositions K, the second-order penalty factor, and the fidelity coefficient. In some implementations, for practical engineering scenarios of broadband harmonic analysis in power systems, the number of mode decompositions K can be matched and set according to the total number of fundamental, harmonic, and interharmonic components contained in the broadband signal to be analyzed, with a typical value range of 3 to 8. The typical value of the second-order penalty factor is 2000, and the typical value of the fidelity coefficient is 0. This ensures the separation effect of components in different frequency bands while taking into account the computational efficiency of subsequent processing. It should be noted that variational mode decomposition constructs and solves a constrained variational problem to decompose a single-frame signal sequence into multiple narrowband IMF components with independent center frequencies and finite bandwidths. Each IMF component is an amplitude-frequency modulated (AM-FM) signal that meets the requirements of a single-component signal.

[0036] S3. Perform SG filtering on each IMF component to obtain the filtered mode, which specifically includes: Perform a fast Fourier transform on each IMF component to obtain the spectrum of the corresponding IMF component, and calculate the center frequency and bandwidth of the IMF component based on the spectrum. Based on the center frequency and bandwidth of the IMF component, determine the filtering parameters of the corresponding IMF component; Based on the determined filtering parameters, in length Within the window, perform p-order polynomial least squares fitting to obtain the filtered IMF components; Amplitude compensation is performed by comparing the energy of the filtered modes to obtain the final pure IMF component.

[0037] Specifically, the time-domain discrete sequence of the k-th IMF component obtained from VMD decomposition The time-domain to frequency-domain conversion is performed using Fast Fourier Transform (FFT), yielding the frequency-domain spectrum sequence corresponding to the IMF component, expressed as:

[0038] in, Here, f is the Fast Fourier Transform operator, and f represents the discrete frequency points. Let be the frequency domain spectrum sequence of the k-th IMF component.

[0039] Based on the obtained spectral sequence, the power spectral density sequence of this IMF component is further calculated, and the expression is:

[0040] In the formula, For modulo operation, Let f be the power spectral density value of the k-th IMF component at the corresponding discrete frequency point f. It should be noted that the power spectral density can accurately characterize the distribution of the signal energy of the IMF component at different frequency points. Compared with the feature calculation based on the amplitude spectrum directly, the calculation method based on the power spectrum can further improve the anti-interference capability of frequency domain feature extraction and better adapt to the broadband signal processing scenario of power system with on-site noise.

[0041] Based on the power spectral density sequence obtained above, the center frequency of the IMF component is calculated using the centroid method, which characterizes the core frequency band position of the IMF component and provides a core reference for subsequent filter parameter matching. The formula for calculating the center frequency is as follows:

[0042] Based on the center frequency and power spectral density sequence obtained from the aforementioned calculations, the root mean square equivalent bandwidth of the IMF component is calculated. This characterizes the bandwidth distribution of the IMF component and provides a direct quantitative basis for the adaptive adjustment of the filter window length. The formula for calculating the bandwidth is as follows:

[0043] This formula calculates the second moment of the power spectrum relative to the center frequency to obtain the statistically significant equivalent bandwidth of the IMF component. It can accurately characterize the frequency band dispersion of the IMF component. The larger the bandwidth value, the wider the frequency band distribution of the IMF component and the richer the transient characteristics of the signal. The corresponding matching SG filter window length is shorter, so as to completely preserve the transient characteristics of the signal while filtering out noise. The smaller the bandwidth value, the more concentrated the frequency band distribution of the IMF component and the higher the signal smoothness. The corresponding matching SG filter window length is longer, so as to achieve a better smoothing and noise reduction effect. This enables differentiated and adaptive filtering processing of IMF components in different frequency bands, solving the problem of insufficient filtering or filtering that is easy to occur with fixed filtering parameters.

[0044] Based on the center frequency and bandwidth of the IMF component obtained from the aforementioned calculations, the SG filter parameters for the corresponding IMF component are determined. These SG filter parameters include the filter window length and the order of the fitted polynomial. The formula for the filter window length is:

[0045] The formula for fitting the order of the polynomial is:

[0046] In the formula, For the first k The SG filter window length for each mode, To fit the order of the polynomial to SG, and An empirical coefficient for adjusting the window length; The lower and upper limits of the allowed window length, It is a very small positive number, used to avoid the denominator being zero.

[0047] Based on the aforementioned determined filtering parameters, SG filtering is performed on the corresponding IMF components, in length... Perform within the window Polynomial least squares fitting is used to obtain the filtered modes, where the expression for polynomial least squares fitting is:

[0048] The final formula for calculating the filtered IMF components is as follows:

[0049] In the formula These are the convolution coefficients of the SG filter.

[0050] After completing the SG filtering process, the amplitude compensation of the filtered IMF component is performed by comparing the modal energy before and after filtering to eliminate the amplitude attenuation problem introduced by the filtering operation, thus obtaining the final pure IMF component. The calculation formula for amplitude compensation is as follows:

[0051] Signal energy is calculated using the sum of squares of the time-domain sequence to ensure energy conservation of the IMF components before and after filtering, thus avoiding systematic deviations in subsequent harmonic amplitude parameter identification due to amplitude attenuation caused by filtering. It should be noted that the purified IMF components after amplitude compensation effectively filter out noise interference from the original broadband signal while fully preserving the time-domain characteristics, amplitude information, and phase information of each frequency band component. This provides a high-precision signal foundation for subsequent improvements to HHT time-frequency analysis and full-dimensional harmonic parameter identification.

[0052] S4. The filter mode is subjected to improved HHT processing to identify the reconstruction parameters. The target broadband signal is reconstructed based on the identified reconstruction parameters to complete the harmonic analysis of the broadband signal.

[0053] In some optional implementations, the SG filtering results are subjected to improved HHT processing, specifically including: Based on the endpoint mirror extension method, a mirror-extended signal sequence is constructed, and the extension length is correlated with the window length of the aforementioned SG filter to match the filtering characteristics of different frequency band components. structure Length of extended sequence :

[0054]

[0055] In the formula, Extension length, The sequence after mirror extension.

[0056] It should be noted that the endpoint effect of the traditional HHT transform can cause severe distortion in the envelope fitting at both ends of the signal, thereby contaminating the calculation results of instantaneous frequency and instantaneous phase. This invention, by associating the extension length with the SG filter window length, matches different extension lengths for IMF components of different frequency bands. This can suppress the endpoint effect while avoiding redundant calculations introduced by excessive extension, and fully adapt to the processing requirements of each component of broadband signals.

[0057] Improved EMD sieving is performed on the mirror-extended signal sequence. First, the initialization operation of the outer sieving loop is performed, and the mirror-extended signal sequence is assigned the initial residual signal. ,Right now:

[0058] After initialization, the outer screening loop is entered. In each round of the outer screening loop, all local maxima and local minima of the current residual signal are extracted, and the total number of local extrema is counted. This total number is compared with a preset threshold of 2. If the total number is greater than or equal to 2, the inner screening loop is entered to extract the single component IMF with the corresponding index. If the total number is less than 2, it is determined that there are no effective decomposable single component signals in the current residual signal, and the entire outer screening loop is terminated.

[0059] During the inner layer screening cycle, the inner layer initialization operation is performed first, and the residual signal corresponding to the outer layer screening cycle in the current round is initialized. Assigning values ​​to the initial candidate signal ,Right now:

[0060] Then extract the current candidate signal. All local maxima and local minima form a corresponding set of extreme points. Cubic spline interpolation is then applied to both sets of local maxima and local minima to obtain an upper envelope of equal length to the candidate signal. and lower envelope Calculate the mean envelope values ​​of the upper and lower envelopes at each discrete sampling point:

[0061] The candidate signal update for this iteration is completed by subtracting the mean envelope value of the corresponding sampling point from the current candidate signal. The update formula is as follows:

[0062] Where j is the number of iterations in the inner sieving cycle.

[0063] After updating the candidate signals for this inner iteration, the screening convergence index SD corresponding to this iteration is calculated. Screening stop is determined based on the screening convergence index. The calculation formula for the screening convergence index SD is as follows:

[0064] The calculated sieving convergence index SD is compared with the preset convergence threshold. The comparison is performed; if the screening convergence index SD is less than the preset convergence threshold... If the current inner-layer screening cycle is terminated, the candidate signals obtained in this iteration will be processed. The expression for the i-th IMF component extracted in the current outer layer screening cycle is: IMF .

[0065] If the screening convergence index SD is greater than or equal to the preset convergence threshold Then update the residual:

[0066] Increment the number of inner layer sieving iterations j by 1, and use the candidate signal updated in this iteration as the initial candidate signal for the next inner layer iteration. Repeat the operations of extreme point extraction, envelope fitting, envelope mean calculation, candidate signal update and convergence index comparison in the inner layer sieving loop until the convergence termination condition of the inner layer sieving loop is met.

[0067] Preferred convergence threshold The typical value range is 0.2 to 0.3, so as to ensure the single-component characteristics of the IMF component while avoiding signal waveform distortion caused by excessive screening.

[0068] After completing the improved EMD screening process, the output IMF component is truncated back to the original signal window length, and the extended sequence portions at both ends are removed to obtain an IMF component of the same length as the original single-frame signal sequence. This eliminates the influence of the extended sequence on subsequent parameter identification and ensures the consistency of signal timing.

[0069] It should be noted that the truncated IMF components not only suppress the envelope fitting distortion caused by the endpoint effect through mirror extension, but also completely match the time window and number of sampling points of the original single-frame signal, thus fully preserving the time domain characteristics and frequency domain information of the signal.

[0070] Furthermore, based on the improved EMD sieving results, a Hilbert transform is performed to obtain the analytical signals corresponding to each IMF component. The calculation formula for the analytical signals is as follows:

[0071] Based on the analytical signal obtained from the solution, the instantaneous amplitude sequence, instantaneous phase sequence, and instantaneous frequency sequence of the corresponding signal components are demodulated. The calculation formulas for the instantaneous amplitude sequence and the instantaneous phase sequence are as follows:

[0072] In the formula, Let be the instantaneous amplitude of the k-th IMF component in the j-th frame signal at the n-th sampling point. Let be the instantaneous phase of the k-th IMF component in the j-th frame signal at the n-th sampling point. This is the analytic signal obtained after Hilbert transforming the corresponding IMF components. For the modulo operation of complex numbers, For complex number argument operations, k is the index of the IMF component, j is the frame index of the single frame signal, and n is the discrete sampling index.

[0073] It should be noted that the instantaneous amplitude characterizes the time-varying amplitude of the signal component at the corresponding sampling time, and the instantaneous phase characterizes the time-varying phase of the signal component at the corresponding sampling time. This formula can be used to directly demodulate the instantaneous time-domain characteristics of the signal from the analytical signal, providing a direct data source for the identification of the core harmonic parameters of amplitude, frequency, and phase. In practical engineering implementation, the phase angle calculation is usually implemented using the four-quadrant arctangent function, which can cover the full phase range of [0, 2π], avoid the calculation error caused by phase truncation, and fully adapt to the full phase identification requirements of broadband signals in power systems.

[0074] The instantaneous frequency sequence is obtained as follows: Based on the instantaneous phase sequence obtained from the aforementioned solution, the phase derivative in the continuous domain is approximated through discrete difference operations to obtain the instantaneous frequency sequence of the corresponding IMF component, thereby characterizing the time-varying characteristics of the signal frequency. The formula for obtaining the instantaneous frequency sequence is as follows:

[0075] In the formula, Let be the instantaneous frequency of the k-th IMF component in the j-th frame signal at the n-th sampling point. The signal sampling period, Let be the signal sampling frequency, and satisfy . =1 / π is the mathematical constant of a circle.

[0076] Based on the instantaneous amplitude sequence, instantaneous phase sequence, and instantaneous frequency sequence obtained from the aforementioned solution, the harmonic reconstruction parameters are identified. These reconstruction parameters include amplitude, phase, and frequency parameters. The specific identification process is as follows: the median of the instantaneous amplitude sequence within a single-frame signal window is taken as the amplitude parameter corresponding to the IMF component; the median of the instantaneous frequency sequence within a single-frame signal window is taken as the frequency parameter corresponding to the IMF component; and the instantaneous phase at the center sampling point of the single-frame signal window is taken as the phase parameter corresponding to the IMF component. The identification formulas for the amplitude and frequency parameters are as follows:

[0077]

[0078] In the formula, To obtain the median operator, This is the amplitude parameter corresponding to the k-th IMF component. This is the frequency parameter corresponding to the k-th IMF component.

[0079] The formula for identifying the phase parameter is:

[0080] in is the index of the center sampling point of a single frame signal window, and N is the total number of sampling points in the original single frame signal sequence.

[0081] It should be noted that using the median operator to identify amplitude and frequency parameters can effectively suppress outliers and noise interference in the instantaneous parameter calculation process. Compared with the traditional mean calculation method, it has stronger anti-interference ability and robustness. Furthermore, using the instantaneous phase of the sampling point at the center of the window as the phase identification result can minimize the pollution of phase calculation by the endpoint effect, further improve the identification accuracy of phase parameters, and fully adapt to the broadband signal analysis scenario of power system with on-site noise.

[0082] Furthermore, based on the amplitude parameters and instantaneous amplitude sequences obtained from the aforementioned identification, the attenuation factor of the corresponding IMF component is solved using the least squares method, thus completing the full-dimensional identification of the reconstructed parameters. The following matrix is ​​first defined during the calculation:

[0083] In the formula, is the attenuation factor corresponding to the k-th IMF component.

[0084] In the discrete domain, let According to the model Take the logarithm:

[0085] The solution obtained using the least squares method is as follows:

[0086] The attenuation factor of the corresponding signal component can be obtained based on the solution results. .

[0087] After completing the full-dimensional identification of the amplitude, frequency, phase, and attenuation factors of all IMF components, a parameter component set for all signal components is constructed. Based on the identified full-dimensional reconstruction parameters, combined with the unified expression of the broadband signal constructed in step S1, a reconstruction sequence of the target broadband signal is constructed to complete the reconstruction of the target broadband signal.

[0088] After signal reconstruction is completed, quantitative evaluation indicators are established using root mean square error (RMSE) and signal-to-noise ratio (SNR) to verify the signal reconstruction accuracy and harmonic analysis effect. The relevant calculation formulas are as follows:

[0089]

[0090] In the formula, For the reconstructed signal, The original signal is N, and the signal length is N.

[0091] The following experimental example is given to verify the effectiveness of the invention method: Establish a wideband signal with a sampling frequency of 3000Hz and a sampling time of 0.3s, such as... Figure 2 As shown, its wideband signal expression is as follows:

[0092] In the formula: The original signal, Added 20dB of Gaussian white noise.

[0093] The signal was decomposed into four Integrated Mode Factors (IMFs) using VMD. Then, these IMFs were subjected to SG filtering. Finally, the improved HHT identification parameters were used to reconstruct the signal, and the effectiveness of the invention was verified using root mean square error (RMSE) and signal-to-noise ratio (SNR). Table 1 shows the RMSE and SNR of different methods. The traditional HHT method has poor noise reduction, while the improved HHT method has a higher SNR, reaching 23.33687, with a RMSE of only 0.1806.

[0094] Table 1. Root mean square error values ​​for different methods

[0095] The improved HHT was used to identify the amplitude, frequency, phase, and attenuation factor information of the signal. The results are shown in Table 2 and as follows: Figure 3 Fitted plot.

[0096] Table 2 Fitting errors for each parameter

[0097] Table 2 shows that the average error in amplitude is 0.643%, and the average error in frequency is 0.353%. Figure 3 As can be seen, the curve obtained by fitting almost overlaps with the original signal curve, further verifying the accuracy of the method.

[0098] The above are merely specific embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. The scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A wideband signal harmonic analysis method based on VMD-SG filtering-improved HHT, characterized in that, The broadband signal harmonic analysis method includes the following steps: The discrete sampling sequence of the broadband signal to be analyzed is obtained, and the discrete sampling sequence is subjected to sliding frame division processing to obtain multiple sets of continuous single-frame signal sequences. Variational mode decomposition is performed on each frame of the single-frame signal sequence to obtain k IMF components; SG filtering is performed on each IMF component to obtain the filtered mode; The filter mode is subjected to improved HHT processing to identify the reconstruction parameters. The target broadband signal is then reconstructed based on the identified reconstruction parameters to complete the harmonic analysis of the broadband signal.

2. The method according to claim 1, wherein the method is a wideband signal harmonic analysis method based on VMD-SG filtering-improved HHT. The discrete sampling sequence of the broadband signal to be analyzed is obtained, and the discrete sampling sequence is subjected to sliding frame division processing to obtain multiple sets of continuous single-frame signal sequences, specifically including: Based on the characteristics of broadband signals, construct the expression for broadband signals: at a predetermined sampling frequency The wideband signal to be analyzed is sampled to obtain a discrete sampling sequence of the wideband signal: The discrete sampling sequence y[n] is subjected to sliding frame processing. Based on the preset window length N and frame step size H, the m-th frame signal sequence is generated, as shown in the following expression: + In the formula, M is the total number of signal components. Let be the amplitude of the i-th signal component. Let i be the frequency of the i-th signal component. Let be the initial phase of the i-th signal component, n(t) be Gaussian white noise, n[n] be a discrete-domain Gaussian white noise sequence, and n represent the n-th sampling point. Let be the attenuation factor of the i-th signal component. The sampling period.

3. The broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT according to claim 1, characterized in that, Variational mode decomposition is performed on each frame of the single-frame signal sequence to obtain k IMF components, as shown in the following equation: In the formula, Let r[n] be the k-th IMF component obtained from the decomposition, and r[n] be the residual.

4. The broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT according to claim 3, characterized in that, Each IMF component is subjected to SG filtering to obtain the filtered mode, specifically including: Perform a fast Fourier transform on each IMF component to obtain the spectrum of the corresponding IMF component, and calculate the center frequency and bandwidth of the IMF component based on the spectrum. Based on the center frequency and bandwidth of the IMF component, determine the filtering parameters of the corresponding IMF component; Based on the determined filtering parameters, in length Within the window, perform p-order polynomial least squares fitting to obtain the filtered IMF components; Amplitude compensation is performed by comparing the energy of the filtered modes to obtain the final pure IMF component.

5. A broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT according to claim 4, characterized in that, The SG filtering results are then subjected to improved HHT processing, specifically including: Constructing the mirrored signal sequence based on the endpoint mirroring extension method. ; For the extended signal sequence Perform improved EMD sieving process and determine sieving stop based on sieving convergence index. Output the current IMF component when the index requirements are met; otherwise, update the residual signal. The output IMF component is truncated back to the original signal window length to obtain an IMF component of the same length as the original single-frame signal sequence; The HHT transform is performed on the IMF components that are of the same length as the original single-frame signal sequence to obtain the analytical signals of each IMF component.

6. The broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT according to claim 5, characterized in that, For the extended signal sequence Implement improved EMD screening processes, specifically including: The residual signal is initialized by assigning the signal sequence after mirror extension to the initial residual signal of the outer layer. An outer-layer sieving loop is performed on the initial residual signal. During the outer-layer sieving loop, all local maxima and local minima of the current residual signal are extracted, the total number of local extrema is counted, and the total number of local extrema is compared with a threshold. If the threshold is exceeded, an inner-layer sieving loop is performed. The final residual signal obtained during the outer layer screening cycle is assigned as the initial candidate signal of the inner layer. The inner layer initial candidate signal is subjected to an inner layer screening loop, and all local maxima and local minima of the current candidate signal are extracted to form a set of local maxima and a set of local minima; Interpolation fitting is performed on the set of local maxima and the set of local minima respectively to obtain upper and lower envelopes of the same length as the candidate signal, and the envelope mean of the upper and lower envelopes at each discrete sampling point is calculated respectively. Candidate signals are updated based on the envelope mean, and the screening convergence index is calculated based on the updated candidate signals.

7. A broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT according to claim 6, characterized in that, The screening stop determination is based on the screening convergence index, specifically including: The calculated screening convergence index is compared with the preset convergence threshold. If the screening convergence index is less than the preset convergence threshold, the current inner screening cycle is terminated, and the candidate signal obtained by iteration is used as the current IMF component. If the screening convergence index is greater than or equal to the preset convergence threshold, the residual signal is updated and the improved EMD screening process is repeated.

8. A broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT according to claim 7, characterized in that, The reconstructed parameters include amplitude parameters, phase parameters, and frequency parameters, and their acquisition process is as follows: Based on the demodulation of the analytical signal, the instantaneous amplitude sequence, instantaneous phase sequence, and instantaneous frequency sequence of the corresponding signal components are obtained; The median of the instantaneous amplitude sequence within a single frame signal window is taken as the amplitude parameter corresponding to the IMF component; the instantaneous phase sequence is differentially processed to obtain the instantaneous frequency sequence, and the median of the instantaneous frequency sequence within a single frame signal window is taken as the frequency parameter corresponding to the IMF component; the instantaneous phase at the center sampling point of the single frame signal window is taken as the phase parameter corresponding to the IMF component.

9. A broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT according to claim 8, characterized in that, Based on the amplitude, frequency, and phase parameters of each IMF component obtained by identification, and combined with the attenuation factor of the corresponding IMF component, a reconstruction sequence of the target broadband signal is constructed to complete the reconstruction of the target broadband signal.

10. A broadband signal harmonic analysis method based on VMD-SG filtering and improved HHT according to claim 9, characterized in that, The attenuation factor is solved using the least squares method based on the amplitude parameters of the IMF components.

Images