一种高准确度非整周期采样谐波分析测量方法

摘要

为了获得高准确的谐波测量结果，本发明提供了一种谐波测量方法，对被测量信号进行采样，对采样结果进行DFT或FFT处理，得到离散频谱，从离散频谱中选取若干谱线修正长范围频谱泄漏、短范围频谱泄漏、负频点泄漏等造成的测量误差，得到高准确度的谐波测量结果。本发明可应用于谐波的精密测量领域。

说明书

技术领域

本发明涉及一种信号测量方法，特别是对包含谐波成分的信号的测量方法。

背景技术

随着电力电子装置、半导体器件等非线性负荷的广泛使用，电力系统中的谐波污染越来越严重。为了保证电能质量以及电力系统的安全运行，应对谐波进行准确检测和分析，从而为治理谐波提供更科学的依据。另外，目前国内外出现了多种谐波测量仪器，为了确保这些仪器量值的统一，也有必要研究高准确度的谐波测量方法。

用于谐波测量的数字信号处理方法有相关分析法、基于离散傅立叶变换(DFT)的快速傅立叶变换(FFT)法、小波变换法、神经网络法及遗传算法等。目前，FFT应用较广。在利用FFT测量含有谐波的信号时，若数据采集系统的采样频率满足采样定理，而且采样样本覆盖的时间是被测信号周期的整数倍，即相当于数据采集系统做到了整周期采样，此时利用FFT可以获得准确度很高的谐波参数。但是在实际中，由于被测信号的周期未知，加上数据采集系统硬件的限制，整周期采样通常很难实现。在非整周期采样条件下，由于长范围泄漏、短范围泄漏、负频点泄漏的影响，利用FFT测量谐波参数时将出现较大的误差。为了提高测量准确度，现有的方法是在忽略负频点泄漏的基础上，通过对采样样本施加窗函数来减小长范围泄漏，采用插值方法来减小短范围泄漏。由于这些方法忽略了负频点泄漏，而采用窗函数也无法将长范围泄漏降至最低，因而其准确度有限。现有的测量方法得到结果的准确度约在10^-5～10^-4之间。

发明内容

为了获得高准确的谐波测量结果，本发明提供了一种谐波测量方法，同时减小长范围泄漏、短范围泄漏和负频点泄漏，从而提高了谐波测量结果的准确度。

本发明的技术方案如下：

高准确度谐波测量方法，包括：

A、对包含谐波的被测量信号进行采样的步骤，和

B、对得到的采用样本进行DFT或FFT处理的步骤，还包括如下步骤：

C、从经过步骤B处理得到的离散频谱中选取谱线，选取的规则是：

选取谱线的数量为K+1，K为大于1的自然数，K的值为基波数量和要测量的谐波数量之和；选取的谱线为第p₁根谱线，第p₂根谱线，......，第p_(K+1)根谱线，其中p₁，p₂，…，p_K+1为谱线在离散频谱中的序号；第p₁根谱线的实部为虚部为第p₂根谱线的实部为虚部为第p₃根谱线的实部为虚部为......，第p_(K+1)根谱线的实部为虚部为

第p₁根谱线选取离散频谱中谱峰最大的谱线，对应于基波分量；

第p₂根谱线选择第p₁根谱线两侧紧邻的两根谱线中幅值较大谱线；

从第p₃根谱线开始直到第p_(K+1)根谱线，每根谱线的选取方法是选取对应各谐波的谱峰；或用谐波次数乘以p₁得到数值s，离散频谱中的第s根谱线作为选取的相应谱线；

D、建立以下矩阵：

通过初等行列变化，使F′变为如下形式

$F=\left(\begin{array}{cc}{F}_1& F_2\\ 0& \left(\begin{array}{cc}{F}_{\mathrm{KK}}& F_{K(K+1)}\\ F_{(K+1)K}& F_{(K+1)(K+1)}\end{array}\right)\end{array}\right)$

求解方程F_KKF_(K+1)(K+1)-F_(K+1)KF_K(K+1)＝0，得到τ；

E、计算b₁，…，b_K

$\left(\begin{array}{c}{b}_1\\ .\\ .\\ .\\ b_K\end{array}\right)={\left(\begin{array}{ccc}\frac{\sin \left(\frac{2\pi p_1}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\pi p_1}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\\ .& & .\\ .& ...& .\\ .& & .\\ \frac{\sin \left(\frac{2\pi p_K}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\pi p_K}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\end{array}\right)}^{-1}·\left(\begin{array}{c}{I}_{p_1}\\ .\\ .\\ .\\ I_{p_K}\end{array}\right)$

和a₁，…，a_K

$\left(\begin{array}{c}{a}_1\\ .\\ .\\ .\\ a_K\end{array}\right)={\left(\begin{array}{ccc}\frac{\sin \left(\frac{2\mathrm{\pi \tau }}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\\ .& & .\\ .& ...& .\\ .& & .\\ \frac{\sin \left(\frac{2\mathrm{\pi \tau }}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\end{array}\right)}^{-1}·\left(\begin{array}{c}{R}_{P_1}-\underset{k=1}{\Sigma }{b}_k\\ .\\ .\\ .\\ R_{P_K}-\underset{k=1}{\Sigma }{b}_k\end{array}\right)$

得到：第k次谐波的频率为第k次谐波的幅值A_k为第k次谐波的初相位为

步骤A所述对被测量信号进行采样为非整周期采样，即采用周期与被测信号周期之间为非整数倍关系。

步骤A所述对被测量信号进行采样还可以是整周期采样，即采用周期与被测信号周期之间为整数倍关系。

本发明的技术效果：

以下通过理论推导的方式对本发明达成的技术效果进行说明。

时域内连续的谐波信号通常可表示为如下形式

其中，t为时间；f为基波信号的频率；k为谐波次数，k＝1时表示基波；A_k为第k次谐波的幅值；为第k次谐波的初相位。

忽略模数转换过程中的量化误差，以及测量过程中的各种随机误差，利用采样频率为f_s(邻近两次采样的时间间隔为T_s＝1/f_s)的数据采集系统得到N个样本

其中n为大于等于0且小于等于N-1的整数。

u(n)的FFT的结果为

$U_m=\frac{1}{N}\underset{n=0}{\overset{N-1}{\Sigma }}u\left(n\right)e^{-j\left(\frac{2\pi }{N}\right)\mathrm{nm}}---\left(3\right)$

其中，m为大于等于0且小于等于N-1的整数。

在整周期采样条件下，即kNf＝pf_s，p为整数，则第k次谐波的频率为pf_s/N，幅值和初相位可利用u(n)的离散频谱中的第p根谱线获得，因为

在非整周期采样条件下，kNf＝(p+ε)f_s，|ε|＜1且ε≠0，此时第p根谱线的值变为

$=\frac{1}{2\mathrm{jN}}[U_D+U_C-U_F]$

其中，τ＝Nf/f_s。由短范围频谱泄漏造成的测量误差存在于U_D中，由长范围频谱泄漏造成的误差存在于U_C中，由负频点频谱泄漏造成的误差存在于U_F中。由于上述泄漏的影响，利用上式很难准确得到第k次谐波的参数。现有的技术或者忽略长范围频谱泄漏和短范围频谱泄漏造成的误差，或者忽略长范围频谱泄漏造成的误差，或者忽略负频点频谱泄漏造成的误差，因此得到的第k次谐波的参数的准确度有限。

在非整周期采样条件下，令根据式(5)可知，第p根谱线的虚部为

$I_p=\underset{k}{\Sigma }\frac{b_k}{2}\frac{2\sin \left(\frac{2\mathrm{\pi p}}{N}\right)}{\cos \left(\frac{2\mathrm{\pi p}}{N}\right)-\cos \left(\frac{2\mathrm{\pi k\tau }}{N}\right)}---\left(6\right)$

第p根谱线的实部为

$R_p=\underset{k}{\Sigma }[b_k+\frac{a_k}{2}\frac{2\sin \left(\frac{2\mathrm{\pi k\tau }}{N}\right)}{\cos \left(\frac{2\mathrm{\pi p}}{N}\right)-\cos \left(\frac{2\mathrm{\pi k\tau }}{N}\right)}]---\left(7\right)$

在利用数据采集系统采集时域连续信号时，数据采集系统的采样频率f_s以及获取的样本数N均已知。对采样数据进行FFT后，第p根谱线的虚部I_p和实部R_p也已知，而需要测量的基波和谐波的数量K也已知。那么，在式(6)中，只剩下b_k、τ是未知数。为了准确测量τ，按前述发明技术方案中所述方法选取(K+1)根谱线来修正长范围频谱泄漏、短范围频谱泄漏、负频点频谱泄漏等造成的测量误差。假设这(K+1)根谱线分别为第p₁，p₂，…，p_K+1根谱线，根据式(6)有

在式(8)中，由于b₁，b₂，…，b_K总计K个未知数，而方程总计(K+1)个，因而通过初等行列变换，便可消去b₁，b₂，…，b_K得到一个关于τ的方程，求解该方程便可求得τ。该过程可利用下述矩阵化简描述。根据式(8)，可构造如下矩阵

对式(9)中的矩阵进行初等行列变换，可得如下形式

$F=\left(\begin{array}{cc}{F}_{11}& F_{12}\\ 0& \left(\begin{array}{cc}{F}_{\mathrm{KK}}& F_{K(K+1)}\\ F_{(K+1)K}& F_{(K+1)(K+1)}\end{array}\right)\end{array}\right)---\left(10\right)$

令$\frac{F_{\mathrm{KK}}}{F_{(K+1)K}}=\frac{F_{K(K+1)}}{F_{(K+1)(K+1)}},$则可得到$g(N,p_1,...,p_K,I_{p_1},...,I_{p_{K+1}},\tau )=0,$求解该方程可得到多个τ解，选择距离p₁最近的解为最终的τ解，再根据τ＝Nf/f_s，便可求出信号的基波频率f＝τf_s/N。

根据式(6)，有

$\left(\begin{array}{ccc}\frac{\sin \left(\frac{2\pi p_1}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\pi p_1}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{\sin \left(\frac{2\pi p_K}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\pi p_K}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\end{array}\right)\left(\begin{array}{c}{b}_1\\ .\\ .\\ .\\ b_K\end{array}\right)\left(\begin{array}{c}{I}_{p_1}\\ .\\ .\\ .\\ I_{p_K}\end{array}\right)---\left(11\right)$

求出b₁，…，b_K后，根据式(7)，有

$\left(\begin{array}{ccc}\frac{\sin \left(\frac{2\mathrm{\pi \tau }}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{\sin \left(\frac{2\mathrm{\pi \tau }}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\end{array}\right)\left(\begin{array}{c}{a}_1\\ .\\ .\\ .\\ a_K\end{array}\right)=\left(\begin{array}{c}{R}_{P_1}-\underset{k}{\Sigma }{b}_k\\ .\\ .\\ .\\ R_{P_K}-\underset{k}{\Sigma }{b}_k\end{array}\right)---\left(12\right)$

由上式可求出a₁，…，a_K。再根据可求出第k次谐波的幅值A_k和初相位分别为

$A_k=\frac{2N\sqrt{a_k^2+b_k^2}}{1-\cos \left(2\mathrm{\pi k\tau }\right)}$

从上述推导过程可见，本发明方法能够同时消除长范围泄漏效应、短范围泄漏效应和负频点泄漏效应，具有较高的准确度。本发明获得的谐波参数的准确度可达10^-7以上。

附图说明

图1为本发明方法的流程图。

图2为给出的测量实例得到的离散频谱。

图2中的标识说明如下：

1、第18根谱线；2、第19根谱线；3、第53根谱线；4、第87根谱线。

具体实施方式

以下结合附图对本发明的技术方案进行详细说明。

如图1所示的流程图，本发明的高准确度谐波测量方法具体步骤如下：

首先需要利用数据采集系统对输入的被测谐波信号(即包含谐波的被测量信号)进行采样，获得N个采样样本，N为自然数，f_s为采样频率。

然后对采集到的N个采样样本进行离散傅立叶变换(DFT)或快速傅立叶变换(FFT)，得到被测量信号的离散频谱。

第三步，从离散频谱中选取K+1根谱线，选取谱线的规则为：

K为大于1的自然数，K的值为基波数量和要测量的谐波数量之和；假定选取的谱线为第p₁根谱线，第p₂根谱线，......，第p_(K+1)根谱线，其中p₁，p₂，…，p_K+1为谱线在离散频谱中的序号；第p₁根谱线的实部为虚部为第p₂根谱线的实部为虚部为第p₃根谱线的实部为虚部为......，第p_(K+1)根谱线的实部为虚部为

第p₁根谱线选取离散频谱中谱峰最大的谱线，第p₁根谱线实际是选取基波谱峰，一般基波的成分都大于谐波的成分，所以基波的谱峰在离散频谱中最高。第p₂根谱线选择第p₁根谱线两侧紧邻的两根谱线中幅值较大谱线。从第p₃根谱线开始直到第p_(K+1)根谱线，每根谱线的选取方法是以下方法之一：

1、选取对应各谐波的谱峰作为相应谱线。即在对应谐波的谱线及邻近该谱线的若干根谱线选择幅值最大的谱线作为相应谱线。

2、用谐波次数乘以p₁得到数值s，离散频谱中的第s跟谱线作为选取的相应谱线。

上述两个选择谱线的方法优选第1种方法，如果利用第1种方法不能选择出谱峰，则采用第2种方法。

记录每根谱线实部的大小(i＝1，…，K+1)和虚部的大小(i＝1，…，K+1)。

第四步，构建如下矩阵

并将其化简成如下形式

$F=\left(\begin{array}{cc}{F}_1& F_2\\ 0& \left(\begin{array}{cc}{F}_{\mathrm{KK}}& F_{K(K+1)}\\ F_{(K+1)K}& F_{(K+1)(K+1)}\end{array}\right)\end{array}\right)$

第五步，求解方程F_KKF_(K+1)(K+1)-F_(K+1)KF_K(K+1)＝0，选择距离p₁最近的解为最终的τ解τ。得到第k次谐波的频率kf为

第六步，计算b₁，…，b_K

$\left(\begin{array}{c}{b}_1\\ .\\ .\\ .\\ b_K\end{array}\right)={\left(\begin{array}{ccc}\frac{\sin \left(\frac{2\pi p_1}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\pi p_1}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\\ .& & .\\ .& ...& .\\ .& & .\\ \frac{\sin \left(\frac{2\pi p_K}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\pi p_K}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\end{array}\right)}^{-1}·\left(\begin{array}{c}{I}_{p_1}\\ .\\ .\\ .\\ I_{p_K}\end{array}\right)$

和a₁，…，a_K

$\left(\begin{array}{c}{a}_1\\ .\\ .\\ .\\ a_K\end{array}\right)={\left(\begin{array}{ccc}\frac{\sin \left(\frac{2\mathrm{\pi \tau }}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}{\cos \left(\frac{2\pi p_1}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\\ .& & .\\ .& ...& .\\ .& & .\\ \frac{\sin \left(\frac{2\mathrm{\pi \tau }}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi \tau }}{N}\right)}& ...& \frac{\sin \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}{\cos \left(\frac{2\pi p_K}{N}\right)-\cos \left(\frac{2\mathrm{\pi K\tau }}{N}\right)}\end{array}\right)}^{-1}·\left(\begin{array}{c}{R}_{P_1}-\underset{k}{\Sigma }{b}_k\\ .\\ .\\ .\\ R_{P_K}-\underset{k}{\Sigma }{b}_k\end{array}\right)$

得到第k次谐波的幅值

得到第k次谐波的初相位

至此，获得了需要测量的谐波信号的参数。本发明的方法可以适用于非整周期采样和整周期采样的情况。

以下应用本发明方法对一个仿真实例进行测算，以验证本发明的方法。

利用MATLAB软件产生一个仿真实例，具体设置为：基波频率为f＝50.5Hz，含有三次谐波分量、五次谐波分量的正弦信号，其表达式为：

y(t)＝5sin(2πft+1)+0.5sin(6πft-0.3)+0.08sin(10πft+1.2)

利用采样频率为f_s＝1500Hz的数据采集系统采得N＝512个样本进行FFT处理，所得到的离散幅频特性如图2所示。根据本发明的方法选取4根谱线(即选择K+1根谱线，K值为谐波数量2与基波数量1的和)，第1根谱线选择离散频谱中谱峰最大的谱线，即图2中第18根谱线(标识为1的谱线)。第2根谱线选择第1根谱线左右两根相邻谱线中幅值较大的一个，即图2中第19根谱线(标识为2的谱线)。第3根谱线选择：三次谐波对应的谱峰为第53根谱线(标识为3的谱线)。第4根谱线选择：五次谐波对应的谱峰为第87根谱线(标识为4的谱线)。

选择谱线的具体数据如下表：

  谱线序号  离散频谱实部  离散频谱虚部  18  4.5188  0.7865  19  -1.3563  -0.2507  53  -0.4171  -0.1637  87  0.0683  0.0105

构造矩阵F′，并将其化简成和F一样的形式，得到关于τ的方程，求解可得τ＝17.237333333333311，于是基波的频率为50.499999999999936。

计算b₁，b₂，b₃可得：b₁＝-0.001152243419121，b₂＝-0.000274816801420，b₃＝-0.000018501183705。

计算a₁，a₂，a₃可得：a₁＝-0.006524189713665，a₂＝-0.000716996787663，a₃＝-0.000084464919214。

最后可得，基波的幅值为4.999999999999942，初相位为1.000000000000070；

三次谐波的幅值为0.500000000000019，初相位为-0.299999999999807；

九次谐波的幅值为0.079999999999997，初相位为1.200000000000319。

由上可见，本发明方法的准确度可达10^-8以上。现有测量方法均达不到该准确度。