阅读说明：本技术 一种电力系统中谐波和间谐波的分析方法与装置 (Method and device for analyzing harmonic waves and inter-harmonic waves in electric power system ) 是由 蔡晓峰 张鸿博 窦小磊 王佳欣 于 2019-11-13 设计创作，主要内容包括：本发明涉及一种电力系统中谐波和间谐波的分析方法与装置,属于电网谐波分析领域,方法包括获取电力系统信号,采用Root-MUSIC算法求出信号的谐波和间谐波频率,根据谐波和间谐波中相邻两个频率之间的频率间隔和设定的间隔阀值对谐波和间谐波进行分组；对于加窗处理后的信号,利用离散傅里叶变换进行谱分析,得到加窗处理后信号的频谱；当分组中谐波和间谐波数量M≥2时,选取与分组中各谐波和间谐波距离最近的M条谱线,通过构建谱线方程组求解谐波和间谐波的幅值和相位。本发明估计出谐波和间谐波的频率后,根据频率间隔对(间)谐波进行分组求解,运算量较小,能有效克服频率相近的间谐波间以及谐波和间谐波之间的主瓣干扰,从而获得较高精度的幅值和相位。(The invention relates to a method and a device for analyzing harmonic waves and inter-harmonic waves in an electric power system, belonging to the field of electric network harmonic wave analysis, wherein the method comprises the steps of obtaining signals of the electric power system, solving the frequencies of the harmonic waves and the inter-harmonic waves of the signals by adopting a Root-MUSIC algorithm, and grouping the harmonic waves and the inter-harmonic waves according to the frequency interval between two adjacent frequencies in the harmonic waves and the inter-harmonic waves and a set interval threshold value; performing spectrum analysis on the windowed signal by using discrete Fourier transform to obtain a frequency spectrum of the windowed signal; and when the number M of the harmonic waves and the inter-harmonic waves in the grouping is more than or equal to 2, selecting M spectral lines closest to the harmonic waves and the inter-harmonic waves in the grouping, and solving the amplitudes and phases of the harmonic waves and the inter-harmonic waves by constructing a spectral line equation set. After the frequencies of the harmonic waves and the inter-harmonic waves are estimated, the (inter) harmonic waves are grouped and solved according to the frequency intervals, the calculation amount is small, and the main lobe interference among the inter-harmonic waves with similar frequencies and among the harmonic waves and the inter-harmonic waves can be effectively overcome, so that the amplitude and the phase with high precision are obtained.)

1. A method for analyzing harmonics and inter-harmonics in an electrical power system, comprising the steps of:

acquiring a signal of a power system, solving the frequency of harmonic waves and inter-harmonic waves of the signal by adopting a Root-MUSIC algorithm, and grouping the harmonic waves and the inter-harmonic waves according to the frequency interval between two adjacent frequencies in the harmonic waves and the inter-harmonic waves and a set interval threshold value;

windowing the signal by adopting a window function, and performing spectrum analysis on the windowed signal by utilizing discrete Fourier transform to obtain a frequency spectrum of the windowed signal;

when the number of the harmonic waves and the inter-harmonic waves in the group is M and M is more than or equal to 2, selecting M spectral lines closest to the harmonic waves and the inter-harmonic waves in the group, and constructing and solving the following spectral line equation set to obtain the amplitudes and phases of the harmonic waves and the inter-harmonic waves;

X＝WA

in the system of equations, the matrix X, W, A is a vector or matrix of complex numbers, where X is a vector of complex values corresponding to the M spectral lines, W is a matrix of spectral function values corresponding to the M spectral lines,

2. The method of analyzing harmonics and inter-harmonics in an electric power system according to claim 1, wherein the amplitude and phase of the harmonics or inter-harmonics having only a single frequency in the group are calculated as follows:

A＝2|X_w(k_1r)|/|W(2π(k_1r-k₁)/N)|

θ＝arg[2jX_w(k_1r)/W(2π(k_1r-k₁)/N)]

wherein A is the amplitude of the harmonic or inter-harmonic to be solved, and X_w(k_1r) As an angular frequency omega from the original signal₁Nearest k (th)_1rValue of the lines, k₁Is an angular frequency omega₁Corresponding line coordinate, W is the kth_1rThe spectrum function value of the spectral line, N is the number of sampling points of the signal after windowing processing by utilizing discrete Fourier transform; theta is the phase of the harmonic or the inter-harmonic to be solved, and arg represents the phase of the complex number to be solved.

3. The method according to claim 1, wherein the set interval threshold is determined according to a main lobe width corresponding to the window function when the window function is a combined cosine window function.

4. The method according to claim 1, wherein when M spectral lines closest to each harmonic and inter-harmonic in the group are the same spectral line K, the remaining M-1 spectral lines are selected from adjacent spectral lines of the spectral line K.

5. An apparatus for processing harmonics and inter-harmonics in an electrical power system, comprising a processor for executing instructions to perform the steps of:

acquiring a signal of a power system, solving the frequency of harmonic waves and inter-harmonic waves of the signal by adopting a Root-MUSIC algorithm, and grouping the harmonic waves and the inter-harmonic waves according to the frequency interval between two adjacent frequencies in the harmonic waves and the inter-harmonic waves and a set interval threshold value;

windowing the signal by adopting a window function, and performing spectrum analysis on the windowed signal by utilizing discrete Fourier transform to obtain a frequency spectrum of the windowed signal;

when the number of the harmonic waves and the inter-harmonic waves in the group is M and M is more than or equal to 2, selecting M spectral lines closest to the harmonic waves and the inter-harmonic waves in the group, and constructing and solving the following spectral line equation set to obtain the amplitudes and phases of the harmonic waves and the inter-harmonic waves;

X＝WA

in the system of equations, the matrix X, W, A is a vector or matrix of complex numbers, where X is a vector of complex values corresponding to the M spectral lines, W is a matrix of spectral function values corresponding to the M spectral lines,

6. The apparatus for processing harmonics and inter-harmonics in an electric power system according to claim 5, wherein the amplitude and phase of the harmonics or inter-harmonics having only a single frequency in the group are calculated as follows:

A＝2|X_w(k_1r)|/|W(2π(k_1r-k₁)/N)|

θ＝arg[2jX_w(k_1r)/W(2π(k_1r-k₁)/N)]

wherein A is the amplitude of the harmonic or inter-harmonic to be solved, and X_w(k_1r) As an angular frequency omega from the original signal₁Nearest k (th)_1rValue of the lines, k₁Is an angular frequency omega₁Corresponding line coordinate, W is the kth_1rThe spectrum function value of the spectral line, N is the number of sampling points of the signal after windowing processing by utilizing discrete Fourier transform; theta is the phase of the harmonic or the inter-harmonic to be solved, and arg represents the phase of the complex number to be solved.

7. The apparatus for processing harmonics and inter-harmonics in an electric power system according to claim 5, wherein when the window function is a combined cosine window function, the set interval threshold is determined according to a main lobe width corresponding to the window function.

8. The apparatus for processing harmonics and interharmonics in an electric power system according to claim 5, wherein when M spectral lines closest to each of harmonics and interharmonics in the group are the same spectral line K, the remaining M-1 spectral lines are selected from adjacent spectral lines of the spectral line K.

Technical Field

The invention belongs to the field of power grid harmonic analysis, and particularly relates to a method and a device for analyzing harmonic waves and inter-harmonic waves in a power system.

Background

With the wide application of new energy power generation and power electronic technology in power systems in recent years, a large amount of harmonics and inter-harmonics appear in a power grid, and meanwhile, various precise power utilization equipment sensitive to the power quality is also popularized continuously, so that the problems of the harmonics and the inter-harmonics are more prominent, accurate power grid (inter) harmonic analysis is a premise of the technologies such as power quality monitoring, harmonic treatment, power metering and the like, and the following introduction is made to several common methods for harmonic and inter-harmonic analysis in the prior art:

FFT (Fast Fourier Transform) is a common method for (inter) harmonic analysis at present, but the problems of spectral leakage and fence effect occur when sampling is not synchronized. The windowing interpolation FFT algorithm can inhibit spectrum leakage and barrier effect caused by asynchronous sampling to a certain extent, but when the frequency of harmonic waves is close to that of inter-harmonic waves or inter-harmonic waves with different frequencies, the serious spectrum interference influences the accurate estimation of the windowing interpolation FFT algorithm on the parameters of the harmonic waves and the inter-harmonic waves. Although the excellent window function can be used for suppressing the spectrum leakage, the main lobe of the window function with excellent performance is wide, and inter-harmonics with similar frequencies or main lobes between harmonics and inter-harmonics are easily overlapped, so that the interference between the main lobes cannot be avoided by analyzing the inter-harmonics in a conventional interpolation mode.

The Multiple Signal Classification (MUSIC) is a modern spectral estimation method for estimating Signal space parameters, and its basic idea is to perform characteristic decomposition on an autocorrelation matrix formed by sampled data to obtain a Signal subspace and a noise subspace, and then estimate the frequency of a Signal by utilizing the orthogonality of the Signal subspace and the noise subspace, and it has high resolution, estimation accuracy and stability. The Root-MUSIC algorithm is an improvement on MUSIC, and replaces the spectral peak search on a frequency domain axis with the polynomial Root calculation, so that the calculation amount is reduced, and meanwhile, the precision is improved. However, both the MUSIC and Root-MUSIC algorithms can only detect the frequency of the intermediate harmonic of the signal, and cannot estimate the amplitude and phase of the intermediate harmonic. The authors in the publication of "power automation equipment" in journal "published a paper" inter-harmonic parameter estimation based on real-valued Root-MUSIC and Prony algorithm ", which firstly uses Root-MUSIC algorithm to estimate the frequency of inter-harmonic harmonics, and then uses Prony algorithm to estimate the amplitude and phase angle, but Prony algorithm is susceptible to noise, and can only solve all inter-harmonic harmonics together, but cannot selectively solve specific harmonics or inter-harmonic signals, resulting in more inter-harmonic harmonics and large algorithm computation when the number of sampling points is large.

Disclosure of Invention

The invention aims to provide a method and a device for analyzing harmonic waves and inter-harmonic waves in a power system, which are used for solving the problems of low precision and large calculation amount of harmonic wave and inter-harmonic wave parameter estimation in the prior art.

Based on the above purpose, a technical scheme of an analysis method for harmonic waves and inter-harmonic waves in a power system is as follows:

acquiring a signal of a power system, solving the frequency of harmonic waves and inter-harmonic waves of the signal by adopting a Root-MUSIC algorithm, and grouping the harmonic waves and the inter-harmonic waves according to the frequency interval between two adjacent frequencies in the harmonic waves and the inter-harmonic waves and a set interval threshold value;

windowing the signal by adopting a window function, and performing spectrum analysis on the windowed signal by utilizing discrete Fourier transform to obtain a frequency spectrum of the windowed signal;

when the number of the harmonic waves and the inter-harmonic waves in the group is M and M is more than or equal to 2, selecting M spectral lines closest to the harmonic waves and the inter-harmonic waves in the group, and constructing and solving the following spectral line equation set to obtain the amplitudes and phases of the harmonic waves and the inter-harmonic waves;

X＝WA

in the system of equations, the matrix X, W, A is a vector or matrix of complex numbers, where X is a vector of complex values corresponding to the M spectral lines, W is a matrix of spectral function values corresponding to the M spectral lines,

m is 1,2, …, M, where a and θ are the amplitude and phase of the harmonics to be found and the inter-harmonics in the packet, respectively.

Based on the above purpose, a technical solution of a processing device for harmonic and inter-harmonic in an electric power system is as follows:

the system comprises a processor and a harmonic analysis module, wherein the processor is used for executing instructions to realize the analysis method of the harmonic and the inter-harmonic in the power system.

The two technical schemes have the beneficial effects that:

according to the method, the frequencies of the harmonic waves and the inter-harmonic waves are estimated by using the Root-MUSIC algorithm, then the frequency distribution density of the harmonic waves and the inter-harmonic waves is utilized, the (inter) harmonic waves are solved in groups according to the frequency interval between two adjacent frequencies in the harmonic waves and the inter-harmonic waves and a set interval threshold value, and the order of a matrix W in each group is very small, so that the amplitude and the phase of the harmonic waves and the inter-harmonic waves can be quickly obtained. Compared with the prior art for solving large equation sets, the matrix order after each group is grouped is far smaller than the order of the matrix during solving in the prior art, and the calculation amount is smaller; and the main lobe interference between inter-harmonics with similar frequencies and between harmonics and inter-harmonics can be effectively overcome, so that amplitude and phase with higher precision are obtained. In addition, compared with the prior art, the method combines the Root-MUSIC algorithm with the method for constructing the spectral line equation set, and has better anti-noise performance on the estimation of the amplitudes of the harmonic waves and the inter-harmonic waves.

For the group with only one harmonic or only one inter-harmonic, because the windowing processing is carried out by adopting the window function, the frequency spectrum leakage generated by the long-distance frequency component is greatly inhibited, and therefore the following amplitude and phase calculation formulas of the harmonic and the inter-harmonic are directly adopted:

A＝2X_w(k_1r)|/|W(2π(k_1r-k₁)/N)

θ＝arg[2jX_w(k_1r)/W(2π(k_1r-k₁)/N)]

wherein A is the amplitude of the harmonic or inter-harmonic to be solved, and X_w(k_1r) As an angular frequency omega from the original signal₁The most recentk_1rValue of the lines, k₁Is an angular frequency omega₁Corresponding line coordinate, W is the kth_1rThe spectrum function value of the spectral line, N is the number of sampling points of the signal after windowing processing by utilizing discrete Fourier transform; theta is the phase of the harmonic or the inter-harmonic to be solved, and arg represents the phase of the complex number to be solved.

When the window function adopts a combined cosine window function, such as a nuttalll window, the set interval threshold is determined according to the width of the main lobe corresponding to the window function.

Further, if the frequencies of the harmonics and/or inter-harmonics within a group are similar, then the M spectral lines closest to each harmonic and inter-harmonic in the group may be the same. Based on this, when M spectral lines closest to each harmonic and inter-harmonic in the group are the same spectral line K, the remaining M-1 spectral lines are selected from adjacent spectral lines of the spectral lines K.

Drawings

FIG. 1 is a schematic illustration of the grouping of harmonics and inter-harmonics in an embodiment of the method of the present invention;

FIG. 2 is a magnitude spectrum of a fast Fourier transform of a signal in an embodiment of the method of the present invention;

FIG. 3 is a flow chart of a method of harmonic and inter-harmonic analysis in a method embodiment of the present invention.

Detailed Description

The following further describes embodiments of the present invention with reference to the drawings.

The method comprises the following steps:

the embodiment provides a method for analyzing harmonics and inter-harmonics in an electric power system, as shown in fig. 3, including the following steps:

firstly, a Root-multiple signal classification algorithm is used for estimating the (inter) harmonic frequencies of the signals, and the (inter) harmonics are grouped according to the frequency intervals among the (inter) harmonic frequencies and a set interval threshold.

Specifically, let harmonic and inter-harmonic sampling signals (power signals) of the power system be expressed as:

wherein x (n) is a sampling signal; n is the number of sampling points; t is_sIs a sampling period; p is the number of harmonics and inter-harmonics; f. of_iThe frequency of the ith harmonic or inter-harmonic component; a. the_i、θ_iThe amplitude and initial phase of the ith harmonic or inter-harmonic component, respectively; v (n) is zero mean and variance σ²White gaussian noise.

Intercepting a signal sampling sequence with the length of N, and constructing a data array X:

firstly, the covariance matrix R of the matrix X is solved_xAnd to the matrix R_xAnd (3) decomposing the characteristic values to obtain m characteristic values which are arranged from small to large as: lambda [ alpha ]₁≈λ₂≈…≈λ_m-2p<λ_m-2p+1<…<λ_m. Wherein λ₁…λ_m-2pThe m-2p eigenvalues are significantly smaller than lambda_m-2p+1…λ_mThe 2p characteristic values can judge the number of the harmonic waves and the inter-harmonic waves in the signal to be p. The eigenvectors corresponding to the larger 2p eigenvalues are expanded into a signal subspace of U_SThe expansion of the eigenvectors corresponding to the smaller m-2p eigenvalues into the noise subspace is U_N。

In order to extract spatial parameter information using the noise feature vector, a Root-MUSIC type polynomial is constructed:

wherein f (z) is a polynomial, and p (z) ═ 1, z, …, z^m-1]^TZ is a complex variable, superscript H denotes conjugate transpose, U_NAnd (4) expanding a noise subspace for the eigenvectors corresponding to the m-2p eigenvalues. This equation yields a zero point corresponding to the frequency of the signal, but the equation is not completely a polynomial of z, and a power term of the conjugate of z exists, so that equation f (z) is 0The root process becomes complex, based on which the above equation can be modified:

the formula is a polynomial of a Root-MUSIC algorithm, the highest degree of the polynomial is 2(m-1), namely, (m-1) pairs of roots exist, p pairs of roots in the (m-1) pairs of roots are just distributed on a unit circle, the p pairs of roots correspond to p frequencies of collected signals, but in practical application, due to the existence of errors in various aspects, the p pairs of roots are not necessarily just positioned on the unit circle, and at the moment, only the p pairs of roots closest to the unit circle need to be obtained. Let z be the root closest to the unit circle, the frequency of the corresponding sampled signal (which is a sinusoidal signal) is f-arg (z/(2 π T)_s) Arg denotes the phase of the complex number, T_sRepresenting a sampling time interval.

According to the Euler formula, a real sine can be represented as the sum or difference of two complex sine, so that the Root-MUSIC algorithm is applied to obtain the frequency of the sampling signal in a positive-negative paired form, and the positive frequency is taken as the frequency estimation value of the harmonic wave and the interharmonic wave.

The Root-MUSIC algorithm can detect each frequency of harmonic and inter-harmonic in a signal, and in order to calculate the amplitudes and phases of the harmonic and inter-harmonic, the harmonic and inter-harmonic are grouped, namely, a boundary line is drawn between two adjacent frequency components with the frequency interval above an interval threshold value 6 delta f, and the harmonic and inter-harmonic are divided into a plurality of groups.

The reason that 6 Δ f (Δ f is the spectral line interval corresponding to the discrete fourier transform, i.e., the frequency resolution) is adopted as the division width (i.e., the set interval threshold) is that 6 Δ f corresponds to the main lobe width of the cosine window function of 6 terms added to the windowed FFT in the following step two, if two frequency components are separated by more than 6 Δ f, it is indicated that the two frequency components are both outside the main lobe of the window function of the other party, and because of the excellent performance of the cosine window of 6 terms, the leakage spectrum outside the main lobe is greatly suppressed, so the leakage interference between different groups can be ignored. For example, assume a signal contains f₁-f₆6 frequency components of (1), shown in FIG. 1According to the frequency interval relationship in the figure, f₁-f₃Are divided into a group f₄-f₅Are divided into a group f₆Are grouped.

Secondly, the collected signals are subjected to spectrum analysis by utilizing the windowed FFT. The method comprises the steps of windowing signals by adopting a window function, and performing spectrum analysis on the windowed signals by utilizing discrete Fourier transform to obtain frequency spectrums of the windowed signals.

In particular, for a frequency f₁A single frequency signal x (t) with amplitude A and initial phase theta, and a sampling rate f_sObtaining a discrete signal x (n) after the analog-to-digital conversion:

x(n)＝Asin(2πf₁n/f_s+θ)＝Asin(ω₁n+θ) (5)

wherein, ω is₁＝2πf₁/f_sThe digital angular frequency corresponding to the signal x (n).

Defining the time domain form of a combined cosine window with the length of N as W (N), and recording the signal subjected to windowing by adopting a window function as x if the continuous frequency spectrum of the combined cosine window is W (2 pi f)_w(n) if x is ignored_w(n) negative frequency point-omega₁Side lobe effect of the processing peak, x_w(n) at positive frequency ω₁Nearby continuous spectral function X_w(ω) can be expressed as:

wherein j is a complex symbol; the above formula is subjected to N-point discrete sampling to obtain X_wExpression X of N-point discrete Fourier transform of (omega)_w(k Δ ω) is:

where Δ ω is a discrete frequency interval, N is a data truncation length, the digital frequency interval Δ ω is 2 pi/N, and the corresponding analog frequency interval is Δ f_sK is spectral line coordinate, k is 0,1,2, …, N-1, k₁Is frequency omega₁Correspond toLine coordinate of (a), k₁＝ω₁/Δω＝f₁/Δf＝Nf₁/f_s。

Due to frequency ω of x (n)₁It is difficult to locate exactly at the frequency point of the discrete spectral line, and therefore k₁Are generally not integers; suppose distance ω₁The nearest spectral line is the kth_1rStrip, then there are:

Ae^jθ＝2jX_w(k_1r)/W(2π(k_1r-k₁)/N) (9)

wherein A is the amplitude of the harmonic or inter-harmonic to be solved, and X_w(k_1r) As an angular frequency omega from the original signal₁Nearest k (th)_1rValue of the lines, k₁Is an angular frequency omega₁Corresponding line coordinate, W is the kth_1rA spectral function of the spectral lines, wherein N is the number of sampling points for performing discrete Fourier transform on the windowed signal; theta is the phase of the harmonic or inter-harmonic to be solved.

Combining the above two formulas, the signal amplitude a is obtained as:

A＝2|X_w(k_1r)|/|W(2π(k_1r-k₁)/N)| (10)

the initial phase θ of the signal is:

θ＝arg[2jX_w(k_1r)/W(2π(k_1r-k₁)/N)](11)

wherein, | | represents a modulus of the complex number, and arg represents a phase of the complex number.

In order to better suppress the spectrum leakage under asynchronous sampling, before performing discrete fourier transform, a windowing process is performed on an original signal in a time domain by using a window function, in this embodiment, a combined cosine window is used, and a calculation formula of an N-point combined cosine window is as follows:

where i is 0,1,2,.., N-1, D is the number of terms of the window function, a_dAre coefficients of a window function.

The N-point combined cosine window spectrum W (ω) can be expressed as:

wherein:

for ease of analysis, let ω ═ 2 π δ/N simplify to obtain:

according to formula (6) and formula (15), let δ be k-k₁The discrete spectral lines obtained by spectral analysis of the available windowed FFT can be expressed as:

the cosine window adopted in this embodiment is a 6-term cosine window, and the coefficients of the corresponding window functions are respectively:

a₀＝0.24609375，a₁＝0.41015625，a₂＝0.234375，a₃＝0.087890625，a₄＝0.01953125，a₅＝0.001953125。

table 1 contrasts and analyzes the side lobe characteristic parameters of the 6 cosine windows and the common window function, and it can be seen that the 6 cosine windows have good performances of small side lobe peak value and high attenuation speed.

TABLE 1

And thirdly, for dense grouping of harmonic waves and inter-harmonic waves, obtaining harmonic wave and inter-harmonic wave parameters, namely amplitude and phase, by constructing a spectral line equation set and solving. The specific process is as follows:

when a packet contains multiple frequency components, the result of spectral analysis using a windowed FFT is a superposition of the spectra of the individual frequency components, for example x (n) is a superposition of 3 frequency components:

x(n)＝A₁sin(ω₁n+θ₁)+A₂sin(ω₂n+θ₂)+A₃sin(ω₃n+θ₃) (17)

then according to equation (7) there is:

in the formula k₁＝ω₁/Δω,k₂＝ω₂/Δω,k₃＝ω₃A,/Δ ω, find h₁、h₂、h₃The three spectral lines combine the upper formula (18) to write three equations in a column form to form an equation set, and the equation set is written in a matrix form:

in the formula, h₁,h₂,h₃Representing three selected spectral lines, and the position coordinates of the three spectral lines are h in sequence₁、h₂、h₃，X(h₁),X(h₂),X(h₃) Corresponding to the values of three spectral lines, A₁,A₂,A₃For three amplitudes to be solved in the group, theta₁,θ₂,θ₃For the three phases to be sought in the packet.

In the formula (19), W is a matrix, and each element in the matrix is a frequency spectrum formed by combining N points with a cosine window

That is, the above equation (15) is calculated, but the δ value corresponding to each element is different.Equation (19) is abbreviated: x is WA, where X is 2j X [ X (k)₁) X(k₂) X(k₃)]^T。

In the above formula, the elements in the vector X can be obtained by windowed FFT, and the matrix W can also be calculated because the Root-MUSIC algorithm has already obtained the frequency values contained in the signal, so that only the vector a is unknown, and a ═ W can be obtained^-1X，W^-1Is an inverse matrix of W. The modulus of the element in A corresponds to the amplitude of each frequency component, and the phase angle corresponds to the initial phase angle of each frequency component.

In the above description, the x (n) includes 3 frequency components, and when x (n) includes any M frequency components (M is an integer of 2 or more), the following steps are performed:

x(n)＝A₁sin(ω₁n+θ₁)+A₂sin(ω₂n+θ₂)…A_Msin(ω_Mn+θ_M) (20)

then according to equation (7) there is:

in the formula k₁＝ω₁/Δω,k₂＝ω₂/Δω,…,，k_M＝ω_MA,/Δ ω, find h₁、h₂、…,,h_MThe total M spectral lines are combined with the upper formula (20), M equations can be written in a column mode to form an equation system, and the equation system is written in a matrix form:

in the formula, h₁,h₂,…,h_MRepresenting the selected M spectral lines, and the position coordinates of the M spectral lines are h in sequence₁、h₂、…,h_M，X(h₁),X(h₂),…,X(h_M) Corresponding to the values of M spectral lines, A₁,A₂,…,A_MFor M amplitudes, θ, to be found in the packet₁,θ₂,…,θ_MFor the M phases to be sought in the packet.

It should be noted that the grouping in this embodiment is a grouping of harmonics and/or inter-harmonics of at least two frequency components, and the grouping may have only harmonics or only inter-harmonics, or both harmonics and inter-harmonics, so the spectral order in each grouping is the same as the number of harmonics and inter-harmonics grouped, and is at least two (i.e., a 2 × 2 matrix).

In this embodiment, for the group of harmonics and inter-harmonics having only one frequency component, since windowing processing is performed by using a window function with excellent performance, the spectrum leakage generated by the remote frequency component is greatly suppressed, so that the harmonic and inter-harmonic parameters can be directly calculated according to the equations (10) and (11), and for the group containing two or more frequency components, the above method for constructing the spectral line equation set (i.e., the method in step three) is used for solving.

In this embodiment, in selecting the spectral lines, the spectral line closest to each (inter) harmonic frequency in the group is selected first, because the spectral line amplitude near the signal frequency point after the fourier transform of the time-domain signal is large, and the noise immunity is strong. If the inter-harmonics are close in frequency and the lines closest to them may be the same line, e.g. line K, then the remaining lines are selected from the neighbouring lines of line K. In this way, the solution of the equation system is not needed for each frequency component, but only when needed, and the solution of the large equation system is avoided by adopting the grouping solution, so that the computation amount is reduced compared with the algorithm (such as the Prony algorithm) in the prior art. The Prony algorithm can only solve all harmonics and inter-harmonics together, the number of the harmonics and inter-harmonics to be analyzed is large, the number of sampling points is large, the algorithm computation is large, and the method cannot selectively solve the harmonics and the inter-harmonics in groups.

In addition, the equation set in the embodiment is constructed and solved in a frequency domain, and because the 6-term cosine window is adopted for processing, the 6-term cosine window has excellent spectrum leakage suppression capability, the spectrum leakage of other frequency components with long distance is greatly suppressed, and conditions are provided for obtaining a high-precision inter-harmonic estimation result.

Compared with the operation quantity of the Root-MUSIC + Prony algorithm in the prior art, the two algorithms need to solve the frequency of harmonic waves between harmonics through the Root-MUSIC algorithm, so the operation quantity of the two algorithms is the same in frequency estimation, but the operation quantity of the two algorithms is greatly different in amplitude and phase estimation, the Prony algorithm needs to construct a Van der Menu matrix according to a frequency estimation result, and then the solution is carried out according to a least square method, namely:

a＝(Z^HZ)^-1Z^Hx

wherein, a is a vector formed by the amplitude and the phase of the signal to be solved, x is a vector formed by sampling points, Z is a vandermonde matrix, H represents conjugate transpose, the order p thereof is equal to the number of inter-harmonic harmonics in the signal, and as 1 real sine can be represented as the sum or difference of 2 complex sinusoids, the order p of Z is 2 times of the number of inter-harmonic harmonics in the signal for the real signal. . According to the power grid harmonic measurement standard, 25 harmonics generally need to be measured, and in addition, the requirement of measuring inter-harmonics, the order p of Z is not suitable to be smaller than 50, so that the inversion operation of a 50-order matrix is involved, and in addition, multiple matrix multiplication operations are involved, specifically:

Z^Hz is the multiplication of the p × N and N × p matrices (a p × p matrix is obtained, each element of the matrix needs N multiplications for a total of p × p × N multiplications), Z^HInverse of Z and then multiplying by Z^HIs multiplication of a p x p matrix and a p x N matrix (a p x N matrix is obtained, where each element of the matrix requires p multiplications for a total of p x N x p multiplications), and finally (Z^HZ)^-1Z^HAnd the vector x multiplication is the multiplication of a p multiplied by N matrix and an N multiplied by 1 vector (p multiplied by N), N is the number of sampling points, and finally, the amplitude and the phase of harmonic waves among harmonics are estimated by a Prony algorithm, wherein the sum of the 2 multiplied by p multiplied by N + p multiplied by N and the p multiplied by p matrix inversion are required.

The method of the present invention requires a fast Fourier transform operation of N/2 × log₂N multiplications, assuming 102 samples (51 harmonics can be analyzed theoretically) per cycle and 1024 samples in total, the fast Fourier transform only needs 5 XN multiplications, which is much smaller than 2 Xp × p × N + p × N (p × N) of Prony method>50) And (4) secondary multiplication. In addition, the method of the invention carries out grouping solution on the harmonic waves and the inter-harmonic waves according to the frequency intervals, the scale of the constructed spectral line equation set is generally smaller,the solution operation amount is far less than p multiplied by p (p)>50) And (5) calculating the inverse of the order matrix. Therefore, the method of the present invention is far less computationally intensive than the Prony algorithm in estimating the amplitudes and phases of harmonics and inter-harmonics.

The following simulations were performed according to the method of the invention as follows:

the simulation signal comprises harmonic components and inter-harmonic components, and the expression is as follows:

wherein p is the number of harmonics, A_iIs the i-th harmonic amplitude, f₁At fundamental frequency, θ_iIs the initial phase of the i-th harmonic; q is the number of inter-harmonics, f_jFrequency of interharmonic j, A_jAmplitude of the inter-harmonics number j, θ_jIs the initial phase of the j-th inter-harmonic.

According to the requirement of the power grid harmonic measurement standard, 40-th harmonic needs to be measured, but the 40-th harmonic has too many parameters, limited by space, only some low-order harmonics and inter-harmonics are considered below, and the parameters of each frequency component are shown in the following table:

TABLE 2

The sampling frequency is set to 5120Hz, and the total sampling frequency is 1024 points for 10 rated power frequency periods (0.2 s).

The amplitude spectrum of the fast fourier transform of the signal is shown in fig. 2, since 53Hz inter-harmonic is very close to the fundamental wave, and the amplitude is much smaller than the fundamental wave, the two main lobes are overlapped, and the 53Hz inter-harmonic is seriously interfered by the fundamental wave main lobe, which cannot be observed from the amplitude spectrum. Meanwhile, due to the existence of 53Hz inter-harmonics, the measurement of the fundamental frequency is greatly influenced, and six items are based according to the article in the journal of protection and control of electric power system by the authorThe fundamental frequency of the calculation is 50.2562Hz and the error reaches 0.0562Hz according to the windowed interpolation principle provided by the cosine window four-spectral line interpolation FFT high-precision harmonic detection algorithm, because the frequency is normally i multiplied by f in the harmonic analysis of the fast Fourier transform₀Since the frequency of the higher harmonic wave is a deviation of the frequency of the fundamental wave, it is not clear that the amplitude and phase correction of the higher harmonic wave are more affected, and therefore, it is not preferable to analyze the harmonic wave and the inter-harmonic wave by the interpolation FFT algorithm for a signal having inter-harmonic wave, particularly, inter-harmonic wave in the vicinity of the fundamental wave. On the other hand, the 66Hz inter-harmonic wave appears on the amplitude spectrum, but is interfered by the fundamental main lobe, and the interpolation FFT algorithm cannot be used to obtain a high-precision parameter estimation result.

The frequency estimation result of Root-MUSIC algorithm in this embodiment is shown in table 3, where the error in the table is relative error, aE ± n represents a × 10^±n(the same below), as can be seen from table 3, the Root-MUSIC algorithm has high frequency resolution, can quite accurately analyze the frequency of inter-harmonic waves, and lays a foundation for analyzing the amplitude and phase of the inter-harmonic waves; however, the array element number M directly affects the accuracy of the Root-MUSIC algorithm, and a larger M value can bring higher accuracy, but the required operation time is long. Simulation shows that when the sampling frequency is 5120Hz and the sampling point is 1024 points, the number M is 400, so that the high operation speed can be obtained while the accuracy is ensured, and the anti-noise capability is good, therefore, the number M of the array elements in the Root-MUSIC algorithm is 400.

The frequency estimation result of the Root-MUSIC algorithm under the noise-free condition is shown in the following table:

TABLE 3

In order to compare and analyze the amplitude and phase measurement precision of harmonic inter-harmonics, a Root-MUSIC + Prony algorithm (called algorithm 1) and a method (called algorithm 2) of the invention, which are provided by a new article in journal electric power automation equipment (called inter-harmonic parameter estimation based on real value Root-MUSIC and Prony algorithms), are adopted to respectively perform simulation analysis of harmonics and inter-harmonics, and it needs to be noted that the article simulates and compares several common inter-harmonic analysis algorithms, verifies that the Root-MUSIC + Prony algorithm has higher precision, so that other methods are not considered, the simulation result is shown in table 4, and the amplitude and phase estimation result under the noise-free condition is shown in the following table:

TABLE 4

Table 4 shows that, under the noise-free condition, the estimation of the amplitude and the phase by the algorithm 1 and the algorithm 2 achieves very high precision, and the amplitude measurement precision of the algorithm 2 for most harmonics and inter-harmonics is higher, although the amplitude measurement precision of the algorithm 2 for 2 nd and 4 th harmonics is slightly lower than that of the algorithm 1, the difference is not large; in the phase measurement precision, the precision of the algorithm 1 is superior to that of the algorithm 2, so that the two algorithms have components with higher precision than the other one and components with lower precision than the other one, and the advantages of the two algorithms are not obvious.

The actually acquired power signals often contain noise, and the results of signal simulation under the condition of meaningful noise are shown in tables 5 and 6, wherein the table 5 is the simulation results of the algorithm 1 and the algorithm 2 when the signal-to-noise ratio is 60dB, and the table 6 is the simulation results of the algorithm 1 and the algorithm 2 when the signal-to-noise ratio is 40 dB.

TABLE 5

TABLE 6

Tables 5 and 6 show that the accuracy of both algorithms is reduced under the condition of noise, and the accuracy of the algorithms is reduced along with the reduction of the signal to noise ratio (because the measurement errors have certain randomness under the noise condition, the errors in the tables are the average value of 20 measurements). Further analyzing tables 5 and 6, it can be easily found that algorithm 2 has higher measurement accuracy for harmonic and inter-harmonic amplitudes than algorithm 1 under the conditions of signal-to-noise ratios of 60dB and 40dB, and shows better anti-noise performance. However, from the comparison of the phase errors of the two algorithms in tables 5 and 6 (the term indicated by bold font in tables 5 and 6 is that the precision of algorithm 1 is better than that of algorithm 2, and the term indicated by non-bold font is that the precision of algorithm 2 is better than that of algorithm 1), the two algorithms respectively have components with precision higher than that of the other side and components with precision lower than that of the other side, and the advantages of the two algorithms are not obvious. However, in the analysis of power harmonics and inter-harmonics, the amplitude of the harmonics and inter-harmonics is often more concerned, and in these applications, the algorithm 2 is favorable for obtaining more accurate amplitude measurement results.

In summary, the method of the present invention has the following advantages:

1. according to the harmonic and inter-harmonic analysis method based on the Root-MUSIC algorithm and the spectral line equation set, the frequencies of the harmonic and inter-harmonic are estimated by the Root-MUSIC algorithm, and then the harmonic and inter-harmonic are solved in groups according to the frequency interval between the harmonic and inter-harmonic, so that the problem of solving a large equation set is avoided, and the calculation amount is small.

2. The spectral line equation set in the embodiment is established on the basis of adding 6-term cosine window FFT, and because the 6-term cosine window has excellent spectrum leakage inhibition capability, the spectrum leakage of other frequency components with longer distance is greatly inhibited, and conditions are provided for obtaining high-precision harmonic and inter-harmonic estimation results.

3. For the grouping with dense inter-harmonics, the method obtains the inter-harmonic parameters by constructing a spectral line equation set and solving, can effectively overcome the main lobe interference between the inter-harmonics with similar frequencies, and obtains higher analysis precision.

4. Compared with the existing Root-MUSIC + Prony algorithm inter-harmonic analysis method, the harmonic and inter-harmonic analysis method provided by the invention has better anti-noise performance on estimation of harmonic and inter-harmonic amplitudes.

The embodiment of the device is as follows:

the embodiment provides a processing device for harmonic waves and inter-harmonic waves in a power system, which comprises a processor and is used for executing instructions to realize the following steps:

acquiring signals of a power system, solving harmonic and inter-harmonic frequencies of the signals by adopting a Root-MUSIC algorithm, and grouping the harmonic and inter-harmonic frequencies according to a frequency interval set between two adjacent frequencies in the harmonic and inter-harmonic frequencies;

windowing is carried out on the signals by adopting a window function, and spectrum analysis is carried out on the signals subjected to windowing by utilizing discrete Fourier transform to obtain frequency spectrums of the signals subjected to windowing;

when the number of the harmonic waves and the inter-harmonic waves in the group is M and M is more than or equal to 2, selecting M spectral lines closest to the harmonic waves and the inter-harmonic waves in the group, and obtaining the amplitudes and phases of the harmonic waves and the inter-harmonic waves by constructing and solving the following spectral line equation set;

X＝WA

in the equation set, the matrix X, W, A is a complex matrix, where X is the spectral function value of the selected M spectral lines, W is the matrix formed by the spectral function values of the corresponding M spectral lines,m is 1,2, …, M, where a and θ are the amplitude and phase of the harmonics to be found and the inter-harmonics in the packet, respectively.

The processing device in the above embodiment is actually a computer solution based on the method flow of the present invention, that is, a software framework, and the above device is a processing procedure corresponding to the method flow. The above-described method is not described in detail, since it is sufficiently clear and complete that it has been described in method embodiments.

In addition, the processor in this embodiment may be a computer, a microprocessor such as an ARM, or a programmable chip such as an FPGA, a DSP, or the like.

The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. For example, in the present embodiment, a window function with a more rapid side lobe attenuation may be used, and accordingly, in the inter-harmonic grouping, the frequency interval required for the grouping may be adjusted according to the main lobe width corresponding to the window function.

Therefore, any modification, equivalent replacement, improvement or the like made within the spirit and principle of the present invention should be included in the scope of the claims of the present invention.