CN102323476B - 采用谱估计和混沌理论的电力系统谐波和间谐波测量方法 - Google Patents

Abstract

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

Description

采用谱估计和混沌理论的电力系统谐波和间谐波测量方法

技术领域

本发明专利涉及一种采用谱估计和混沌理论的电力系统谐波和间谐波测量方法。

技术背景

电力系统理想的电压、电流波形是正弦波，但由于电力系统中存在各种非线性元件，使电压和电流波形发生畸变而产生谐波和间谐波(即具有非整数倍基波频率的信号分量)。谐波(间谐波)的存在，不仅会造成电网的功率损耗增加、设备寿命缩短、保护功能失常，还会引起变电站局部并联或串联谐振等。近年来，随着电力电子技术的飞速发展，越来越多的谐波源负载和装置在电力系统、工业、交通及家庭中的应用日益广泛，谐波已被认为是电网的一大公害。因此，采取有效措施抑制系统中的谐波和间谐波，对提高电力系统的安全与经济运行具有重要意义。

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

发明内容

本发明针对以上方法的不足，提供一种采用谱估计和混沌理论的电力系统谐波和间谐波测量方法，该方法能够以较高的频率分辨率和幅值测量精度高效地测出电力系统谐波和间谐波的频率和幅值，且对背景噪声的抗干扰能力强。

为实现上述目的，本发明的技术方案是：

一种采用谱估计和混沌理论的电力系统谐波和间谐波测量方法，包括以下步骤：

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

式中，A_i、f_i、

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

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

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

式中，T为采样周期，s为非负整数；

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

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

式中，

为噪声序列的方差，ω为角频率，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 \left({\mathrm{\ω}}_{i}t\right))---\left(4\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    (15)

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

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

式中，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.---\left(17\right)$

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

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

统计量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，…}    (19)

……

{x_t，x_2t，x_3t，…}

根据式(18)，对式(19)中每个子序列的S(m，N，r，t)采取求平均的策略，为：

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

令N→∞，有：

$S(m,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\mathrm{\&Sigma;}}}[C_s(m,r,t)-C_s^m(1,r,t)],(m=\mathrm{2,3},...)---\left(21\right)$

如果时间序列是独立同分布的，那么对固定的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)}    (22)

式(22)度量了关于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}{\mathrm{\&Sigma;}}\underset{j}{\mathrm{\&Sigma;}}S(m,r_j,t)---\left(23\right)$

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

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

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

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

$m=\frac{{\mathrm{\&tau;}}_w}{\mathrm{\&tau;}}+1---\left(26\right)$

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

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

$d_j\left(0\right)=\min \left|\right|Y_j-Y_{\hat{j}}\left|\right|---\left(27\right)$

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

$|j-\hat{j}|>P---\left(28\right)$

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

$d_j\left(i\right)=|Y_{j+i}-Y_{\hat{j}+i}|,i=\mathrm{1,2},...,\min (M-j,M-\hat{j})---\left(29\right)$

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

d(t)＝Ce^βt                    (30)

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

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

对上式两边取对数，有：

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

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

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

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

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

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

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

式中：n为0，1，2，……的整数。

根据谱估计理论，式(5)可以转化为：

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

式中，

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

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

式中，ω为角频率；

为噪声序列的方差，数值上等于阶次为p时的最小预测误差功率ρ_p。式(7)表明，功率谱描述了信号功率随角频率的变化，功率谱的谱峰所对应的频率即为谐波和间谐波的频率值。为求得各谐波和间谐波的频率，只需求得自回归谱估计模型的参数

和α_p，k。定义前、后向预测误差分别为：

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

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

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

$\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]---\left(10\right)$

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

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

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

得到反射系数：

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

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

$\left\{\begin{array}{c}{a}_{p,k}=a_{p-1,k}+{\mathrm{\&mu;}}_p{a}_{p-1,p-k},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.---\left(13\right)$

以上算法的具体步骤为：

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

(2)由

得k＝1时的参数：a_1，1＝μ₁，

(3)由μ₁和式(10)求出f₁(n)和e₁(n)，再由式(14)求出μ₂。

(4)依据式(10)、(12)、(13)的递推关系，求出k＝2时的a_2，1、a_2，2和ρ₂。

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

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

利用测得的谐波和间谐波频率ω_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 \left({\mathrm{\&omega;}}_{i}t\right))---\left(14\right)$

式中，c为阻尼比，γ为策动力的幅值，ω_i为谐波和间谐波的角频率。利用式(14)所示的混沌检测振子检测谐波的具体思路为：首先调节策动力幅值γ，使振子处于从混沌状态向大尺度周期状态过渡的临界状态，得到此时的策动力γ₁；然后，将包含各次谐波和间谐波在内的电力系统信号，作为策动力的一部分而并入系统，由于混沌振子对与自身策动力同频的信号极其敏感，混沌振子立即进入大尺度周期稳定状态；同时，由于振子对噪声具有免疫特性，即振子的相变只与加入系统的同频信号有关，而与噪声无关，因此，通过减小系统自身的策动力γ，使振子再次处于大尺度周期到混沌的临界状态，根据变化后的策动力幅值γ₂，即可求得角频率为ω_i的间谐波幅值为γ₂-γ₁。

上述幅值检测过程的关键是确定混沌振子混沌状态和大尺度周期状态转变的临界状态。为准确确定混沌振子的临界状态，提高谐波和间谐波幅值的测量精度，本发明使用最大李亚普诺夫指数λ作为混沌振子的相变在量上的判断依据，从而准确确定振子临界状态的策动力阈值。最大李亚普诺夫指数判定混沌系统运动状态的依据是：λ＞0时，系统处于混沌状态；λ＜0时，系统处于大尺度周期状态；λ＝0或λ≈0时，系统处于混沌与大尺度周期的临界状态，此时的策动力幅值即为系统发生相变的临界阈值。

最大李亚普诺夫指数计算主要包括两步：相空间重构与李亚普诺夫指数计算。

①相空间重构

将时间维上的观测数据扩展到抽象的三维甚至更高维的相空间中去，这就是时间序列的相空间重构。设数据{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    (15)

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

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

式中，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.---\left(17\right)$

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

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

统计量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，…}    (19)

……

{x_t，x_2t，x_3t，…}

根据式(18)，对式(19)中每个子序列的S(m，N，r，t)采取求平均的策略，为：

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

令N→∞，有：

$S(m,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\mathrm{\&Sigma;}}}[C_s(m,r,t)-C_s^m(1,r,t)],(m=\mathrm{2,3},...)---\left(21\right)$

如果时间序列是独立同分布的，那么对固定的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)}    (22)

式(22)度量了关于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}{\mathrm{\&Sigma;}}\underset{j}{\mathrm{\&Sigma;}}S(m,r_j,t)---\left(23\right)$

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

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

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

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

$m=\frac{{\mathrm{\&tau;}}_w}{\mathrm{\&tau;}}+1---\left(26\right)$

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

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

$d_j\left(0\right)=\min \left|\right|Y_j-Y_{\hat{j}}\left|\right|---\left(27\right)$

这里，

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

$|j-\hat{j}|>P---\left(28\right)$

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

$d_j\left(i\right)=|Y_{j+i}-Y_{\hat{j}+i}|,i=\mathrm{1,2},...,\min (M-j,M-\hat{j})---\left(29\right)$

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

d(t)＝Ce^βt    (30)

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

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

对上式两边取对数，有：

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

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

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

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

本发明的优点与效果

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

附图说明

图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所测得的基波、谐波及间谐波的频率测量值可见，该法的频率测量精度较高。

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

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

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

Claims (4)

1.一种采用谱估计和混沌理论的电力系统谐波和间谐波测量方法，其特征在于，包括以下步骤：

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

式中，A_i、f_i、

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

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

y(x)＝y(s·T)s＝0,1,2,…        （2）式中，T为采样时间间隔，s为非负整数；

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

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

为噪声序列的方差，ω为角频率，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 \left({\mathrm{\&omega;}}_{i}t\right))---\left(4\right)$ 式中，ω_i为测得的基波、谐波及间谐波角频率，c为阻尼比，γ为策动力的幅值；为提高幅值测量精度，利用最大李亚普诺夫指数λ及混沌状态和大尺度周期状态相变的临界条件作为混沌检测振子相变检测的依据，将不能引起所构造的混沌检测振子由大尺度周期状态向混沌状态相变的信号，作为“虚假”信号剔除，即得到精确的电力系统谐波和间谐波的频率和幅值。

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

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

4.如权利要求3所述的方法，其特征在于：所述步骤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   （15）式中，M＝N-(m-1)t；m为嵌入维数；t为时延序列，若时间序列的采样间隔为τ_s，相空间重构中的延迟时间τ_d表示为τ_d＝tτ_s；对于式(15)的重构相空间，其关联积分函数定义为：

$C(m,N,r,t)=\frac{2}{M(M-1)}\underset{1\&le;i\&le;j\&le;M}{\mathrm{\&Sigma;}}\mathrm{\&Theta;}(r-d_{\mathrm{ij}}),r>0---\left(16\right)$ 式中，d_ij＝||Y_i-Y_j||，r为半径，

为Heaviside阶跃函数：

$\mathrm{\&Theta;}\left(a\right)=\left\{\begin{array}{cc}0& a<0\\ 1& a\&GreaterEqual;0\end{array}\right.---\left(17\right)$ 由关联积分函数得到描述非线性时间序列相关性的统计量：

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

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

$\begin{array}{c}\{x_1,X_{t+1},\&CenterDot;\&CenterDot;\&CenterDot;,X_{2t+1},\&CenterDot;\&CenterDot;\&CenterDot;\}\\ \{x_2,x_{t+2},\&CenterDot;\&CenterDot;\&CenterDot;,x_{2t+2},\&CenterDot;\&CenterDot;\&CenterDot;\}\\ \&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\\ \{x_t,x_{2t},x_{3t},\&CenterDot;\&CenterDot;\&CenterDot;\}\end{array}---\left(19\right)$

根据式(18)，对式(19)中每个子序列的S(m,N,r,t)采取求平均的策略，为：

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

令N→∞，有：

$S(m,r,t)=\frac{1}{t}\underset{s=1}{\overset{t}{\mathrm{\&Sigma;}}}[C_s(m,r,t)-C_s^m(1,r,t)],(m=2,3,\&CenterDot;\&CenterDot;\&CenterDot;)---\left(21\right)$

如果时间序列是独立同分布的，那么对固定的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)}   （22）式(22)度量了关于r的最大偏差；由于S(m,r,t)的零点对所有m、r几乎相等，而ΔS(m,t)的最小值对所有的m也几乎相等，所以，延迟时间τ_d对应着这些局部最大时间t中的一个；当σ/2≤r≤2σ时，σ为BDS统计量的方差，BDS统计量为公式（18）表述的统计量，对所有S(m,r,t)和ΔS(m,t)分别求平均得：

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

$\mathrm{\&Delta;}\stackrel{\_}{S}\left(t\right)=\frac{1}{m}\underset{m}{\mathrm{\&Sigma;}}\mathrm{\&Delta;S}(m,t)---\left(24\right)$ 其中，j为r_j的数目；取下述统计量：

$S_{\mathrm{cor}}\left(t\right)=\mathrm{\&Delta;}\stackrel{\_}{S}\left(t\right)+\left|\stackrel{\_}{S}\left(t\right)\right|---\left(25\right)$ 的最小值作为嵌入窗τ_w的最优值，于是，最佳嵌入维数m'为：

$m\&prime;=\frac{{\mathrm{\&tau;}}_w}{{\mathrm{\&tau;}}_d}+1---\left(26\right)$

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

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

$d_j\left(0\right)=\min \left|\right|Y_j-Y_{\hat{j}}\left|\right|---\left(27\right)$

这里，

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

$|j-\hat{j}|>P---\left(28\right)$ 式中，P为时间序列的平均周期；然后，对相空间中每个点Y_j，计算出邻点对的i个离散时间步后的距离d_j(i)：

$d_j\left(i\right)=|Y_{j+i}-Y_{\hat{j}+1}|,i=\mathrm{1,2}\&CenterDot;\&CenterDot;\&CenterDot;,\min (M-j,M-\hat{j})---\left(29\right)$

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

d(t)＝Ce^λt   （30）

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

$d_j\left(i\right)=d_j\left(0\right)e^{\mathrm{\&lambda;}\left({\mathrm{i\&tau;}}_s\right)}---\left(31\right)$

对上式两边取对数，有：

lnd_j(i)＝lnd_j(0)+λ(i·τ_s)   （32）

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

$y\left(i\right)=\frac{1}{i\&CenterDot;{\mathrm{\&tau;}}_s}\&lang;\ln d_j\left(i\right)\&rang;---\left(33\right)$

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

Images