CN103245830A - 一种结合ar谱估计与非线性优化的间谐波检测方法 - Google Patents

Abstract

一种结合AR谱估计与非线性优化的间谐波检测方法，其特征应用AR模型谱估计计算谐波初始参数来建立信号模型：选取待分析的采样信号序列，对其进行AR模型谱估计，给出精度较低的谐波参数(频率、幅值和相位)信息；再取模型，提出求解信号模型参数的非线性优化算法，对相位值进行计算。本发明的有益效果是：解决了非线性优化算法的参数建模问题，克服了AR谱估计不能计算谐波幅值和相位的缺点，进一步提高了频率的计算精度；降低迭代算法对初值的敏感度，提高了迭代稳定性和计算效率；且有一定的抗噪能力。

Description

一种结合AR谱估计与非线性优化的间谐波检测方法

技术领域

本发明涉及一种电网信号间谐波分析领域，具体的说，涉及一种结合AR谱估计与非线性优化的间谐波检测方法。

背景技术

电力电子等非线性设备在电力系统中的广泛应用，谐波和间谐波日益增多，因而分析谐波和间谐波对电力系统的经济安全运行有着重要意义。目前的谐波分析主要通过对电网信号采样和数字化处理实现的，如快速傅里叶变换(FFT)和加窗插值类FFT方法。但这些方法在检测谐波时不可避免会产生频谱泄漏和栅栏效应，影响检测精度，用于间谐波(非整数次谐波)检测时效果更差，原因在于间谐波常出现在整数次谐波附近的频点且幅值较小，易被整数次谐波的旁瓣淹没，导致其频率分辨率大大降低，而低分辨率造成的谱线识别错误又会进一步地降低幅值与相位的检测精度。

现有技术的缺点是：难以获得电力线路中各次谐波的频率、幅值和相位的准确值，无法描述出含谐波干扰的电路波形。

发明内容

本发明针对现有技术的不足，提供了一种结合AR谱估计与非线性优化的间谐波检测方法，用于计算谐波参数的非线性优化算法，实现电网信号的频率、幅值和相位的准确计算。

为了实现上述发明目的，本发明提供一种结合AR谱估计与非线性优化的间谐波检测方法。

一种结合AR谱估计与非线性优化的间谐波检测方法，其特征在于,按以下步骤进行：

步骤一：AR谱初步估计

1.1）：根据谐波参数的初步信息建立该谐波的参数模型，取参数模型如下：

$y\left(t\right)\underset{m}{\mathrm{\Σ}}{A}_m\sin (2\mathrm{\π}{\mathrm{fr}}_{m}t+{\mathrm{\θ}}_m)---\left(1\right)$

式中y(t)为待分析信号，A_m为各次谐波和间谐波的幅值，fr_m为对应的频率，θ_m为对应的相位，m为对应的谐波个数，t为时间；

1.2）：对求解信号模型参数进行非线性优化算法：将所述模型（1）作如下处理，将幅值、相位(A_m,θ_m)等价转化为求解两个幅值(C_m,D_m)：

A_msin(2πfr_mt+θ_m)=A_msin(2πfr_mt)cos(θ_m)+A_mcos(2πfr_mt)sin(θ_m)

　　　　　　　　=C_msin(2πfr_mt)+D_mcos(2πfr_mt)    (2)

其中，C_m=A_mcos(θ_m)，D_m=A_msin(θ_m)；

1.3）：针对步骤一中的参数模型（1），将式(2)代入式(1)，建立最小二乘方程如下：

式中n为采样点数；目标函数

表征了信号计算结果和实际采样结果的误差，令y(t_n)为信号第n点的计算值、y₁(t_n)为信号第n点实际采样值，将参数f_m、C_m和D_m统一用向量x表示，对应于各采样时刻的误差向量用f(x)表示：

f(x)=[f₁,f₂，...，f_n]

　　　=[y(t₁)-y₁(t₁),(t₂)-y₁(t₂)，...，(t_n)-y₁(t_n)]

代入（3）式可得：

式中[f(x)]^T表示f(x)的转置，则参数的求解转化为求如下最值的问题：

步骤二：共轭梯度法一次迭代，也就是参数初值进行预先处理，使参数迭代值接近全局最优解。

2.1）对参数初值进行迭代，迭代方程为：

x^k+1=x^k+λ_kp_k    (5)

其中，k为迭代次数，λ_k为最优步长，p_k为迭代方向，

2.2）p_k满足下式（7）

置初始近似值x⁰∈Rⁿ（Rⁿ为实数域），一次迭代误差限ε₁>0

2.3）：计算，令0→k，p₀是迭代方向矩阵p_k的初值；

2.4）：求最优步长λ_k满足，计算迭代值x^k+1=x^k+λ_kp_k以及

2.5）：计算

，当迭代到

时，转步骤三，否则，转2.6；

2.6）：计算

k+1→k，转2.4；

步骤三：当

后，再用阻尼最小二乘法二次迭代：即采用阻尼最小二乘法快速寻找到全局最优解，即最终的间谐波参数计算值。

3.1）：为使沿着其梯度方向减小，必有

的法方程等于零：

可得：

式中：G(x^b)=Df(x^b)^TDf(x^b)，g(x^b)=Df(x^b)^Tf(x^b)，p(x^b)=-G(x^b)^-1g(x^b)；

3.2）：考虑阻尼时有：

其中μ_b为第b次迭代的阻尼因子，I是与G(x^b)同维数的单位对角阵；则最终迭代公式为：

当μ_b>0总可以保证

是收敛的，原则上μ_b=10^-4～10^-2；

3.3）：计算得出f(x^b),Df(x^b),

令0→j（j为3.7步骤所需的一个判断变量）；

3.4）：解方程组

求得

3.5）：计算x^b+1=x^b+p(x^b)及

3.6）：若

则取

（缩放因子v>1，可取2,5,10），转3.3），否则转3.7）；

3.7）：若j=0，则取μ_b=vμ_b,1→j，转3.3），否则转3.8）；

3.8）：若||p(x^b)||≤ε₂，则x^b+1为极小值点x^*的近似，算法停止，输出x^b+1；否则取x^b+1→x^b，b+1→b，转3.3）。

本发明所述步骤一中的参数模型式（1）中的m次数可根据功率谱中局部最大值个数确定，然后依次确定功率谱中每一个局部最大值处对应的频率和局部最大值的开方，把它们分别作为参数f_m、A_m的初值，然后用0～2π内的随机值作为相位参数θ_m的初值。

本发明所述步骤三中的3.2）阻尼因子μ_b，当μ_b>0总可以保证

因而是收敛的，当μ_b太大时会收敛速度下降，若μ_b太小则收敛域过小，原则上μ_b=10^-4～10^-2，通常取μ_b=10^-2。

本发明AR模型谱估计有较好的频率分辨率，能检测出间谐波频率，但其幅值参数计算结果与实际值相距甚远且不能给出相位参数，所以需结合非线性优化算法来搜寻准确的谐波参数；应用AR谱估计初步估算信号的谐波参数，给出频率和幅值估值，相位初始参数在0～2π内随机选取；然后基于前述参数建立信号模型，再应用非线性优化算法求取最终的模型参数，给出间谐波检参数检测结果。在非线性优化算法中，为克服阻尼最小二乘法对初值的敏感而选择共轭梯度法对参数初值进行预先处理，使参数迭代值接近全局最优解，再用阻尼最小二乘法快速寻找到全局最优解，即最终的间谐波参数计算值。

式（1）中各次谐波的幅值A_m和频率fr_m是可以准确获得，但相位θ_m为假设值，后续步骤就是为使相位θ_m与真实值相吻合。

步骤三：当

后，再用阻尼最小二乘法二次迭代：即采用阻尼最小二乘法快速寻找到全局最优解，即最终的间谐波参数计算值。

3.1）：由于Df(x)为f(x)的偏导矩阵：

$\mathrm{Df}{\left(x\right)}^T=\left[\begin{array}{cccc}\frac{{\∂f}_1}{\∂x_1}& \frac{{\∂f}_2}{{\∂x}_1}& ...& \frac{{\∂f}_n}{{\∂x}_1}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \frac{{\∂f}_1}{{\∂x}_m}& \frac{{\∂f}_2}{{\∂x}_m}& ...& \frac{{\∂f}_n}{{\∂x}_m}\end{array}\right]$

则可得到：

为了使

沿着其梯度方向减小，必有

的法方程等于零：

将f(x)在x^b(b为迭代次数)处作泰勒公式线性展开得：

$f\left(x\right)=f\left(x^b\right)+{\left[f^{\′}\left(x^b\right)\right]}^T(x-x^b)+{\left[\frac{f^{\′\′}\left(x^b\right)}{2!}\right]}^T{(x-x^b)}^2+\·\·\·{\left[\frac{f^n\left(x^b\right)}{n!}\right]}^T{(x-x^b)}^n+R_n\left(x\right)$

又[f″(x^b)]^T(x-x^b)²以后各项的值都远小于前两项，故近似有：

f(x)≈f(x^b)+Df(x^b)^T(x-x^b)    (9)

令l_b(x)=Df(x^b)^T(x-x^b)+f(x^b)，代入法方程（8）得到：

0=Df(x)^Tf(x)≈Df(x^b)^Tl_b(x)

　　=Df(x^b)^T[Df(x^b)(x-x^b)+f(x^b)]

　　=Df(x^b)^TDf(x^b)(x-x^b)+Df(x^b)^Tf(x^b)

则迭代公式为：

$x^{b+1}=x^b-\frac{\mathrm{Df}{\left(x^b\right)}^{T}f\left(x^b\right)}{\mathrm{Df}{\left(x^b\right)}^T\mathrm{Df}\left(x^b\right)}---\left(10\right)$

为表示方便，上式可进一步写为：

式中：G(x^b)=Df(x^b)^TDf(x^b)，g(x^b)=Df(x^b)^Tf(x^b)，p(x^b)=-G(x^b)^-1g(x^b)；

3.2）：考虑阻尼时有：

其中μ_b为第b次迭代的阻尼因子，I是与G(x^b)同维数的单位对角阵；则最终迭代公式为：

当μ_b>0总可以保证

是收敛的；

3.3）计算得出f(x^k),Df(x^k),

令0→j（j为3.7步骤所需的一个判断变量）；

3.4）解方程组

求得

3.5）计算x^b+1=x^b+p(x^b)及

3.6）若

则取

（缩放因子v>1，可取2,5,10），转3.3），否则转3.7）；

3.7）若j=0，则取μ_b=vμ_b,1→j，转3.3），否则转3.8）；

3.8）若||p(x^b)||≤ε₂，则x^b+1为极小值点x^*的近似，算法停止，输出x^b+1；否则取x^b+1→x^b，b+1→b，转3.3）。

本发明是通过应用AR模型谱估计计算谐波初始参数来建立信号模型的方法：将待分析的采样信号序列进行AR模型谱估计，给出精度较低的谐波参数(频率、幅值和相位)信息；然后基于初步的谐波信息建立信号的参数模型，再取模型进行式（1）计算。模型(1)为非线性模型，需采用非线性优化算法；同时，为减小模型的非线性，将参数模型进行处理，将幅值、相位(A_m,θ_m)等价转化为求解两个幅值(C_m,D_m)进行计算。并且将式中参数fr_m、C_m和D_m统一用向量x表示，则对应于各采样时刻的误差向量用f(x)表示，进而将非线性优化问题转化最小值问题进行求解。

其中目标函数

表征了信号计算结果和实际采样结果的误差，差别越小表明参数计算精度越高。

所述参数模型式（1）中的m次数可根据功率谱中局部最大值个数确定，然后依次确定功率谱中每一个局部最大值处对应的频率和局部最大值的开方，分别作为参数fr_m、A_m的初值，然后用0～2π内的随机值作为相位参数θ_m的初值。

由于谐波的出现使待分析信号波形产生了畸变，所以通过式(1)中参数的确定可通过选取AR功率谱中的局部最大值实现，然后用0～2π内的随机值作为相位参数θ_m的迭代初值；

所述参数模型式（1）中的m次数可根据功率谱中局部最大值个数确定，然后依次确定功率谱中每一个局部最大值处对应的频率和局部最大值的开方，把它们分别作为参数fr_m、A_m的初值，然后用0～2π内的随机值作为相位参数θ_m的初值。

所述步骤三中的3.2的阻尼因子μ_k，当μ_b>0总可以保证

因而是收敛的，当μ_b太大时会收敛速度下降，若μ_b太小则收敛域过小，原则上μ_b=10^-4～10^-2，通常取μ_b=10^-2；

本发明的显著效果是：

利用AR谱估计给出模型的结构和初始参数，解决了非线性优化算法的参数建模问题，也给出了参数的初值；结合非线性优化算法，克服了AR谱估计不能计算谐波幅值和相位的缺点，且进一步提高了频率的计算精度。

非线性优化算法的设计中，首先选用共轭梯度法对参数初值进行预处理，使之接近全局最优解；随后再使用阻尼最小二乘法能使参数迭代的谐波参数计算值)，降低迭代算法对初值的敏感度，提高了迭代稳定性和计算效率；

AR谱估计与非线性优化结合方法的参数计算精度高出加Hanning窗插值法1～4个数量级，且有一定的抗噪能力。

附图说明

图1为本发明流程图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明：

如图1所示，本发明实现的流程图按以下步骤进行：

（1）输入采样数据y(n)并作AR谱估计运算给出精度较低的谐波参数信息；

（2）基于初步谐波信息，通过选取功率谱局部最大值建立参数模型；

（3）应用非线性优化算法确定准确的谐波参数，步骤包括：

1）置初始近似值x⁰∈Rⁿ，误差限ε₁>0；

2）计算

令0→k；

3）求最优步长t_k满足

计算迭代值x^k+1=x^k+λ_kp_k以及

4）计算

若

转步骤6），否则，转步骤5）；

5）计算

k+1→k，转步骤3；

6）置ε₂>0，阻尼因子μ₀=10^-2，缩放因子v>1（可取2,5,10），0→k；

7）计算f(x^k),Df(x^k),

0→j；

8）解方程组

求得p(x^k)=p(μ_k)；

9）计算x^k+1=x^k+p(x^k)及

10）若

则取

转步骤7），否则转步骤11）；

11）若j=0，则取μ_k=vμ_k,1→j，转步骤7），否则转步骤12）；

12）若||p(x^k)||≤ε₂，则x^k+1为极小值点x^*的近似，算法停止，输出x^k+1；否则取x^k+1→x^k，k+1→k，转步骤7）；

算例1：取含基波与相近频率间谐波的电网信号模型如下：

$y\left(t\right)=\sin (2\mathrm{\π}50t+150*\frac{\mathrm{\π}}{180})+0.01\sin (2\mathrm{\π}55t+56*\frac{\mathrm{\π}}{180})+v\left(t\right)$

其中，v(t)是信噪比为60dB的高斯白噪声；取采样频率f_s=1kHz，采样点数N=512个，AR模型的阶数p=176；在非线性优化法中，取误差限ε₁=1e-2，ε₂=1e-6，缩放因子v=2，阻尼因子μ₀=10^-2。

表1谐波与间谐波参数计算精度相对误差表

注：加窗指加Hanning窗插值FFT算法，AR指AR谱估计与非线性优化结合方法。

算例2：取含多个谐波与间谐波的电网信号模型如下：

$y\left(t\right)=\sin (2\mathrm{\π}50t+150*\frac{\mathrm{\π}}{180})+0.01\sin (2\mathrm{\π}120.2t+56*\frac{\mathrm{\π}}{180})+0.01\sin (2\mathrm{\π}150t+120*\frac{\mathrm{\π}}{180})$

$+0.015\sin (2\mathrm{\π}217.9t+85*\frac{\mathrm{\π}}{180})+0.0012\sin (2\mathrm{\π}245t+23*\frac{\mathrm{\π}}{180})+v\left(t\right)$

其中，v(t)是信噪比为60dB的高斯白噪声；取采样频率f_s=1kHz，采样点数N=512个，AR模型的阶数p=176；在非线性优化法中，取误差限ε₁=1e-2，ε₂=1e-6，缩放因子v=2，阻尼因子μ₀=10^-2。

表2谐波与间谐波参数计算精度相对误差表

注：加窗法指加Hanning窗插值FFT算法，结合法指AR谱估计与非线性优化结合方法。

从表1和2的结果可见，噪声下的谐波和间谐波的检测中，AR谱估计和非线性优化的结合方法能准确地检测出信号的谐波与间谐波分量，且参数计算精度较加Hanning窗插值法提高了约1～4个数量级。

本发明的技术方案可应用于电网谐波分析、电能计量和电能质量监测。

Claims (3)

1.一种结合AR谱估计与非线性优化的间谐波检测方法，其特征在于,按以下步骤进行：

步骤一：AR谱初步估计

1.1）：根据谐波参数的初步信息建立该谐波的参数模型，取参数模型如下：

$y\left(t\right)=\underset{m}{\mathrm{\Σ}}{A}_m\sin (2\mathrm{\π}{\mathrm{fr}}_{m}t+{\mathrm{\θ}}_m)---\left(1\right)$

式中y(t)为待分析信号，A_m为各次谐波和间谐波的幅值，fr_m为对应的频率，θ_m为对应的相位，m为对应的谐波个数，t为时间；

1.2）：对求解信号模型参数进行非线性优化算法：将所述模型（1）作如下处理，将幅值、相位(A_m,θ_m)等价转化为求解两个幅值(C_m,D_m)：

A_msin(2πfr_mt+θ_m)=A_msin(2πfr_mt)cos(θ_m)+A_mcos(2πfr_mt)sin(θ_m)

　　　　　　　　=C_msin(2πfr_mt)+D_mcos(2πfr_mt)    (2)

其中，C_m=A_mcos(θ_m)，D_m=A_msin(θ_m)；

1.3）：针对步骤一中的参数模型（1），将式(2)代入式(1)，建立最小二乘方程如下：

式中n为采样点数；目标函数(x₁,x₂,...x_m)表征了信号计算结果和实际采样结果的误差，令y(t_n)为信号第n点的计算值、y₁(t_n)为信号第n点实际采样值，将参数f_m、C_m和D_m统一用向量x表示，对应于各采样时刻的误差向量用f(x)表示：

f(x)=[f₁,f₂，...，f_n]

　　　=[y(t₁)-y₁(t₁),(t₂)-y₁(t₂)，...，(t_n)-y₁(t_n)]

代入（3）式可得：

式中[f(x)]^T表示f(x)的转置，则参数的求解转化为求如下最值的问题：

步骤二：共轭梯度法一次迭代，也就是参数初值进行预先处理，使参数迭代值接近全局最优解：

2.1）对参数初值进行迭代，迭代方程为：

x^k+1=x^k+λ_kp_k    (5)

其中，k为迭代次数，λ_k为最优步长，p_k为迭代方向，

2.2）p_k满足下式（7）

置初始近似值x⁰∈Rⁿ（Rⁿ为实数域），一次迭代误差限ε₁>0

2.3）：计算

令0→k，p₀是迭代方向矩阵p_k的初值；

2.4）：求最优步长λ_k满足计算迭代值x^k+1=x^k+λ_kp_k以及

2.5）：计算当迭代到时，转步骤三，否则，转2.6；

2.6）：计算

k+1→k，转2.4；

步骤三：当

后，再用阻尼最小二乘法二次迭代：即采用阻尼最小二乘法快速寻找到全局最优解，即最终的间谐波参数计算值：

3.1）：为使(x)沿着其梯度方向减小，必有

(x)的法方程等于零：

可得：

式中：G(x^b)=Df(x^b)^TDf(x^b)，g(x^b)=Df(x^b)^Tf(x^b)，p(x^b)=-G(x^b)^-1g(x^b)；

3.2）：考虑阻尼时有：

其中μ_b为第b次迭代的阻尼因子，I是与G(x^b)同维数的单位对角阵；则最终迭代公式为：

当μ_b>0总可以保证

是收敛的，原则上μ_b=10^-4～10^-2；

3.3）：计算得出f(x^b),Df(x^b),(x^b),

(x^b),G(x^b)，令0→j（j为3.7步骤所需的一个判断变量）；

3.4）：解方程组求得

3.5）：计算x^b+1=x^b+p(x^b)及

3.6）：若

则取

1→j，（缩放因子v>1，可取2,5,10），转3.3），否则转3.7）；

3.7）：若j=0，则取μ_b=vμ_b,1→j，转3.3），否则转3.8）；

3.8）：若||p(x^b)||≤ε₂，则x^b+1为极小值点x^*的近似，算法停止，输出x^b+1；否则取x^b+1→x^b，b+1→b，转3.3）。

2.根据权利要求1所述一种结合AR谱估计与非线性优化的间谐波检测方法，其特征在于：所述步骤一中的参数模型式（1）中的m次数可根据功率谱中局部最大值个数确定，然后依次确定功率谱中每一个局部最大值处对应的频率和局部最大值的开方，把它们分别作为参数f_m、A_m的初值，然后用0～2π内的随机值作为相位参数θ_m的初值。

3.根据权利要求1所述一种结合AR谱估计与非线性优化的间谐波检测方法，其特征在于：所述步骤三中的3.2）阻尼因子μ_b，当μ_b>0总可以保证

因而是收敛的，当μ_b太大时会收敛速度下降，若μ_b太小则收敛域过小，原则上μ_b=10^-4～10^-2，通常取μ_b=10^-2。

Images