一种基于滑动窗DFT的高精度谐波参数估计方法与流程

本发明涉及一种针对平稳周期信号的谐波参数估计方法，具体指一种基于滑动窗傅里叶变换(swdft)和prony算法的谐波参数高精度估计方法，本方法用于估计谐波频率、幅度和相位参数，属于电力系统信号谐波分析技术领域。

背景技术：

随着电力电子技术的发展以及工业生产的扩大，电力系统中电子元器件及非线性负荷日益增多。非线性负荷在带来巨大经济效益的同时也向电网中注入了大量的高次谐波，谐波污染变得越来越严重，在加剧电信号畸变程度的同时，给电力系统以及用户造成了诸多危害，如电机过热、电力系统继电保护异常、输电线路损耗增加等。目前关于谐波的治理大多采取向系统注入相应谐波成分反向分量的措施和安装滤波器等方法，因此，对电网信号进行谐波分析，对电能计量及电能质量分析与治理具有重大的研究意义和实用价值。

离散傅里叶变换(dft)是电力系统中进行谐波分析最常用的方法，该类算法具有运算简单、计算效率高的特点，对于平稳无噪声的信号有较好地估计结果。但其存在频率混叠、频谱泄露和栅栏效应等固有缺点，影响了谐波分析的精度。而加窗、插值等改进措施能够在一定程度上降低频谱泄露和栅栏效应造成的估计误差，但当各谐波分量的频率间隔较近时，此类方法并不能有效工作。为了提高估计精度，迭代dft算法应运而生，但迭代的引入会产生附加的计算量，造成巨大的计算负担。除了dft之外，还有一些基于时域的谐波分析方法，如卡尔曼滤波方法，然而这些方法的实现需要建立在已知一些待处理信号数据特征的基础上，且状态矩阵也需要事先准确定义，否则将达不到预先的估计效果。

技术实现要素：

发明目的：针对上述现有技术，提供一种基于滑动窗dft的高精度谐波参数估计方法，可以大大提高电网信号频率、幅值和相位等参数的估计精度。

技术方案：一种基于滑动窗dft的高精度谐波参数估计方法，包括以下步骤：

步骤1：采集电力系统中的多频信号，采样后得到时刻n的离散信号x(n)满足其中am，ωm，φm分别表示第m次谐波的幅值、数字频率和相位，m为信号所含谐波的个数，s(n)为多频谐波信号部分，q(n)为加性高斯白噪声部分，观测的时域数据范围选取n＝0,1,…,l-1，l为离散点个数；

由余弦和指数之间的变换关系将多频谐波信号部分改写后，得到指数信号模型并以此指数模型进行后续分析与估计；其中，ai、ωi和φi为指数信号模型中第i个频率分量的幅值、数字频率和相位，p为多频指数分量个数；

步骤2：对具有指数形式的时域采样信号x(n)做n点滑动窗傅里叶变换，得到序列x(m)；

步骤3：对变换后的序列x(m)使用复最小二乘准则扩展后的prony算法估计模型中各分量的频率；

步骤4：在频率估计值的基础上利用扩展的prony算法估计指数信号模型各次分量对应的幅值和相位；

步骤5：由信号的指数模型估计值，恢复初始电力系统信号的谐波成分参数。

进一步的，所述步骤2中利用滑动窗傅里叶变换得到序列x(m)包括如下步骤：

2.1)将各分量的频率表示为ωi＝2π/n(ki+δi)，其中ki∈{-n/2,…,0,…,n/2}，|δi|≤0.5分别代表各未知频率以2π/n规格化的整数和小数部分，n为滑动窗傅里叶变换的窗长度，通过频率粗估计得出各次分量对应的谱线值ki；

2.2)对待分析的信号分别计算ki次谱线处的滑动窗傅里叶变换结果，得到序列

2.3)将p个频率分量处的变换结果相加，得到所需的序列x(m)。

进一步的，所述步骤3中采用复最小二乘准则扩展后的prony算法进行频率估计的步骤为：

3.1)根据prony算法采用的估计模型n＝0,1,…,l-1，滑动窗傅里叶变换序列仍然满足s(m)＝a1s(m-1)+a2s(m-2)+…+aps(m-p)，从而得到估计参数所使用特征方程a1z^-1+a2z^-2+…+apz^-p-1＝0，其特征根即包含了要估计的频率参数，考虑噪声干扰时，线性关系即为x(m)＝a1x(m-1)+a2x(m-2)+…+apx(m-p)；其中，ci＝aie^jφi是包含时域信号幅度和相位信息的参数；zi为方程的特征根；ai,i＝1,2,…,p表示差分方程的系数；s(m)为时域纯净多频信号s(n)的在每个频率分量对应谱线ki处滑动窗傅里叶变换结果的叠加序列；x(m)为考虑噪声干扰时，含噪多频信号x(n)在p个频率对应谱线处的swdft变换进行叠加的结果；

3.2)定义向量xm＝[x(m),x(m+1),…,x(m+l-n)]^t，则xm＝a1xm-1+a2xm-2+…+apxm-p，相应的估计误差为根据复最小二乘准则，搜索总体均方误差最小值，此时估计得到特征多项式的系数i＝1,2,…,p；

3.3)将方程系数的估计值代入求解信号的特征多项式，由方程的特征根zi得到各信号分量的频率估计值

进一步的，所示步骤4中估计各分量对应的幅值和相角，步骤如下：

4.1)定义包含各频率分量幅值和相位参数信息的参数为由prony算法的估计模型代入特征根和采样的时域信号，获得复最小二乘准则下的估计值其中，和分别为s(m)中各分量的幅值信息和相位参数，且与变换前的信号分量相对应；

4.2)根据估计值结合滑动窗傅里叶对各谐波成分幅值和相位产生影响的数学表达式，计算得出指数信号模型中各频率分量的幅值和相位。

进一步的，步骤5中对所估计的指数模型参数，对其保留频率估计值为正的分量，并将对应的幅度估计值乘以2，相位估计值维持不变，得到电力系统信号原始的谐波成分参数。

有益效果：本发明中，将含有谐波污染的电力系统信号的参数估计问题通过余弦和指数之间的数学转换关系，导出了信号的多频指数模型，并以此为基础进行本方案的具体谐波参数估计过程。进一步而言，将信号看作若干个指数分量的叠加，即信号分量的个数为实际信号的两倍，其中一半对应的频率为各次谐波频率的负数值。滑动窗dft(swdft)保留了信号指数模型的线性关系，且抑制了噪声的干扰,增强了各分量的信噪比(snr)，对变换后的序列采用prony模型进行参数估计，则信噪比的提升可以带来估计精度的提高。此外，为了增强算法的鲁棒性和估计性能，本发明中采用的prony算法是建立在复最小二乘(cls)准则之上的。因此，基于swdft和prony算法的谐波参数估计方法可以获得较高的估计精度，且算法的鲁棒抗噪性能和计算复杂度均具有明显的优势。

与现有的技术相比，本发明具有以下优点：1.充分利用了swdft变换结果和原始时域序列的数学关系及其带来的信噪比增益，提升了参数估计方法的抗噪性能和估计精度。2.与传统的prony算法相比，采用cls准则的扩展prony算法在降低方程求解误差的同时，增强了谐波参数估计方法的鲁棒性。

附图说明

图1为电网谐波信号频率估计均方误差随信噪比的变化图；

图2为电网信号基频分量幅度和相位估计误差随信噪比的变化图，其中图2(a)为幅度估计误差曲线，图2(b)为相位估计误差曲线；

图3为电网信号3次谐波成分幅度和相位估计误差随信噪比的变化图，其中图3(a)为幅度估计误差曲线，图3(b)为相位估计误差曲线；

图4为电网信号5次谐波成分幅度和相位估计误差随信噪比的变化图，其中图4(a)为幅度估计误差曲线，图4(b)为相位估计误差曲线。

具体实施方式

下面结合附图对本发明做更进一步的解释。

一种基于滑动窗傅里叶变换(swdft)的高精度谐波参数估计方法，噪声背景下电力系统的多频信号为：

式中，s(t)表示时间t时包含基频成分的纯净谐波信号部分，q(t)代表均值为0方差为的加性白噪声，m为谐波的个数，fm、am和φm分别表示第m个谐波对应的频率、幅值以及相位。

以周期ts对该信号进行离散采样，得到l点离散序列x(n)：

式中ωm＝2πfmts，s(n)为时刻n时模型中纯净的多频谐波信号，q(n)为模型中的加性高斯白噪声部分。由于每一个正弦分量可以写成两个频率互为相反数的复指数信号相加，(2)式的信号模型可以改写为

其中ai、ωi和φi为指数信号模型中第i个频率分量的幅值、数字频率和相位，p为(3)式所示模型中多频指数分量的个数。与(2)式比较可以发现，指数模型中的信号分量变为了原来的两倍，即p＝2m，产生了相同的负频率分量，相应的幅值变为原来的1/2，而信号的相位维持不变。因此，通过对(3)式所示模型进行估计，然后保留频率估计值为正的部分分量，并将对应的幅值全部扩大为估计值的两倍，则可以得到初始电力系统信号的谐波估计参数。

注意到(3)式所示信号中的任意分量的频率均可以表示为ωi＝2π/n(ki+δi)，其中ki∈{-n/2,…,0,…,n/2}，|δi|≤0.5为各未知频率以2π/n规格化的整数和小数部分，且n为可选参量，对应于后续步骤中swdft变换的滑动窗长度。上述噪声背景下的多频信号，多频指数分量的个数p和规格化的整数数值ki,i＝1,…,p可以由fft谱峰搜索、奇异值分解(svd)定阶等频率粗估计方案准确获得。

对于模型中纯净的多频谐波信号部分s(n)，其n点swdft的k次谱线处的变换结果为：

其中m表示滑动窗的起始时刻，距信号初始时刻延迟m个样点；s(n+m)即为n+m时刻对应的多频信号成分，且式中

由(4)式，swdft变换的结果仍然满足(3)式的指数模型，且swdft变换仅改变了各分量的幅度和相位。根据prony算法，其采用的数学模型为一组p个具有任意幅值、相位、频率和衰减因子的指数函数，离散时间的函数形式为

我们使用作为s(n)的近似，式中方程的特征根αi和ts分别为衰减因子与采样周期。

prony算法的关键在于推导发现(6)式的拟合是一个常系数线性差分方程的齐次解，即满足递推的线性预测差分方程式

s(n)＝a1s(n-1)+a2s(n-2)+…+aps(n-p)(7)

ai,i＝1,2,…,p表示差分方程的系数或线性预测系数，相应的特征多项式为

因此，由(4)式，变换后的序列sk(m)仍然符合prony算法的线性预测结构，即sk(m)＝a1sk(m-1)+…+apsk(m-p)，我们可以使用prony算法估计各分量的参数。具体而言，prony算法通过观测的时间序列，由矩阵关系求出特征多项式的系数，并以此计算特征根，得到信号的频率参数，然后再以特征多项式的解为基础得到信号的幅度和相位估计值。

结合(5)式发现，swdft变换之后每个频率分量的幅度随着其频率与k次谱线对应角频率2πk/n之间的距离的增大而产生衰减。对于第i个频率分量，频率规格化后的整数值为ki，在k＝ki时幅度衰减最小，即纯净多频信号s(n)在ki次谱线处的swdft结果能够使信号中第i个频率分量的幅度衰减最小。

因此，若考虑噪声的干扰，对于实际的含噪多频电网信号x(n)，在第i个频率分量对应的谱线k＝ki处，其n点swdft的变换结果为其中m表示滑动窗相对初始信号时刻的延迟样点数，x(n+m)为n+m时刻对应的电网信号。根据(4)和(5)式，第i个频率分量的信噪比也将较变换前发生变化：在swdft之前的时域信号中，该分量的信噪比为变换后噪声部分q(n)的功率变为了因此变换后的信噪比为

如(9)式所示，对于含有噪声干扰的信号，其ki次谱线处的swdft结果增强了角频率2πki/n附近频率分量的信噪比，相应的参数估计性能也会因此得到提升。如果要精确地估计出各分量的频率参数，我们则需要对信号做p次swdft计算和估计，而这将带来巨大的计算负担。为了在保证估计精度的同时尽可能提高计算效率，接下来考察p路谱线变换结果的相加情形，显然纯净谐波部分的序列仍然满足线性关系：

s(m)＝a1s(m-1)+a2s(m-2)+…+aps(m-p)(10)

s(m)为时域多频信号s(n)的在每个频率分量对应谱线ki处滑动窗傅里叶变换结果的叠加序列。

为了分析噪声干扰环境中该序列swdft前后的信噪比变化，将s(m)详细展开：

其中kj表示信号中第j个信号分量频率的规格化整数值。

因此，对考虑噪声干扰下的变换序列而言，为待分析信号在ki次谱线处的滑动窗傅里叶变换结果，第i次频率分量的信噪比相比原始时域信号变为

通过数学推导可以发现关系式snri≤snri,dft≤snri,swdft，该不等式表明使用序列x(m)进行参数估计的性能将稍低于使用而显著优于使用原始的时域采样信号。因此，本发明中的算法将采用各次信号分量对应谱线处变换结果叠加的序列x(m)进行参数估计。

此外，为了增强算法的鲁棒性，降低prony算法中求解矩阵方程时的误差，采用复最小二乘(cls)结构进行相应参数的估计，具体过程如下：

(10)式在噪声情形中可扩展得到

定义向量xm＝[x(m),x(m+1),…,x(m+l-n)]^t，上式可写为

xm＝a1xm-1+a2xm-2+…+apxm-p(15)

相应的估计误差为

则其自共轭矩阵为cls结构求解(8)式所示特征方程的目标是最小化估计的总体均方误差

令则

通过求解式(17)计算出特征多项式的系数，然后代入特征方程解出p个特征根，相应的频率估计值为

式中∠(·)表示求角度运算。注意到swdft变换序列x(m)因引入n点长的滑动窗变换，序列相比时域信号x(n)，长度由l变为l-n+1，因此由prony算法，包含幅度和相位信息的参数满足

类似地，估计幅值和相角时仍然采用cls准则。令x＝[x(0),x(1),…,x(l-n)]^t，和b＝[b1,b2,…,bp]^t，则采用cls结构的估计方程为

zi为含有频率信息的特征根向量，b为包含幅度和相位信息的参数所扩展成的向量。

由(21)式计算出估计值结合(11)式即可得到信号指数模型对应的幅值和相位

为了恢复出初始电力系统信号的谐波成分参数，只需选取上述估计结果中频率为正值的信号分量，对应的谐波幅度修正为幅度估计值的两倍，相位则维持不变。

以下结合一个实施算例的仿真结果对本发明做进一步阐述。

取谐波电网信号模型为：

x(t)＝20cos(2πf0t)+4cos(3×2πf0t+π/2)+3cos(5×2πf0t+π/3)+q(t)

其中，基础频率为f0＝49.8hz，q(t)表示信号的干扰噪声，信号中除了基频分量外，还包含3次和5次谐波。取采样频率fs＝800hz，采样点数l＝250，swdft变换的窗长度选择n＝round(fs/f)＝16，round(·)代表四舍五入取整。

在第一组模拟中，我们利用本发明提出的算法对电网信号的谐波频率进行估计，为了明确提出算法的优势，我们将此算法与直接利用时域信号进行prony算法(cls准则下)估计的频率估计结果进行比较，并采用均方误差(mse)作为性能衡量数值，均方误差越大则说明算法的估计性能越差。图1画出了上述两种算法的频率估计误差随信噪比的变化曲线。如图所示，相同颜色和图形标注的实线和虚线代表分别代表本发明公开的算法和时域prony算法的频率估计均方误差，仿真表明，本发明所提出算法能够精确地估计出谐波分量的频率参数，其中对基础频率的估计结果尤为准确，该方案的估计误差显著低于时域prony算法。

在频率参数的基础上，继续利用prony算法对各频率分量的幅度和相位进行估计，基频、3次谐波和5次谐波成分的幅度和相位估计误差分别如图2-4所示。对比每张仿真图中两种算法对同一信号分量进行估计的均方误差发现，所提出方案在谐波成分幅值和相位估计方面的估计误差均小于时域prony算法。综上所述，本发明所提出的谐波估计算法在参数估计方面性能优异，能够准确地给出各次谐波成分频率、幅度和相角的估计结果。

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

最后要说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施算例的仿真结果对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。