CN105510707A - 一种电力系统谐波和间谐波测量方法 - Google Patents

Abstract

本发明提出了一种电力系统谐波和间谐波测量方法。本方法的主要实现步骤是：首先应用谱估计初步测量电力系统波形中所含谐波和间谐波的数量及频率；然后基于测得的谐波和间谐波的频率构造出混沌检测振子，对其幅值进行测量，为提高检测精度，应用最大李雅普诺夫指数作为混沌检测振子处于混沌状态和大尺度周期状态相变的临界条件的判断依据。本发明能够以较高的频率分辨率和幅值测量精度高效地测出电力系统谐波和间谐波的频率和幅值，且对背景噪声的抗干扰能力强。

Description

一种电力系统谐波和间谐波测量方法

技术领域

本发明涉及电力系统谐波分析领域，更具体地，涉及一种电力系统谐波和间谐波测量方法。

背景技术

目前，谐波测量的主要方法有：快速傅立叶变换、小波变换、瞬时无功功率理论等。快速傅立叶变换是当前应用最广泛的一种谐波测量方法，但该方法存在频谱混叠、栅栏效应及频谱泄露等问题，使得测量结果在一定程度上无法满足电力系统的要求；虽然通过加窗插值修正算法可以较好地提高测量精确度，减少栅栏效应带来的误差，但往往算法比较复杂，编程实现比较繁琐，且实时性差。近年来提出的小波变换方法对信号具有自适应性及良好的时频局部化特征，在谐波测量上进行了有益的尝试，但小波变换是线性变换，也存在频带混叠、频谱泄露，对脉冲干扰的抑制作用不够理想及暂态电能质量信号特征随尺度增加逐渐被削弱等问题。基于瞬时无功功率理论的谐波检测方法虽具有较好的实时性且可全部采用模拟电路实现，但该方法仅适用于三相电压波形对称且无畸变、电流不含零序分量的场合。

发明内容

本发明提供一种电力系统谐波和间谐波测量方法，该方法能够以较高的频率分辨率和幅值测量精度高效地测出电力系统谐波和间谐波的频率和幅值，且对背景噪声的抗干扰能力强。

为了达到上述技术效果，本发明的技术方案如下：

一种电力系统谐波和间谐波测量方法，包括以下步骤：

1)建立含M个谐波和间谐波的电力系统波形：

式中，A_i、f_i、分别为第i个谐波或间谐波的幅值、频率和初始相位，N(t)为随机噪声；

2)对电力系统波形进行采样，得到采样序列：

y(s)＝y(s·T)

式中，T为采样周期，s＝0，1，2，…为非负整数；

3)利用谱估计对采样序列进行功率谱分析：

$P_y\left(\mathrm{\ω}\right)={\mathrm{\σ}}_p^2/{\[1+\underset{k=1}{\overset{p}{\Σ}}{a}_{p,k}{e}^{-j\mathrm{\ω}k}\]}^2$

式中，为噪声序列的方差，ω为角频率，a_p,k(k＝1，2，…，p)为自回归谱估计模型的参数；利用上式初步得到基波、谐波及间谐波的数量及频率；

4)利用测得的基波、谐波及间谐波的频率，构造如下混沌检测振子，对其幅值进行测量：

$\stackrel{\·\·}{x}\left(t\right)=-{\mathrm{\ω}}_{i}c\stackrel{\·}{x}+{\mathrm{\ω}}_i^2(x-x^3+\mathrm{\γ}cos({\mathrm{\ω}}_{i}t\left)\right)$

式中，ω_i为测得的基波、谐波及间谐波角频率，c为阻尼比，γ为策动力的幅值；为提高幅值测量精度，利用最大李亚普诺夫指数λ及混沌状态和大尺度周期状态相变的临界条件作为混沌振子相变检测的依据，将不能引起所构造的混沌检测振子由大尺度周期状态向混沌状态相变的信号，作为“虚假”信号剔除，即得到精确的电力系统谐波和间谐波的频率和幅值。

进一步地，所述步骤4)中，混沌检测振子检测谐波的具体方法为：首先调节策动力幅值γ，使振子处于从混沌状态向大尺度周期状态过渡的临界状态，得到此时的策动力γ₁；然后，将包含各次谐波和间谐波在内的电力系统信号，作为策动力的一部分而并入系统，混沌振子立即进入大尺度周期稳定状态；同时，通过减小系统自身的策动力γ，使混沌振子再次处于大尺度周期到混沌的临界状态，根据变化后的策动力幅值γ₂，即得到角频率为ω_i的间谐波幅值为γ₂-γ₁。

进一步地，所述步骤4)中使用最大李亚普诺夫指数λ判定混沌系统运动状态的依据是：λ＞0时，系统处于混沌状态；λ＜0时，系统处于大尺度周期状态；λ＝0或λ≈0时，系统处于混沌与大尺度周期的临界状态，此时的策动力幅值即为系统发生相变的临界阈值。

进一步地，所述步骤4)中最大李亚普诺夫指数λ计算包括两步：相空间重构与李亚普诺夫指数计算：

①相空间重构

将时间维上的观测数据扩展到抽象的三维甚至更高维的相空间中去，这就是时间序列的相空间重构；设数据{xk|k＝1，2，3，…，N}为等时间间隔采样得到的时间序列，将其嵌入到m维欧式空间R_m中，得到一个点或向量的集Y，其元素记作：

Y_i＝(x_i,x_i+t,···,x_i+(m-1)t)

式中，i＝1，2，…，M，M＝N-(m-1)t；m为嵌入维数；t为时延序列，若时间序列的采样间隔为τ_s，相空间重构中的延迟时间τ_d表示为τ_d＝tτ_s；对于Y_i的重构相空间，其关联积分函数定义为：

$C(m,N,r,t)=\frac{2}{M(M-1)}\underset{1\≤i\≤j\≤M}{\Σ}\mathrm{\Θ}(r-d_{ij})$

式中，r>0，d_ij＝||Y_i-Y_j||，r为半径，Θ(a)为Heaviside阶跃函数：

$\mathrm{\&Theta;}\left(a\right)=\left\{\begin{array}{cc}0& a<0\\ 1& a\&GreaterEqual;0\end{array}\right.$

由关联积分函数得到描述非线性时间序列相关性的统计量：

S(m，N，r，t)＝C(m，N，r，t)-C^m(1，N，r，t)

统计量S(m，N，r，t)视为一个非线性依赖关系的无维数的度量，利用该统计量寻找到τ_d和m；

对于混沌时间序列{x_i}，i＝1，2，…，N，将其分成t个不相交的子时间序列，即：

{x₁，x_t+1，…，x_2t+1，…}

{x₂，x_t+2，…，x_2t+2，…}

……

{x_t，x_2t，x_3t，…}

其中每个子序列的S(m，N，r，t)采取求平均的策略，为：

$S(m,N,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\&Sigma;}}\&lsqb;C_s\left(m,\frac{N}{t},r,t\right)-C_s^m\left(1,\frac{N}{t},t\right)\&rsqb;$

令N→∞，有：

$S(m,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\&Sigma;}}\&lsqb;C_s\left(m,r,t\right)-C_s^m\left(1,r,t\right)\&rsqb;$

式中，m＝2，3，···如果时间序列是独立同分布的，那么对固定的m，t，当N→∞时，对于所有的r，均有S(m，r，t)＝0；但实际序列是有限的，并且序列元素间可能相关，因此，一般S(m，r，t)≠0；这样，局部最大时间间隔t取S(m，r，t)的零点或对所有的半径r相互差别最小的时间点，选择对应值最大和最小两个半径r，定义差量为：

ΔS(m，t)＝max{S(m，r_j，t)}-min{S(m，r_j，t)}

式中度量了关于r的最大偏差；由于S(m，r，t)的零点对所有m、r几乎相等，而ΔS(m，t)的最小值对所有的m也几乎相等，所以，延迟时间τ_d对应着这些局部最大时间t中的一个；当σ/2≤r≤2σ时，σ为BDS统计量的方差，对所有S(m，r，t)和ΔS(m，t)分别求平均得：

$\stackrel{\&OverBar;}{S}\left(t\right)=\frac{1}{m\&CenterDot;j}\underset{m}{\&Sigma;}\underset{j}{\&Sigma;}S(m,r_j,t)$

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{S}\left(t\right)=\frac{1}{m}\underset{m}{\&Sigma;}\mathrm{\&Delta;}S(m,t)$

其中，j为r_j的数目；取下述统计量：

$S_{cor}\left(t\right)=\mathrm{\&Delta;}\stackrel{\&OverBar;}{S}\left(t\right)+\left|\stackrel{\&OverBar;}{S}\left(t\right)\right|$

的最小值作为嵌入窗τ_ω的最优值，于是，最佳嵌入维数m为：

$m=\frac{{\mathrm{\&tau;}}_{\mathrm{\&omega;}}}{\mathrm{\&tau;}}+1$

②最大李亚普诺夫指数的计算

在重构相空间后，首先依据下式寻找给定轨道上每个点Y_j的最近临近点

d_j(0)＝min||Y_j-Y_i||

这里，必须是与Y_j不同轨道上的点，即满足：

式中，P为时间序列的平均周期；然后，对相空间中每个点Y_j，计算出该邻点对的i个离散时间步后的距离d_j(i)：

d_j(i)＝|Y_j-Y_i|

i＝1,2，···，由于最大李亚普诺夫指数λ体现的是相邻轨迹的平均分岔按指数规律增大，若用d(t)表示在不同轨迹上相邻两点的距离，则有：

$d\left(t\right)={\mathrm{Ce}}^{\mathrm{\&beta;}\left({\mathrm{i\&tau;}}_s\right)}$

其中，C为初始分岔，对重构轨迹上的某一点Y_j而言，C_j＝d_j(0)；所以有：

$d_j\left(t\right)=d_j\left(0\right)e^{\mathrm{\&beta;}\left({\mathrm{i\&tau;}}_s\right)}$

对上式两边取对数，有：

lnd_j(i)＝lnd_j(0)+β(i·τ_s)

该式表示一簇大致平行线，而每条线的斜率都大致与最大李亚普诺夫指数成正比，利用最小二乘法拟合这些线的“平均线”，便可求得最大李亚普诺夫指数；“平均线”为：

$y\left(i\right)=\frac{1}{i\&CenterDot;{\mathrm{\&tau;}}_s}<\ln d_j\left(i\right)>$

式中，<·>表示对所有j进行平均。

本发明提出的电力系统谐波和间谐波的高精度测量方法，首先应用谱估计初步测量电力系统波形中所含谐波和间谐波的数量及频率；然后基于测得的各谐波和间谐波的频率构造出混沌检测振子，对其幅值进行测量，为提高检测精度，应用最大李雅普诺夫指数作为混沌检测振子处于混沌状态和大尺度周期状态相变的临界条件的判断依据。

利用谱估计的谐波和间谐波频率测量

以采样频率f_s对式所述的电力系统波形进行采样，得：

式中，n为0,1,2，···的整数。

根据谱估计理论，上式可以转化为：

$y\left(n\right)=-\underset{k=1}{\overset{p}{\&Sigma;}}{a}_{p,k}y(n-k)+\mathrm{\&eta;}\left(n\right)$

式中，p为自回归谱估计模型的阶数，a_p,k(k＝1，2，…，p)为自回归谱估计模型的参数。可见，含有谐波和间谐波的电力系统波形信号可看作为自回归模型。因此，y(n)的功率谱可表示为：

$P_y\left(\mathrm{\&omega;}\right)={\mathrm{\&sigma;}}_p^2/{\&lsqb;1+\underset{k=1}{\overset{p}{\&Sigma;}}{a}_{p,k}{e}^{-j\mathrm{\&omega;}k}\&rsqb;}^2$

式中ω为角频率；为噪声序列的方差，数值上等于阶次为p时的最小预测误差功率ρ_p。式中表明，功率谱描述了信号功率随角频率的变化，功率谱的谱峰所对应的频率即为谐波和间谐波的频率值。为求得各谐波和间谐波的频率，只需求得自回归谱估计模型的参数和a_p,k。定义前、后向预测误差分别为：

$f_p\left(n\right)=\underset{k=0}{\overset{p}{\&Sigma;}}{a}_{p,k}y(n-k)$

$e_p\left(n\right)=\underset{k=0}{\overset{p}{\&Sigma;}}{a}_{p,k}y(n-p+k)$

计算各阶前、后向预测误差可利用格型滤波器结构递推：

$\left[\begin{array}{c}{f}_p\left(n\right)\\ e_p\left(n\right)\end{array}\right]=\left[\begin{array}{cc}1& {\mathrm{\&mu;}}_p\\ {\mathrm{\&mu;}}_p& 1\end{array}\right]\left[\begin{array}{c}{f}_{p-1}\left(n\right)\\ e_{p-1}(n-1)\end{array}\right]$

式中，μ_p为格型滤波器的反射系数。令前向和后向预测误差的平均功率为：

${\mathrm{\&rho;}}_p=\frac{1}{2}\underset{n=p}{\overset{N-1}{\&Sigma;}}\&lsqb;|f_p\left(n\right){|}^2+|e_p\left(n\right){|}^2\&rsqb;$

式中，N为电力系统波形的采样数据量，为使ρ_p最小，令得到反射系数：

${\mathrm{\&mu;}}_p=\frac{-2\underset{n=p+1}{\overset{N}{\&Sigma;}}\&lsqb;f_{p-1}\left(n\right)e_{p-1}\left(n-1\right)\&rsqb;}{\underset{n=p+1}{\overset{N}{\&Sigma;}}\&lsqb;{\left(f_{p-1}\left(n\right)\right)}^2+{\left(e_{p-1}\left(n-1\right)\right)}^2\&rsqb;}$

利用Levinson递推公式即下式即可求出自回归谱估计模型的参数：

$\left\{\begin{array}{c}{a}_{p,k}=a_{p-1,k}+{\mathrm{\&mu;}}_p{a}_{p-1,k-1},1\&le;k\&le;p-1\\ a_{p,p}={\mathrm{\&mu;}}_p\\ {\mathrm{\&rho;}}_p=(1-{\mathrm{\&mu;}}_p^2){\mathrm{\&rho;}}_{p-1}\end{array}\right.$

具体步骤为：

1、初始条件：e₀(n)＝f₀(n)＝y(n)；由式上式求出μ₁；

2、由 $\mathrm{\&rho;}\left(0\right)=\frac{1}{N}\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|y\left(n\right)\right|}^2$ 得k＝1时的参数：a_1，1＝μ₁， ${\mathrm{\&rho;}}_1=(1-{\mathrm{\&mu;}}_1^2)\mathrm{\&rho;}\left(0\right);$

3、由μ₁求出f₁(n)和e₁(n)，再求出μ₂；

4、求出k＝2时的a_2，1、a_2，2和ρ₂。

5、重复以上过程，直到k＝p，即可求出所有阶次的自回归谱估计模型参数。

利用混沌理论的谐波和间谐波幅值测量

利用测得的谐波和间谐波频率ω_i，构造如下混沌检测振子：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=-{\mathrm{\&omega;}}_{i}c\stackrel{\&CenterDot;}{x}+{\mathrm{\&omega;}}_i^2(x-x^3+\mathrm{\&gamma;}cos({\mathrm{\&omega;}}_{i}t\left)\right)$

式中，ω_i为测得的基波、谐波及间谐波角频率，c为阻尼比，γ为策动力的幅值；为提高幅值测量精度，利用最大李亚普诺夫指数λ及混沌状态和大尺度周期状态相变的临界条件作为混沌振子相变检测的依据，将不能引起所构造的混沌检测振子由大尺度周期状态向混沌状态相变的信号，作为“虚假”信号剔除，即得到精确的电力系统谐波和间谐波的频率和幅值。

进一步地，所述步骤4)中，混沌检测振子检测谐波的具体方法为：首先调节策动力幅值γ，使振子处于从混沌状态向大尺度周期状态过渡的临界状态，得到此时的策动力γ₁；然后，将包含各次谐波和间谐波在内的电力系统信号，作为策动力的一部分而并入系统，混沌振子立即进入大尺度周期稳定状态；同时，通过减小系统自身的策动力γ，使混沌振子再次处于大尺度周期到混沌的临界状态，根据变化后的策动力幅值γ₂，即得到角频率为ω_i的间谐波幅值为γ₂-γ₁。

进一步地，所述步骤4)中使用最大李亚普诺夫指数λ判定混沌系统运动状态的依据是：λ＞0时，系统处于混沌状态；λ＜0时，系统处于大尺度周期状态；λ＝0或λ≈0时，系统处于混沌与大尺度周期的临界状态，此时的策动力幅值即为系统发生相变的临界阈值。

进一步地，所述步骤4)中最大李亚普诺夫指数λ计算包括两步：相空间重构与李亚普诺夫指数计算：

①相空间重构

将时间维上的观测数据扩展到抽象的三维甚至更高维的相空间中去，这就是时间序列的相空间重构；设数据{x_k|k＝1，2，3，…，N}为等时间间隔采样得到的时间序列，将其嵌入到m维欧式空间R_m中，得到一个点或向量的集Y，其元素记作：

Y_i＝(x_i,x_i+t,···,x_i+(m-1)t)

式中，i＝1，2，…，M，M＝N-(m-1)t；m为嵌入维数；t为时延序列，若时间序列的采样间隔为τ_s，相空间重构中的延迟时间τ_d表示为τ_d＝tτ_s；对于Y_i的重构相空间，其关联积分函数定义为：

$C(m,N,r,t)=\frac{2}{M(M-1)}\underset{1\&le;i\&le;j\&le;M}{\&Sigma;}\mathrm{\&Theta;}(r-d_{ij})$

式中，r>0，d_ij＝||Y_i-Y_j||，r为半径，Θ(a)为Heaviside阶跃函数：

$\mathrm{\&Theta;}\left(a\right)=\left\{\begin{array}{cc}0& a<0\\ 1& a\&GreaterEqual;0\end{array}\right.$

由关联积分函数得到描述非线性时间序列相关性的统计量：

S(m，N，r，t)＝C(m，N，r，t)-C^m(1，N，r，t)

统计量S(m，N，r，t)视为一个非线性依赖关系的无维数的度量，利用该统计量寻找到τ_d和m；

对于混沌时间序列{x_i}，i＝1，2，…，N，将其分成t个不相交的子时间序列，即：

{x₁，x_t+1，…，x_2t+1，…}

{x₂，x_t+2，…，x_2t+2，…}

……

{x_t，x_2t，x_3t，…}

其中每个子序列的S(m，N，r，t)采取求平均的策略，为：

$S(m,N,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\&Sigma;}}\&lsqb;C_s\left(m,\frac{N}{t},r,t\right)-C_s^m\left(1,\frac{N}{t},t\right)\&rsqb;$

令N→∞，有：

$S(m,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\&Sigma;}}\&lsqb;C_s\left(m,r,t\right)-C_s^m\left(1,r,t\right)\&rsqb;$

式中，m＝2，3，···如果时间序列是独立同分布的，那么对固定的m，t，当N→∞时，对于所有的r，均有S(m，r，t)＝0；但实际序列是有限的，并且序列元素间可能相关，因此，一般S(m，r，t)≠0；这样，局部最大时间间隔t取S(m，r，t)的零点或对所有的半径r相互差别最小的时间点，选择对应值最大和最小两个半径r，定义差量为：

ΔS(m，t)＝max{S(m，r_j，t)}-min{S(m，r_j，t)}

式中度量了关于r的最大偏差；由于S(m，r，t)的零点对所有m、r几乎相等，而ΔS(m，t)的最小值对所有的m也几乎相等，所以，延迟时间τ_d对应着这些局部最大时间t中的一个；当σ/2≤r≤2σ时，σ为BDS统计量的方差，对所有S(m，r，t)和ΔS(m，t)分别求平均得：

$\stackrel{\&OverBar;}{S}\left(t\right)=\frac{1}{m\&CenterDot;j}\underset{m}{\&Sigma;}\underset{j}{\&Sigma;}S(m,r_j,t)$

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{S}\left(t\right)=\frac{1}{m}\underset{m}{\&Sigma;}\mathrm{\&Delta;}S(m,t)$

其中，j为r_j的数目；取下述统计量：

$S_{cor}\left(t\right)=\mathrm{\&Delta;}\stackrel{\&OverBar;}{S}\left(t\right)+\left|\stackrel{\&OverBar;}{S}\left(t\right)\right|$

的最小值作为嵌入窗τ_ω的最优值，于是，最佳嵌入维数m为：

$m=\frac{{\mathrm{\&tau;}}_{\mathrm{\&omega;}}}{\mathrm{\&tau;}}+1$

②最大李亚普诺夫指数的计算

在重构相空间后，首先依据下式寻找给定轨道上每个点Y_j的最近临近点

d_j(0)＝min||Y_j-Y_i||

这里，必须是与Y_j不同轨道上的点，即满足：

式中，P为时间序列的平均周期；然后，对相空间中每个点Y_j，计算出该邻点对的i个离散时间步后的距离d_j(i)：

d_j(i)＝|Y_j-Y_i|

式中，i＝1,2，···，由于最大李亚普诺夫指数λ体现的是相邻轨迹的平均分岔按指数规律增大，若用d(t)表示在不同轨迹上相邻两点的距离，则有：

$d\left(t\right)={\mathrm{Ce}}^{\mathrm{\&beta;}\left({\mathrm{i\&tau;}}_s\right)}$

其中，C为初始分岔，对重构轨迹上的某一点Y_j而言，C_j＝d_j(0)；所以有：

$d_j\left(t\right)=d_j\left(0\right)e^{\mathrm{\&beta;}\left({\mathrm{i\&tau;}}_s\right)}$

对上式两边取对数，有：

lnd_j(i)＝lnd_j(0)+β(i·τ_s)

该式表示一簇大致平行线，而每条线的斜率都大致与最大李亚普诺夫指数成正比，利用最小二乘法拟合这些线的“平均线”，便可求得最大李亚普诺夫指数；“平均线”为：

$y\left(i\right)=\frac{1}{i\&CenterDot;{\mathrm{\&tau;}}_s}<\ln d_j\left(i\right)>$

式中，<·>表示对所有j进行平均。

与现有技术相比，本发明技术方案的有益效果是：

本发明采用谱估计和混沌理论对电力系统波形中的谐波和间谐波信号的频率和幅值进行检测，具有测量精度高、检测速度快的优点，不仅能克服快速傅立叶变换方法的频谱混叠、栅栏效应及频谱泄露等问题，且能抑制背景噪声干扰。

附图说明

图1为含基波、谐波和间谐波在内的电力系统波形的功率谱分析。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；

为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；

对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

下面结合附图和实施例对本发明的技术方案做进一步的说明。

实施例1

设包含谐波和间谐波在内的电力系统波形如下式所述：

式中，A_i、f_i、分别为第i个谐波或间谐波的幅值、频率和初始相位，各参数的具体值如表1所示；N(t)是均值为0，标准差为0.5的随机噪声。

步骤1)以固定采样频率1000Hz(采样周期T＝1/1000秒)对上式表示的电力系统波形进行采样，得到采样序列。

y(s)＝y(s·T)s＝0，1，2，…

步骤2)利用谱估计方法对采样序列进行功率谱分析(如图1所示)，初步得到基波及各谐波、间谐波的数量及频率。由图1可以看出，除表1中的基波、谐波和间谐波外，在187.5Hz处还存在一个小的频率峰，对于该“虚假”信号，利用步骤3)中的混沌检测可以剔除。由表1所测得的基波、谐波及间谐波的频率测量值可见，该法的频率测量精度较高。

表1各(间)谐波参数的实际值和测量值

	间谐波1	基波	间谐波2	三次谐波	间谐波3
						实际频率	31.6	50	120.4	149.8	159.2
测量频率	31.1	49.2	120.2	149.2	159.2
						实际幅值Ai	1.2	104	0.78	2.4	2.2
测量幅值	1.1	102	0.811	2.27	2.084

谐波和间谐波的准确、可靠分析对电力系统的安全和经济运行具有重要意义，本发明建立了基于谱估计和混沌理论相结合的电力系统谐波和间谐波测量方法，该方法能够高精度地测量电力系统谐波和间谐波的频率和幅值，且能够克服噪声干扰，不受栅栏效应、频谱泄露等因素的影响。

步骤3)利用测得的基波、谐波及间谐波频率依次构造混沌检测振子，对初步测得的基波、谐波和间谐波的幅值进行测量，为提高间谐波幅值的检测精度，利用最大李亚普诺夫指数作为混沌系统相变检测的依据，测量得到的基波及各间谐波的幅值如表1所示。由于图1中187.5Hz处的频率峰为随机噪声引起，不能引起所构造的混沌检测振子由大尺度周期状态向混沌状态转变的相变，因此，可将该“虚假”信号剔除。

相同或相似的标号对应相同或相似的部件；

附图中描述位置关系的用于仅用于示例性说明，不能理解为对本专利的限制；

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (4)

1.一种电力系统谐波和间谐波测量方法，其特征在于，包括以下步骤：

1)建立含M个谐波和间谐波的电力系统波形：

式中，A_i、f_i、分别为第i个谐波或间谐波的幅值、频率和初始相位，N(t)为随机噪声；

2)对电力系统波形进行采样，得到采样序列：

y(s)＝y(s·T)

式中，T为采样周期，s＝0，1，2，…为非负整数；

3)利用谱估计对采样序列进行功率谱分析：

$P_y\left(\mathrm{\&omega;}\right)={\mathrm{\&sigma;}}_p^2/\&lsqb;1+{\underset{k=1}{\overset{p}{\&Sigma;}}{a}_{p,k}{e}^{-j\mathrm{\&omega;}k}\&rsqb;}^2$

式中，为噪声序列的方差，ω为角频率，a_p,k(k＝1，2，…，p)为自回归谱估计模型的参数；利用上式初步得到基波、谐波及间谐波的数量及频率；

4)利用测得的基波、谐波及间谐波的频率，构造如下混沌检测振子，对其幅值进行测量：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=-{\mathrm{\&omega;}}_{i}c\stackrel{\&CenterDot;}{x}+{\mathrm{\&omega;}}_i^2(x-x^3+\mathrm{\&gamma;}cos({\mathrm{\&omega;}}_{i}t\left)\right)$

式中，ω_i为测得的基波、谐波及间谐波角频率，c为阻尼比，γ为策动力的幅值；为提高幅值测量精度，利用最大李亚普诺夫指数λ及混沌状态和大尺度周期状态相变的临界条件作为混沌振子相变检测的依据，将不能引起所构造的混沌检测振子由大尺度周期状态向混沌状态相变的信号，作为“虚假”信号剔除，即得到精确的电力系统谐波和间谐波的频率和幅值。

2.根据权利要求1所述的电力系统谐波和间谐波测量方法，其特征在于，所述步骤4)中，混沌检测振子检测谐波的具体方法为：首先调节策动力幅值γ，使振子处于从混沌状态向大尺度周期状态过渡的临界状态，得到此时的策动力γ₁；然后，将包含各次谐波和间谐波在内的电力系统信号，作为策动力的一部分而并入系统，混沌振子立即进入大尺度周期稳定状态；同时，通过减小系统自身的策动力γ，使混沌振子再次处于大尺度周期到混沌的临界状态，根据变化后的策动力幅值γ₂，即得到角频率为ω_i的间谐波幅值为γ₂-γ₁。

3.根据权利要求1所述的电力系统谐波和间谐波测量方法，其特征在于，所述步骤4)中使用最大李亚普诺夫指数λ判定混沌系统运动状态的依据是：λ＞0时，系统处于混沌状态；λ＜0时，系统处于大尺度周期状态；λ＝0或λ≈0时，系统处于混沌与大尺度周期的临界状态，此时的策动力幅值即为系统发生相变的临界阈值。

4.根据权利要求1所述的电力系统谐波和间谐波测量方法，其特征在于，所述步骤4)中最大李亚普诺夫指数λ计算包括两步：相空间重构与李亚普诺夫指数计算：

①相空间重构

将时间维上的观测数据扩展到抽象的三维甚至更高维的相空间中去，这就是时间序列的相空间重构；设数据{x_k|k＝1，2，3，…，N}为等时间间隔采样得到的时间序列，将其嵌入到m维欧式空间R_m中，得到一个点或向量的集Y，其元素记作：

Y_i＝(x_i,x_i+t,···,x_i+(m-1)t)

式中，i＝1，2，…，M，M＝N-(m-1)t；m为嵌入维数；t为时延序列，若时间序列的采样间隔为τ_s，相空间重构中的延迟时间τ_d表示为τ_d＝tτ_s；对于Y_i的重构相空间，其关联积分函数定义为：

$C(m,N,r,t)=\frac{2}{M(M-1)}\underset{1\&le;i\&le;j\&le;M}{\&Sigma;}\mathrm{\&Theta;}(r-d_{ij})$

式中，r>0，d_ij＝||Y_i-Y_j||，r为半径，Θ(a)为Heaviside阶跃函数：

$\mathrm{\&Theta;}\left(a\right)=\left\{\begin{array}{cc}0& a<0\\ 1& a\&GreaterEqual;0\end{array}\right.$

由关联积分函数得到描述非线性时间序列相关性的统计量：

S(m，N，r，t)＝C(m，N，r，t)-C^m(1，N，r，t)

统计量S(m，N，r，t)视为一个非线性依赖关系的无维数的度量，利用该统计量寻找到τ_d和m；

对于混沌时间序列{x_i}，i＝1，2，…，N，将其分成t个不相交的子时间序列，即：

{x₁，x_t+1，…，x_2t+1，…}

{x₂，x_t+2，…，x_2t+2，…}

……

{x_t，x_2t，x_3t，…}

其中每个子序列的S(m，N，r，t)采取求平均的策略，为：

$S(m,N,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\&Sigma;}}\&lsqb;C_s\left(m,\frac{N}{t},r,t\right)-C_s^m\left(1,\frac{N}{t},t\right)\&rsqb;$

令N→∞，有：

$S(m,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\&Sigma;}}\&lsqb;C_s\left(m,r,t\right)-C_s^m\left(1,r,t\right)\&rsqb;$

式中，m＝2，3，···如果时间序列是独立同分布的，那么对固定的m，t，当N→∞时，对于所有的r，均有S(m，r，t)＝0；但实际序列是有限的，并且序列元素间可能相关，因此，一般S(m，r，t)≠0；这样，局部最大时间间隔t取S(m，r，t)的零点或对所有的半径r相互差别最小的时间点，选择对应值最大和最小两个半径r，定义差量为：

ΔS(m，t)＝max{S(m，r_j，t)}-min{S(m，r_j，t)}

式中度量了关于r的最大偏差；由于S(m，r，t)的零点对所有m、r几乎相等，而ΔS(m，t)的最小值对所有的m也几乎相等，所以，延迟时间τ_d对应着这些局部最大时间t中的一个；当σ/2≤r≤2σ时，σ为BDS统计量的方差，对所有S(m，r，t)和ΔS(m，t)分别求平均得：

$\stackrel{\&OverBar;}{S}\left(t\right)=\frac{1}{m\&CenterDot;j}\underset{m}{\&Sigma;}\underset{j}{\&Sigma;}S(m,r_j,t)$

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{S}\left(t\right)=\frac{1}{m}\underset{m}{\&Sigma;}\mathrm{\&Delta;}S(m,t)$

其中，j为r_j的数目；取下述统计量：

$S_{cor}\left(t\right)=\mathrm{\&Delta;}\stackrel{\&OverBar;}{S}\left(t\right)+\left|\stackrel{\&OverBar;}{S}\left(t\right)\right|$

的最小值作为嵌入窗τ_ω的最优值，于是，最佳嵌入维数m为：

$m=\frac{{\mathrm{\&tau;}}_{\mathrm{\&omega;}}}{\mathrm{\&tau;}}+1$

②最大李亚普诺夫指数的计算

在重构相空间后，首先依据下式寻找给定轨道上每个点Y_j的最近临近点

d_j(0)＝min||Y_j-Y_i||

这里，必须是与Y_j不同轨道上的点，即满足：

式中，P为时间序列的平均周期；然后，对相空间中每个点Y_j，计算出该邻点对的i个离散时间步后的距离d_j(i)：

d_j(i)＝|Y_j-Y_i|

式中，i＝1,2，···，min(M-j，)由于最大李亚普诺夫指数λ体现的是相邻轨迹的平均分岔按指数规律增大，若用d(t)表示在不同轨迹上相邻两点的距离，则有：

$d\left(t\right)={\mathrm{Ce}}^{\mathrm{\&beta;}\left({\mathrm{i\&tau;}}_s\right)}$

其中，C为初始分岔，对重构轨迹上的某一点Y_j而言，C_j＝d_j(0)；所以有：

$d_j\left(t\right)=d_j\left(0\right)e^{\mathrm{\&beta;}\left({\mathrm{i\&tau;}}_s\right)}$

对上式两边取对数，有：

lnd_j(i)＝lnd_j(0)+β(i·τ_s)

该式表示一簇大致平行线，而每条线的斜率都大致与最大李亚普诺夫指数成正比，利用最小二乘法拟合这些线的“平均线”，便可求得最大李亚普诺夫指数；“平均线”为：

$y\left(i\right)=\frac{1}{i\&CenterDot;{\mathrm{\&tau;}}_s}<\ln d_j\left(i\right)>$

式中，<·>表示对所有j进行平均。