Combined multi-scale electric energy quality signal preprocessing method based on singular spectrum analysis and local mean decomposition

Abstract

The invention discloses a united multi-scale power quality signal preprocessing method based on singular spectrum analysis and local mean decomposition. The united multi-scale power quality signal preprocessing method comprises the following steps: decomposing a singular spectrum component and a residual component from a disturbance signal based on singular spectrum decomposition; decomposing a series of PF components and residual components from the disturbance signal based on local mean decomposition; combining the two groups of components based on scale matching and weighted stacking to obtain a multi-scale mixed component; and analyzing the merged multi-scale mixed components based on S transformation to obtain a corresponding time-frequency graph. According to the method, the dependency relationship of the time sequence can be more powerfully captured, and more obvious signal characteristics and time sequence rules are provided for a subsequent classification model.

Description

Technical Field

[0001] The present invention relates to the field of power quality disturbance signal preprocessing, and in particular to a joint multi-scale power quality signal preprocessing method based on singular spectrum analysis and local mean decomposition. Background Art

[0002] With the rapid growth of new energy sources such as wind and solar energy, new power quality problems have also been brought about. The types of power quality disturbances are increasing, including traditional voltage sags, transient overvoltages, and new disturbances caused by the access of new energy sources. These disturbances are often hidden and complex, and are difficult to effectively identify through traditional monitoring methods. Therefore, the large amount of data, signal noise interference, and limitations of analysis methods are the current difficulties.

[0003] In order to achieve accurate identification of disturbance signals, most existing methods use PQDs identification methods based on deep learning, avoiding the subjectivity of traditional manual feature selection and instead utilizing the powerful automatic learning and feature mining capabilities of neural networks for feature extraction. However, most of these methods do not perform data preprocessing. Directly inputting the original signal into the classification model will contain a large amount of irrelevant or redundant information such as noise. This information will affect the performance of the model and even interfere with the learning features of the model, causing important features to be ignored, thereby reducing the recognition accuracy of the model in practical applications.

[0004] In recent years, in order to achieve data preprocessing of disturbance signals, researchers have explored a variety of multi-scale analysis methods, such as wavelet transform, variational mode decomposition (VMD) and singular spectrum decomposition (SSD), to better capture the local characteristics and global trends of signals. These methods help to improve the performance of signal preprocessing, especially when processing non-stationary and nonlinear signals. However, the limitations of the above methods are also quite obvious. Although wavelet transform can provide multi-scale analysis, it requires manual selection of appropriate wavelet basis functions when processing signals, and the subjectivity of manual selection cannot be avoided; variational mode decomposition (VMD) requires prior knowledge of the characteristics of the signal to set the number of modes and frequency bandwidth; the single use of singular spectrum decomposition (SSD) may cause mode aliasing, affecting the accuracy of subsequent analysis. Summary of the invention

[0005] In order to solve the above technical problems, the present invention provides a joint multi-scale power quality signal preprocessing method based on singular spectrum analysis and local mean decomposition. This method can more effectively capture the dependencies of time series and provide more obvious signal characteristics and timing laws for subsequent classification models.

[0006] The technical solution adopted by the present invention is:

[0007] The combined multi-scale power quality signal preprocessing method based on singular spectrum analysis and local mean decomposition includes the following steps:

[0008] Step 1: Decompose the singular spectrum components and residual components from the disturbance signal based on singular spectrum decomposition;

[0009] Step 2: Decompose a series of PF components and residual components from the disturbance signal based on local mean decomposition;

[0010] Step 3: Combine the two components based on scale matching and weighted superposition to obtain a multi-scale mixed component;

[0011] Step 4: Analyze the combined multi-scale mixed components based on S transform to obtain the corresponding time-frequency diagram.

[0012] In step 1, first, a two-dimensional signal matrix is ​​constructed by splicing a one-dimensional time series with its own lagged version, and the two-dimensional signal matrix is ​​subjected to singular value decomposition to obtain the product of three matrices: a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix. The singular values ​​are arranged in descending order, reflecting the contribution of each singular vector to the original matrix information. The first few largest singular values ​​and their corresponding singular vectors are selected to extract the main trends and periodic components in the time series. These components are also singular spectral components, and each singular spectral component contains part of the signal information.

[0013] The step 1 comprises the following steps:

[0014] Step 1.1: Process the original disturbance signal V(t) into a discrete signal h(n), find the mean z of the initial discrete signal h(n), n = 1, 2, ..., N, N is the length of the sequence after the signal is discretized, and de-mean all discrete sequence point values ​​to obtain a new input signal h′(n);

[0015] h′(n)=h(n)-z

[0016] Step 1.2: Calculate the residual component v at the jth iteration of the singular spectrum decomposition SSD j Power spectral density of (n):

[0017]

[0018] in, represents Fourier transform, f is the frequency, is the residual component v at the jth iteration j (n), and thus estimate the frequency f corresponding to the maximum peak of the power spectrum max :

[0019]

[0020] Among them, arg max f is a function to find the maximum value; if the number of iterations j = 1 and f s is the set sampling frequency, the embedding dimension M is set to N; when the number of iterations j>1 or j=1 and does not satisfy When the embedding dimension M is set to 1.2×(f max / f s );

[0021] Step 1.3: Select the embedding dimension M according to the method in step 1.2, and construct the h′(n) obtained in step 1.1 into an M×(N-M+1) signal matrix H. The matrix H is expressed as:

[0022]

[0023] In the formula, represents the first sliding window of the signal sequence h'(n), i.e., from h'(1) to h'(1+M-1); Represents the second sliding window of the signal sequence h'(n), i.e., from h'(2) to h'(2+M-1); represents the Mth sliding window of the signal sequence h'(n), i.e., from h'(M) to h'(M+M-1); T represents the sequence transposition;

[0024] Construct the jth low-frequency subsequence I from high frequency to low frequency j , I j It is expressed as:

[0025] I j =(i1,i2,…i p )

[0026] Where i1,i2,…i p Respectively represent the row index in the trajectory matrix H;

[0027] Reconstruct the corresponding subsequence by performing a diagonal average operation on the matrix:

[0028] For each diagonal in the matrix H, calculate the average value of its elements. The diagonal in the matrix H can be expressed as H i,j , where i and j are the row index and column index respectively, and i+j=k, k is a constant representing the index of the diagonal. For each diagonal k, calculate its average value

[0029] Among them, D k is the set of all elements on the diagonal k, |D k| is the number of elements on the diagonal k. The average value of all diagonals Arrange them in the order of the diagonal lines to obtain a reconstructed subsequence of the time series.

[0030] Step 1.4: Perform singular value decomposition on the matrix H obtained in step 1.3;

[0031]

[0032] In the formula, U represents an M×M left singular vector matrix, whose column vectors are the left singular vectors of matrix H, which are the orthogonal basis vectors in the singular value decomposition of H; V represents an (N-M+1)×(N-M+1) right singular vector matrix. Its column vectors are the right singular vectors of matrix H, which are also the orthogonal basis vectors in the singular value decomposition of H; T represents transpose; ∑ represents an M×(N-M+1) diagonal matrix, whose diagonal elements are the singular values ​​of matrix H,

[0033] U=[u1,u2,…u M ];

[0034] Where U represents an M×M left singular vector matrix; u1,u2,…u M They respectively represent the column vector in the matrix U, that is, the left singular vector.

[0035] V=[v1,v2,…v W ];

[0036] Where V represents a (N-M+1)×(N-M+1) right singular vector matrix; v1,v2,…v W They respectively represent the nth column vector in the matrix V, that is, the nth right singular vector.

[0037] W = N - M + 1;

[0038] Where W represents the number of right singular vectors, N represents the length of the original sequence, and M represents the embedding dimension;

[0039] Σ=diag(σ1,σ2,…σ r );

[0040] Among them, σ i represents the non-zero singular value of H, i = 1, 2, ..., r, r represents the number of non-zero singular values, and i max =min(M,W), let H i =σ i u i v i T ,u irepresents the i-th column vector in the matrix U, that is, the i-th left singular vector; v i T Represents the transpose of the i-th column vector in the matrix V.

[0041] Then the signal matrix can be expressed as:

[0042] H=H1+H2+…+H q

[0043] In the formula, H1 represents the first subsequence extracted by singular value decomposition, H2 represents the second subsequence extracted by singular value decomposition, and H q represents the qth subsequence extracted by singular value decomposition;

[0044] Step 1.5: After each iteration, a new signal component g is separated from the signal j (n) and a residual component v j+1 (n), after each iteration, the residual component v is calculated j+1 The decomposition process is terminated when the calculated mean square error is less than the preset threshold x.

[0045]

[0046] Where h(n) represents the original signal;

[0047]

[0048] Where mse(j) represents the mean square error between the residual component and the original sequence after the jth iteration, and N represents the sequence length.

[0049] The step 2 decomposes the disturbance signal into a series of PF components and residual components based on local mean decomposition. The singular spectral components and residual components separated from the disturbance signal are obtained in step 1, and step 2 also decomposes the disturbance signal to obtain the singular spectral components and residual components separated from the disturbance signal. Extract the PF(t) component and the monotone residual component The complex signal is decomposed into a series of functions with local mean characteristics through an iterative process.

[0050] Step 2 includes the following steps:

[0051] Step 2.1: Find the initial signal All extreme points n i , calculate the local mean α of two adjacent extreme points i and the local envelope estimate β i ;

[0052]

[0053] Where n i represents the i-th extreme point in the signal, n i+1 Indicates the i+1th extreme point in the signal;

[0054] Step 2.2: Use the exponentially weighted moving average (EWMA) to adjust the local mean α obtained in step 2.1. i and the local envelope estimate β i Perform data smoothing to obtain the local mean function α 11 (t) and the local envelope estimation function β 11 (t);

[0055] a 11 (t) = λa i (t)+(1-λ)a 11 (t-1)

[0056] β 11 (t) = λβ i (t)+(1-λ)β 11 (t-1)

[0057] Among them, λ is the smoothing parameter, a i (t) is the local mean at time point t, a 11 (t-1) is the local mean after smoothing at time t-1; the initial signal Subtract the smoothed local mean function α 11 (t) Get the zero mean signal C 11 (t);

[0058] C 11 (t) = x(t) - α 11 (t)

[0059] The zero-mean signal retains the main features of the original signal, but its mean is 0, which means that the energy of the signal is more evenly distributed in time.

[0060] Step 2.3: From the modulated signal, that is, the zero-mean signal C 11 (t) to extract the original information, namely the modulation function d 11 (t);

[0061]

[0062] Step 2.4: Repeat steps 2.1 and 2.2 to iterate and obtain α ik (t) and β ik (t), and then calculate the pure modulation function ρ 1k The value of (t);

[0063]

[0064] Where, d 1k-1 (t) represents the PF component obtained in k-1 iterations; α 1k (t) represents the local mean function obtained in the kth iteration; C 1k-1 (t) represents the envelope signal obtained in the k-1th iteration; β 1k-1 (t) represents the local envelope estimation function obtained in the k-1th iteration; k represents the number of iterations; and n represents the sequence length.

[0065] Step 2.5: Iterate all local envelope estimation functions β obtained in step 2.4 1k (t) is multiplied to obtain the envelope signal β1(t);

[0066]

[0067] In the formula, β 12 (t) represents the local envelope estimation function obtained in the second iteration; β 1n (t) represents the local envelope estimation function obtained in the nth iteration; β 1q (t) represents the local envelope estimation function; q represents the index.

[0068] Let the envelope signal β1(t) and the first pure modulation function ρ 11 (t) The first product function PF1 can be obtained by multiplying the product;

[0069] PF1=β1(t)ρ 11 (t)

[0070] Step 2.6: From the initial signal The PF1 obtained in step 2.5 is stripped off to obtain a new initial signal Repeat steps 2.1 to 2.6 for k times until Shows monotonicity of function;

[0071]

[0072] In the formula, represents the signal after the first PF component is stripped from the initial signal, that is, the signal after the first iteration, represents the initial signal; PF1(t) represents the first modal component, i.e., the PF component;

[0073] represents the signal after the second PF component is stripped from the signal after the first iteration; PF2(t) represents the second modal component; represents the signal after stripping the kth PF component from the signal after the k-1th iteration; represents the signal after the k-1th iteration; PF k (t) represents the k-1th modal component.

[0074] After k cycles, the initial signal can be The decomposition output is k PF components and a residual monotonic component

[0075]

[0076] Wherein, PF0(t) represents the pth modal component, p represents the index, and p = 1, 2, …, k.

[0077] In the step 3, a set of SSD components and LMD components can be obtained through step 1 and step 2, and the scale characteristics of the two groups of components are analyzed. The main frequency range of each component is determined using a one-dimensional discrete Fourier transform; by setting a frequency standard, the high-frequency, medium-frequency, and low-frequency intervals are divided based on the frequency threshold, and the two groups of components are classified into the set frequency categories. After classification, the components of the same frequency category in the two sets are compared, and scale matching and adjustment are performed. If a high-frequency component decomposed by SSD is similar to a medium-frequency component decomposed by LMD in frequency, but there is a slight difference, one of the components is filtered or resampled to make its frequency range match that of the other component. After scale matching, the two groups of components are weighted superimposed, and the weight is determined based on the contribution of each component to the original signal;

[0078] Step 3 includes the following steps:

[0079] Step 3.1: Analyze the data characteristics of the two components and transform the signal data into the frequency domain through one-dimensional discrete Fourier transform to extract the frequency range of each component:

[0080]

[0081] In the formula, X(e jω ) represents the discrete Fourier transform result of the signal; F[x(n)] represents the discrete Fourier transform result of the signal; e -jωn represents a complex exponential function; x(n) is a discrete sequence, that is, two components obtained by SSD decomposition and LMD decomposition, and the frequency components of different components are extracted according to Fourier transform.

[0082] Step 3.2: By setting the frequency standard, the high frequency, medium frequency and low frequency intervals are divided based on the frequency threshold, and the two groups of components are classified into the set frequency categories. Two thresholds f are set. mid and f low , f mid It is the boundary between mid-frequency and high-frequency, f lowIt is the boundary between low frequency and medium frequency. After classification, the components of the same frequency category in the two sets are compared for scale matching and adjustment. If a high-frequency component decomposed by SSD is similar to a medium-frequency component decomposed by LMD in frequency but slightly different, one of the components is filtered or resampled to match its frequency range with the other component.

[0083] Step 3.3: Determine the sparsification weight coefficient of each component by calculating the signal-to-noise ratio. The signal-to-noise ratio is an indicator that measures the ratio of characteristic information to noise in the signal. The weighted superposition sparsification weight is assigned according to the signal-to-noise ratio. The higher the signal-to-noise ratio of the component, the greater its weight. The calculation formula of the signal-to-noise ratio is as follows:

[0084]

[0085] In the formula, SNR represents signal-to-noise ratio, P signal Represents the power of the signal, P noise represents the power of the noise, s(n) represents the original signal, n(m) represents the noise signal, M represents the length of the signal sequence, N represents the length of the noise, m represents the mean of the signal, and n represents the mean of the noise.

[0086] Step 3.4: Determine the weighted superposition sparsification weight coefficient w by the signal-to-noise of each component calculated in step 3.3 i , and then adjust the sparse weights to ensure that the sum of all weights is 1 to meet the superposition requirements, and use the adjusted weights w i ′For the components x in the same scale in the two groups i (t) is weighted, and then all weighted components are superimposed together to form the final signal y(t). The mathematical expression of this process is as follows:

[0087]

[0088] Among them, Normalized SNR i represents the normalized signal-to-noise ratio of the i-th component, τ is the preset signal-to-noise ratio threshold, and the threshold is set to sparse the weights.

[0089] In step 4, the S transform uses a window whose width varies inversely with the frequency, thereby effectively providing high time resolution for the high-frequency part of the signal and high frequency resolution for the low-frequency part of the signal. Most complex power quality events are non-stationary. The S transform effectively extracts features using a Gaussian window with scalable conversion. The S transform time-frequency analysis is performed on the combined signal y(t), and the S transform is defined as follows:

[0090]

[0091] Among them, S(τ,f) is the S transform of the signal x(t), t is time, f is frequency, w(τ-t,f) is the Gaussian window, and τ represents the parameter controlling the time axis position of the Gaussian window.

[0092] The S transform projects the time-frequency domain waveform of PQD (Power Quality Disturbances) into a two-dimensional space. The time-frequency matrix diagram contains the frequency and amplitude information of the disturbance. Taking the voltage sag disturbance as an example, Figure 6(a) to Figure 6(c) shown.

[0093] The present invention provides a multi-scale power quality disturbance signal preprocessing method based on singular spectrum decomposition and local mean decomposition, and the technical effects are as follows:

[0094] 1) Advantages of step 1 of the present invention: decomposing the singular spectrum component and the residual component from the disturbance signal based on singular spectrum decomposition can overcome the defect of selecting the embedding dimension based on experience, and adaptively decompose the signal according to the distribution from high frequency to low frequency.

[0095] 2) Advantages of step 2 of the present invention: Based on local mean decomposition, the disturbance signal is decomposed into a series of PF components and residual components, which can adaptively extract the PF components corresponding to different signals and effectively reduce the influence of endpoint effects on the decomposition results.

[0096] 3) Advantages of step 3 of the present invention: Based on scale matching and weighted superposition, the two groups of components are merged to obtain a mixed multi-scale component, which can more comprehensively capture the characteristics and dynamic changes of the original signal by merging components from different sources, thereby enhancing the overall performance of the signal. The optimized mixed multi-scale component is more suitable for subsequent disturbance classification tasks.

[0097] 4) Advantages of step 4 of the present invention: Based on the S transform, the combined multi-scale mixed components are analyzed to obtain the corresponding time-frequency diagram, which can integrate the information of all modal components after the joint decomposition, so that the time-frequency characteristics of the power quality disturbance signal can be further enhanced. BRIEF DESCRIPTION OF THE DRAWINGS

[0098] The present invention will be further described below in conjunction with the accompanying drawings and examples;

[0099] Figure 1 This is the flow chart of singular spectrum decomposition.

[0100] Figure 2 This is the local mean decomposition flow chart.

[0101] Figure 3 This is the diagram of the amplitude disturbance separation process.

[0102] Figure 4 This is the result diagram of additional disturbance detection.

[0103] Figure 5 This is the result graph of main frequency point extraction.

[0104] Figure 6(a) is a voltage sag waveform diagram;

[0105] Figure 6(b) is an S-transform feature analysis diagram;

[0106] Figure 6(c) is a three-dimensional display of the s-transform. DETAILED DESCRIPTION

[0107] The power quality disturbance preprocessing method based on singular spectrum decomposition and local mean decomposition includes the following steps:

[0108] Step 1: Decompose the singular spectrum components and residual components from the disturbance signal based on singular spectrum decomposition;

[0109] Step 2: Decompose a series of PF components and residual components from the disturbance signal based on local mean decomposition;

[0110] Step 3: Multi-scale components obtained by combining singular spectrum decomposition and local mean decomposition based on scale matching and weighted superposition;

[0111] Step 4: Analyze the combined multi-scale mixed components based on S transform to obtain the corresponding time-frequency diagram;

[0112] In the step 1, decomposing the singular spectrum component and the residual component from the disturbance signal based on the singular spectrum decomposition comprises the following steps:

[0113] Step S1.1: According to the IEEE standard, seven common single power quality disturbance signal models are established, and the parameters are randomly generated within the range. The signal model and standard parameters of the single power quality disturbance voltage are shown in Table 1.

[0114] Table 1 Single power quality disturbance model

[0115]

[0116] The original disturbance signal V(t) is converted into a discrete signal h(n), and the mean z of the original discrete signal h(n) is calculated, where n=1, 2,…, N, and N is the length of the sequence after the signal is discretized. All discrete sequence point values ​​are de-averaged to obtain a new input signal h′(n).

[0117] h′(n)=h(n)-z

[0118] Step S1.2: Calculate the residual component v at the jth iteration of SSD j The power spectrum density of (n) is used to estimate the frequency f corresponding to the maximum peak of the power spectrum. max If the number of iterations j = 1 and (f s is the set sampling frequency), then the embedding dimension M is set to N; when the number of iterations j>1 or j=1 and does not satisfy When the embedding dimension M is set to 1.2×(f max / f s ).

[0119] Step S1.3: Select the embedding dimension M according to the method in step S1.2, and construct h′(n) obtained in step 1.1 into an M×N signal matrix H. The matrix H is expressed as:

[0120]

[0121] Construct the jth low-frequency subsequence I from high frequency to low frequency j , I j It is expressed as:

[0122] I j =(i1,i2,…i0)

[0123] By changing the matrix H ij Perform diagonal averaging to reconstruct the corresponding subsequence.

[0124] Step S1.4: Perform singular value decomposition on the matrix H obtained in step S1.3.

[0125]

[0126] U=[u1,u2,…u M ]

[0127] V=[v1,v2,…v W ]

[0128] W=N-M+1

[0129] Σ=diag(σ1,σ2,…σ r )

[0130] where σ i (i=1,2,…,r) represents the non-zero singular value of H, and i max =min(M,W), let H i =σ i u i v i T , then the signal matrix can be expressed as:

[0131] H=H1+H2+…+H q

[0132] Step S1.5: After each iteration, a new signal component g is separated from the signal j (n) and a residual component v j+1 (n), after each iteration, the residual component v is calculated j+1 The decomposition process is terminated when the calculated mean square error is less than the preset threshold x.

[0133]

[0134] In the step 2, a series of PF components and residual components are decomposed from the disturbance signal based on local mean decomposition, including the following steps:

[0135] Step S2.1: Find the initial signal All extreme points n i , calculate the local mean α of two adjacent extreme points i and the local envelope estimate β i .

[0136]

[0137] Step S2.2: Use exponentially weighted moving average (EWMA) to calculate the local mean α obtained in step S2.1. i and the local envelope estimate β i Perform data smoothing to obtain the local mean function α 11 (t) and the local envelope estimation function β 11 (t), the initial signal Subtract the smoothed local mean function α 11 (t) Get the zero mean signal C 11 (t),

[0138] C 11 (t) = x(t) - α 11 (t)

[0139] The zero-mean signal retains the main features of the original signal, but its mean is 0, which means that the energy of the signal is more evenly distributed in time.

[0140] Step S2.3: From the modulated signal, that is, the zero-mean signal C 11 (t) to extract the original information, namely the modulation function d 11 (t),

[0141]

[0142] Step S2.4: Repeat steps S2.1 and S2.2 to iterate and obtain α ik (t) and β ik(t), and then calculate the pure modulation function ρ 1k The value of (t).

[0143]

[0144] Step S2.5: All local envelope estimation functions β obtained by iterating step S2.4 are 1k (t) is multiplied to obtain the envelope signal β1(t),

[0145]

[0146] Let the envelope signal β1(t) and the first pure modulation function ρ 11 (t) can be multiplied to obtain the first product function PF1.

[0147] PF1=β1(t)ρ 11 (t)

[0148] Step S2.6: From the initial signal The PF1 obtained in step 2.5 is stripped off to obtain a new initial signal Repeat steps S2.1 to S2.6 for k times until The function exhibits monotonicity.

[0149]

[0150] After k cycles, the initial signal can be The decomposition output is k PF components and a residual monotonic component

[0151]

[0152] In step 3, the two components are combined based on scale matching and weighted superposition to obtain a mixed multi-scale component, including the following steps:

[0153] Step S3.1: Analyze the data characteristics of the two components, and transform the signal data into the frequency domain through one-dimensional discrete Fourier transform to extract the frequency range of each component:

[0154]

[0155] Where x(n) is a discrete sequence, that is, two components obtained by SSD decomposition and LMD decomposition, and the frequency components of different components are extracted according to Fourier transform.

[0156] Step S3.2: By setting the frequency standard, the high frequency (above 2.56kHz), medium frequency (100Hz-2.56kHz), and low frequency (45-100Hz) intervals are divided based on the frequency threshold, and the two groups of components are classified into the set frequency categories. After classification, the components of the same frequency category in the two sets are compared, and scale matching and adjustment are performed. If a high-frequency component decomposed by SSD is similar to a medium-frequency component decomposed by LMD in frequency, but there is a slight difference, one of the components is filtered or resampled to make its frequency range match that of the other component.

[0157] Step S3.3: Determine the sparsification weight coefficient of each component by calculating the signal-to-noise ratio. The signal-to-noise ratio is an indicator that measures the ratio of characteristic information to noise in a signal. The weighted superposition sparsification weight is assigned according to the signal-to-noise ratio. The higher the signal-to-noise ratio of a component, the greater its weight. The calculation formula for the signal-to-noise ratio is as follows:

[0158]

[0159] Step S3.4: Determine the weighted superposition sparsification weight coefficient w based on the signal-to-noise of each component calculated in step S3.3 i , and then adjust the sparse weights to ensure that the sum of all weights is 1 to meet the superposition requirements, and use the adjusted weights w i ′For the components x in the same scale in the two groups i (t) is weighted, and then all weighted components are superimposed together to form the final signal y(t). The mathematical expression of this process is as follows:

[0160]

[0161] Among them, τ is the preset signal-to-noise ratio threshold, and the threshold is set to sparse the weights.

[0162] In the step 4, the combined multi-scale mixed components are analyzed based on the S transform to obtain the corresponding time-frequency diagram, including the following steps:

[0163] Step S4.1: Perform S transform time-frequency analysis on the combined signal y(t) using Matlab. The S transform is defined as follows:

[0164]

[0165] Among them, S(δ,f) is the result of S transformation of the combined signal y(t), t represents time, f represents frequency, and w(δ-t,f) represents the Gaussian window function, δ represents the time shift parameter of the Gaussian window function. The S transformation projects the time-frequency domain waveform of PQD into two-dimensional space, and the real frequency matrix diagram contains the frequency, amplitude, phase, frequency domain and other information of the disturbance.

[0166] In a specific embodiment of the present invention, the proposed method is used to pre-process multiple types of composite disturbance signals, and all disturbance waveforms are generated according to the mathematical model in Table 1.

[0167] Figure 1 This is the singular spectrum decomposition process, corresponding to step 1. First, a two-dimensional signal matrix is ​​constructed by concatenating the one-dimensional disturbance signal with its own lagged version. The signal matrix is ​​subjected to singular value decomposition to obtain the product of three matrices: the left singular vector matrix, the singular value diagonal matrix, and the right singular vector matrix. The singular values ​​are arranged in descending order, reflecting the contribution of each singular vector to the original matrix information. The first few largest singular values ​​and their corresponding singular vectors are selected to extract the main trends and periodic components in the time series.

[0168] Figure 2 This is the local mean decomposition process, corresponding to step 2. In step 1, the singular spectral component and residual component separated from the disturbance signal are obtained, and step 2 also decomposes the disturbance signal, extracts the PF(t) component and the monotonic residual component from the disturbance signal, and decomposes the complex signal into a series of functions with local mean characteristics through an iterative process.

[0169] Figure 3 The original signal of “voltage interruption + harmonics” without preprocessing and the time-frequency diagram after S transformation.

[0170] Figure 4 The original signal of “voltage sag + harmonics” without preprocessing and the time-frequency diagram after S transformation.

[0171] Figure 5 The time-frequency diagram of the signal after single decomposition and joint mode decomposition. Figure 3 and Figure 4 It can be seen from the time-frequency diagram of the original signal without preprocessing that the useful features in the signal, such as the amplitude and frequency components of the disturbance, are significantly enhanced after signal preprocessing. By comparing the time-frequency diagrams obtained by the joint decomposition with those obtained by the single decomposition, it can be seen that the time-frequency diagram obtained by the joint decomposition has more obvious differences in texture shape and distribution, which can more effectively reflect different types of disturbances and provide a more reliable sample basis for the next step of feature recognition.

[0172] From the perspective of signal processing, the present invention uses a joint decomposition method for signal preprocessing, which has the advantage of complementary features. At the same time, SSD and LMD are used to decompose the fault signal to obtain multi-layer modal components that can better reflect the characteristics of the fault signal, effectively suppressing phenomena such as modal aliasing and false components. The one-dimensional decomposition component set is converted into a two-dimensional time-frequency diagram using S transform to obtain more refined time-frequency features, present clearer local details, greatly enhance the characteristics of the original signal, and provide a reliable basis for the feature capture of the subsequent classification model. Therefore, the method proposed in the present invention can achieve effective preprocessing of power quality disturbance signals. The method is simple in principle. When processing complex fault current signals, different modal processing results can be obtained, effectively reducing the possibility of aliased modal subsequence escape under a single decomposition result.

Claims

1. A joint multi-scale power quality signal preprocessing method based on singular spectrum analysis and local mean decomposition is characterized by The following steps are involved: Step 1: Decompose the singular spectrum component and the residual component from the disturbance signal based on the singular spectrum decomposition; Step 2: Decompose a series of PF components and residual components from the disturbance signal based on local mean decomposition; Step 3: Merge the two components based on scale matching and weighted superposition to obtain a multi-scale mixed component; Step 4: Analyze the combined multi-scale mixed components based on S transform to obtain the corresponding time-frequency diagram.

2. The method for preprocessing a multi-scale power quality signal based on singular spectrum analysis and local mean decomposition according to claim 1, characterized in that: In the step 1, a two-dimensional signal matrix is ​​first constructed by concatenating a one-dimensional time series with its own lagged version, and the two-dimensional signal matrix is ​​subjected to singular value decomposition to obtain the product of three matrices: a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix; the singular values ​​are arranged in descending order to reflect the contribution of each singular vector to the original matrix information, and the first few largest singular values ​​and their corresponding singular vectors are selected to extract the main trends and periodic components in the time series. These components are also singular spectral components, and each singular spectral component contains part of the signal information.

3. The method for preprocessing a multi-scale power quality signal based on singular spectrum analysis and local mean decomposition according to claim 2 is characterized in that: The step 1 comprises the following steps: Step 1.1: Process the original disturbance signal V(t) into a discrete signal h(n), find the mean z of the initial discrete signal h(n), n = 1, 2, ..., N, N is the length of the sequence after the signal is discretized, and de-mean all discrete sequence point values ​​to obtain a new input signal h′(n); h′(n)=h(n)-z Step 1.2: Calculate the residual component v at the jth iteration of the singular spectrum decomposition SSD j Power spectral density of (n): in, represents Fourier transform, f is the frequency, is the residual component v at the jth iteration j (n), and thus estimate the frequency f corresponding to the maximum peak of the power spectrum max : Among them, arg max f is a function to find the maximum value; if the number of iterations j = 1 and f s is the set sampling frequency, the embedding dimension M is set to N; when the number of iterations j>1 or j=1 and does not satisfy When the embedding dimension M is set to 1.2×(f max / f s ); Step 1.3: Select the embedding dimension M according to the method in step 1.2, and construct the h′(n) obtained in step 1.1 into an M×(N-M+1) signal matrix H. The matrix H is expressed as: In the formula, represents the first sliding window of the signal sequence h'(n), i.e., from h'(1) to h'(1+M-1); Represents the second sliding window of the signal sequence h'(n), i.e., from h'(2) to h'(2+M-1); represents the Mth sliding window of the signal sequence h'(n), i.e., from h'(M) to h'(M+M-1); T represents the sequence transposition; Construct the jth low-frequency subsequence I from high frequency to low frequency j , I j It is expressed as: I j =(i1,i2,…i p ) In the formula, i1,i2,…i p Respectively represent the row index in the trajectory matrix H; The corresponding subsequence is reconstructed by performing diagonal averaging on the matrix. For each diagonal in the matrix H, the average value of its elements is calculated. The diagonal in the matrix H can be expressed as H i,j , where i and j are the row index and column index respectively, and i+j=k, k is a constant representing the index of the diagonal. For each diagonal k, calculate its average value Among them, D k is the set of all elements on the diagonal k, |D k | is the number of elements on the diagonal k; the average of all diagonals Arrange them in the order of the diagonal lines to obtain a reconstructed subsequence of the time series; Step 1.4: Perform singular value decomposition on the matrix H obtained in step 1.3; In the formula, U represents an M×M left singular vector matrix, whose column vectors are the left singular vectors of the matrix H, which are the orthogonal basis vectors in the singular value decomposition of H; V represents a (N-M+1)×(N-M+1) right singular vector matrix; its column vectors are the right singular vectors of the matrix H, which are also the orthogonal basis vectors in the singular value decomposition of H; T represents transpose; ∑ represents an M×(N-M+1) diagonal matrix, whose diagonal elements are the singular values ​​of the matrix H, U=[u1,u2,…u M ]; Where U represents an M×M left singular vector matrix; u1,u2,…u M They represent the nth column vector in the matrix U, that is, the nth left singular vector; V=[v1,v2,…v W ]; Where V represents a (N-M+1)×(N-M+1) right singular vector matrix; v1,v2,…v W They represent the nth column vector in the matrix V, that is, the nth right singular vector; W = N - M + 1; Where W represents the number of right singular vectors, N represents the length of the original sequence, and M represents the embedding dimension; Σ=diag(σ1,σ2,…σ r ); Among them, σ i represents the non-zero singular value of H, i = 1, 2, ..., r, r represents the number of non-zero singular values, and i max =min(M,W), let H i =σ i u i v i T ,u i represents the i-th column vector in the matrix U, that is, the i-th left singular vector; v i T represents the transpose of the i-th column vector in the matrix V; Then the signal matrix can be expressed as: H=H1+H2+…+H q In the formula, H1 represents the first subsequence extracted by singular value decomposition, H2 represents the second subsequence extracted by singular value decomposition, and H q represents the qth subsequence extracted by singular value decomposition; Step 1.5: After each iteration, a new signal component g is separated from the signal j (n) and a residual component v j+1 (n), after each iteration, the residual component v is calculated j+1 The decomposition process is terminated when the calculated mean square error is less than the preset threshold x. Where h(n) represents the original signal; Where mse(j) represents the mean square error between the residual component and the original sequence after the jth iteration, and N represents the sequence length.

4. The method for preprocessing a multi-scale power quality signal based on singular spectrum analysis and local mean decomposition according to claim 1, characterized in that: The step 2 comprises the following steps: Step 2.1: Find the initial signal All extreme points n i , calculate the local mean α of two adjacent extreme points i and the local envelope estimate β i ; Where n i represents the i-th extreme point in the signal, n i+1 Indicates the i+1th extreme point in the signal; Step 2.2: Use the exponentially weighted moving average (EWMA) to adjust the local mean α obtained in step 2.

1. i and the local envelope estimate β i Perform data smoothing to obtain the local mean function α 11 (t) and the local envelope estimation function β 11 (t); a 11 (t)=λa i (t)+(1-λ)a 11 (t-1) b 11 (t)=λβ i (t)+(1-λ)β 11 (t-1) Among them, λ is the smoothing parameter, a i (t) is the local mean at time point t, a 11 (t-1) is the local mean after smoothing at time t-1; the initial signal Subtract the smoothed local mean function α 11 (t) Get the zero mean signal C 11 (t); C 11 (t)=x(t)-α 11 (t) Step 2.3: From the modulated signal, that is, the zero-mean signal C 11 (t) to extract the original information, namely the modulation function d 11 (t); Step 2.4: Repeat steps 2.1 and 2.2 to iterate and obtain α ik (t) and β ik (t), and then calculate the pure modulation function ρ 1k The value of (t); Where, d 1k-1 (t) represents the PF component obtained in k-1 iterations; α 1k (t) represents the local mean function obtained in the kth iteration; C 1k-1 (t) represents the envelope signal obtained in the k-1th iteration; β 1k-1 (t) represents the local envelope estimation function obtained in the k-1th iteration; k represents the number of iterations; n represents the sequence length; Step 2.5: Iterate all local envelope estimation functions β obtained in step 2.4 1k (t) is multiplied to obtain the envelope signal β1(t); In the formula, β 12 (t) represents the local envelope estimation function obtained in the second iteration; β 1n (t) represents the local envelope estimation function obtained in the nth iteration; β 1q (t) represents the local envelope estimation function; q represents the index; Let the envelope signal β1(t) and the first pure modulation function ρ 11 (t) The first product function PF1 can be obtained by multiplying the product; PF1=β1(t)ρ 11 (t) Step 2.6: From the initial signal The PF1 obtained in step 2.5 is stripped off to obtain a new initial signal Repeat steps 2.1 to 2.6 for k times until Shows function monotonicity; In the formula, represents the signal after the first PF component is stripped from the initial signal, that is, the signal after the first iteration, represents the initial signal; PF1(t) represents the first modal component, i.e., the PF component; represents the signal after the second PF component is stripped from the signal after the first iteration; PF2(t) represents the second modal component; represents the signal after stripping the kth PF component from the signal after the k-1th iteration; represents the signal after the k-1th iteration; PF k (t) represents the k-1th modal component; After k cycles, the initial signal can be The decomposition output is k PF components and a residual monotonic component In the formula, PF p (t) represents the pth modal component, p represents the index, p = 1, 2, ..., k.

5. The method for joint multi-scale power quality signal preprocessing based on singular spectrum analysis and local mean decomposition according to claim 1 is characterized in that: The step 3 comprises the following steps: Step 3.1: Analyze the data characteristics of the two components and transform the signal data into the frequency domain through one-dimensional discrete Fourier transform to extract the frequency range of each component: In the formula, X(e jω ) represents the discrete Fourier transform result of the signal; F[x(n)] represents the discrete Fourier transform result of the signal; e -jωn represents a complex exponential function; x(n) is a discrete sequence, that is, two components obtained by SSD decomposition and LMD decomposition, and the frequency components of different components are extracted according to Fourier transform; Step 3.2: By setting the frequency standard, the high frequency, medium frequency and low frequency intervals are divided based on the frequency threshold, and the two groups of components are classified into the set frequency categories. Two thresholds f are set. mid and f low , f mid It is the boundary between mid-frequency and high-frequency, f low It is the boundary between low frequency and medium frequency. After classification, the components of the same frequency category in the two sets are compared for scale matching and adjustment. If a high-frequency component decomposed by SSD is similar to a medium-frequency component decomposed by LMD in frequency, but slightly different, one of the components is filtered or resampled to match its frequency range with the other component. Step 3.3: Determine the sparsification weight coefficient of each component by calculating the signal-to-noise ratio. The signal-to-noise ratio is an indicator that measures the ratio of characteristic information to noise in the signal. The weighted superposition sparsification weight is assigned according to the signal-to-noise ratio. The higher the signal-to-noise ratio of the component, the greater its weight. The calculation formula of the signal-to-noise ratio is as follows: In the formula, SNR represents signal-to-noise ratio, P signal Represents the power of the signal, P noise represents the power of the noise, s(n) represents the original signal, n(m) represents the noise signal, M represents the length of the signal sequence, N represents the length of the noise, m represents the mean of the signal, and n represents the mean of the noise; Step 3.4: Determine the weighted superposition sparsification weight coefficient w by the signal-to-noise of each component calculated in step 3.3 i , and then adjust the sparse weights to ensure that the sum of all weights is 1 to meet the superposition requirements, and use the adjusted weights w i ′For the components x in the same scale in the two groups i (t) is weighted, and then all weighted components are superimposed together to form the final signal y(t). The mathematical expression of this process is as follows: Among them, Normalized SNR i represents the normalized signal-to-noise ratio of the i-th component, τ is the preset signal-to-noise ratio threshold, and the threshold is set to sparse the weights.

6. The method for preprocessing a multi-scale power quality signal based on singular spectrum analysis and local mean decomposition according to claim 1, characterized in that: In step 4, the S transform effectively extracts features using a scalable transformed Gaussian window; the combined signal y(t) is subjected to S transform time-frequency analysis, and the S transform is defined as follows: Where S(τ,f) is the S transform of the signal x(t), t is time, f is frequency, w(τ-t,f) is the Gaussian window, and τ represents the parameter controlling the time axis position of the Gaussian window; The S transform projects the time-frequency domain waveform of the power quality disturbance PQD into a two-dimensional space, and the time-frequency matrix diagram contains the frequency and amplitude information of the disturbance.

Images