CN103308804A - 基于快速k-s变换电能质量扰动信号时频参数提取方法 - Google Patents

Abstract

基于快速K-S变换电能质量扰动信号时频参数提取方法，包括以下步骤: （1）被测信号采样；（2）快速K-S变换；（3）重复步骤（2），直到分析出所有扰动信号的幅值、相位和起止、突变时刻；通过以上步骤完成电能质量扰动信号时频参数的提取。与现有技术相比，本发明采用Kaiser取代传统S变换的Gauss窗，且能实现随频率变化的Kaiser窗函数宽度、高度的自适应调节，能提高传统S变换时频分析的能量聚集性和频率自适应调节能力；通过特征频率判断，能减少计算量；适用范围广泛，快速K-S变换与电能质量扰动信号时频参数提取方法在振动信号分析、故障诊断、模式识别、噪声测量、生物医学检测领域具有广阔的应用前景。

Description

基于快速K-S变换电能质量扰动信号时频参数提取方法

技术领域

本发明涉及电能质量信号分析领域，尤其是涉及一种基于快速K-S变换电能质量扰动信号时频参数提取方法。

背景技术

根据电压扰动的频谱特征、持续时间、幅值变化等，电力系统中各种扰动引起的电能质量问题主要分为稳态事件和暂态事件两大类。稳态电能质量问题主要包括谐波、电压波动与闪变、频率偏差等；暂态电能质量问题主要指电力系统发生故障或较大负荷变化引起的电压均方根值在短时间内随时间改变的现象，包括电压骤升、电压骤降、电压间断、暂态震荡等。

电能质量扰动信号时频参数包括频率、幅值、相角、起止时刻、突变时刻等。电能质量扰动时频分析是电能质量问题的主要研究内容，是计算电能质量参数、研究电能质量问题产生根源、改善和提高电能质量的前提和关键依据，对保证电网和电气设备安全、经济运行，增强国民经济的总体效益和保障居民正常生活有着重要的意义。

目前，电能质量扰动信号分析的主要方法有快速傅里叶变换（FFT）、短时傅里叶变换、小波变换、神经网络。

FFT只能了解信号的全局特性，不能反映信号频率随时间的变化规律，不具有时频局部性，无法全面描述扰动信号的局部时变特征。

短时傅里叶变换虽然在一定程度上克服了FFT不具有局部分析能力的缺陷，可获得某一时刻的频率信息和某一频率点的信号幅值信息，但其窗函数固定，无法自适应调节信号分析的时域、频域分辨率。

小波变换是具有多分辨率特性的时频局部化分析方法，其缺点是不同尺度的小波函数在频域中存在相互干扰。

神经网络需要进行网络训练，且需要大量先验信息，算法复杂度高，无法满足实时性电能质量扰动信号分析需要。

S变换是由小波变换和短时傅里叶结合发展起来的一种新型时频分析方法，以FFT→Gauss窗→IFFT为算法思路建立时频分析变换，克服了短时傅里叶变换窗口高度和宽度固定的缺陷，计算复杂度低于小波变换和神经网络，大大提高扰动信号时频分析精度和实用性。但S变换其Gauss窗形状调节因子无统一标准，且自适应调节能力较差，对复杂畸变扰动信号难以实现较高的时频能量聚集，且时频复数矩阵信息大，运算复杂，实施性较差，因此，S变换无法满足扰动信号时频特征的准确分析的要求。

以上分析可知，现有的时频分析方法频域或时域信息单一；分辨率低，不能精细刻画信号的时频信息；计算复杂，实用性差，在全面性、实时性、准确性方面不能兼顾，工程应用范围下。

发明内容

本发明所要解决的技术问题是，克服现有技术存在的上述缺陷，提供一种能分析频域和时域信息，分辨率高，能精细刻画信号时频信息，计算简单，实用性好，且能兼顾全面性、实时性、准确性的基于快速K-S变换的电能质量扰动信号时频参数提取方法。

本发明解决其技术问题所采用的技术方案是，基于快速K-S变换电能质量扰动信号时频参数提取方法，包括以下步骤:

（1）被测信号采样：将电能质量扰动信号输入电力系统分析仪内进行采样，得被测信号x(n)；

（2）快速K-S（凯撒-S）变换：对步骤（1）采样好的被测信号进行K-S变换，得到特征频率f_m对应的扰动信号K-S变换谱，分析出特征频率为f_m对应的扰动信号的幅值、相位和起止、突变时刻；

（3）重复步骤（2），直到计算出所有特征频率f_m对应的扰动信号K-S变换谱，完成扰动信号的K-S变换KS(τ,f)，分析出所有扰动信号的幅值、相位和起止、突变时刻；

通过以上步骤完成电能质量扰动信号时频参数的提取。

进一步，所述步骤（2）中，K-S变换的方法包括以下步骤：

（1）傅里叶变换运算：对采样好的被测信号x(n)长度N点进行傅里叶变换运算，得到信号频谱y(n);

（2）特征频率提取:从步骤（1）中计算出的信号频谱y(n)中提取扰动信号特征频率f_m；

（3）自适应Kaiser窗：对步骤（2）提取的扰动信号特征频率f_m ,计算自适应优化Kaiser窗傅里叶变换K(f_m,ω_n)，ω_n代表为频域变量；

（4）频谱平移和频谱相乘：平移信号频谱y(n)到y(n+m)，计算y(m+n)K(f_m,ω_n)， y(n)为信号频谱，m为平移大小，K(f_m,ω_n)为特征频率f_m对应的Kaiser窗傅里叶变换；

（5）傅里叶反变换运算：计算y(m+n)K(f_m,ω_n)傅里叶反变换，得到扰动信号特征频率f_m对应的K-S变换谱ks(f_m,t_n)，t_n对应各采样时刻点；

（6）复数时频矩阵分析：根据K-S变换谱，进行K-S矩阵元素分析，由复数时频矩阵分离得到时频相位矩阵和时频幅值矩阵分析，通过对时频相位矩阵和时频幅值矩阵分析，得到特征频率为f_m的对应的扰动信号的幅值、相位和起止、突变时刻。

进一步，所述步骤（3）中，Kaiser窗为一组窗度和高度能调的窗函数，且能自由选择主瓣宽度和旁瓣高度之间的比重，其时域离散表达式为：

$k\left(n\right)=\frac{I_0\left[\mathrm{\β}\sqrt{1-(1-\frac{2n}{N-1})}\right]}{I_0\left(\mathrm{\β}\right)}$ ；

式中，β为形状参数，n为时域变量，N为窗函数长度 ，I₀为第一类修正零阶Bessel函数，其幂级数展开式：

$I_0\left(\mathrm{\β}\right)=1+\underset{m=1}{\overset{\∞}{\mathrm{\Σ}}}{\left[\frac{{(\mathrm{\β}/2)}^m}{m!}\right]}^2$ ；

式中，β为形状参数，用以调节主瓣宽度和旁瓣电平，m取正整数，Kaiser窗频域表达式为：

$K\left(\mathrm{\ω}\right)=\frac{N}{I_0\left(\mathrm{\β}\right)}\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{(\mathrm{N\ω}/2)}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{(\mathrm{N\ω}/2)}^2}}$ ；

式中，sinh为双曲正弦函数，N为窗函数长度，β为形状参数，ω为频域变量；

基于Kaser窗，定义自适应优化Kaiser窗形状参数为:

$\mathrm{\β}=\frac{\mathrm{\α}}{f}$ ；

式中，α为常数，f为被测信号频率，由上式可见，Kaiser窗口的高度和宽度随被测信号频率变化而变化，其幅值与频率大小成正比，宽度随频率的增大而变小，即能实现窗口宽度和高度随信号频率的自动调节。

进一步，所述K-S变换的处理过程为：

将信号x(t)的连续K-S变换KS(τ,f)定义如下：

$\mathrm{KS}(\mathrm{\τ},f)={\∫}_{-\∞}^{\∞}x\left(t\right)k\left(\frac{\mathrm{\α}(\mathrm{\τ}-t)}{f}\right)e^{-i2\mathrm{\πft}}\mathrm{dt}$ ；

式中，k()为Kaiser窗函数，τ,t为时移因子，f为频率，α为常数，上式为卷积方程，利用卷积定理，将时域卷积运算简化为频域乘积运算；

令长度为N的连续采样离散采样信号x(n)：

x(n)=[x₁…x_n…x_N]；

式中，n对应时域序列，N为采样长度；

采样信号x的傅里叶频谱表示为：

y(n)=[y₁…y_n…y_N]；

式中，n对应频域序列，N为频谱长度；

式中，y_n满足：

$y_n=\frac{2}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i{e}^{-j\frac{2\mathrm{\πi}}{N}}$ ；

式中，N为采样长度， j代表虚数单位，x_i为采样序列，i为正整数；

特征频率f_m对应的谱线y_n满足以下公式：

y_n-1＜y_n＞y_n+1；

|y_n|＞ε；

式中，ε为设定阈值，n为对应频域序列，谱线y_n对应频率f_m为特征频率；

将信号频谱y(n)平移首尾相连构建频谱矩阵H_M×N：

$H_{M\·N}=\left[\begin{array}{ccccc}{y}_2& y_1& ...& y_N& y_1\\ y_1& y_4& ...& y_1& y_2\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ y_M& y_{M+1}& ...& y_{M-2}& y_{M-1}\\ y_{M+1}& y_{M+2}& ...& y_{M-2}& y_M\end{array}\right]$ ；

式中，N为频谱长度，M频谱平移次数，y对应频谱各谱线；

构建不同频率下Kaiser窗函数，记为Kaiser窗函数矩阵，矩阵各行对应各频率下自适应优化Kaiser窗函数：

$k_{M\×N}=\left[\begin{array}{ccccc}k(f_1,t_1)& ...& k(f_1,t_n)& ...& k(f_1,t_N)\\ .& & .& & .\\ .& ...& .& ...& .\\ .& & .& & .\\ k(f_m,t_1)& ...& k(f_m,t_n)& ...& k(f_m,t_N)\\ .& & .& & .\\ .& ...& .& ...& .\\ .& & .& & .\\ k(f_M,t_1)& ...& k(f_M,t_n)& ...& k(f_M,t_N)\end{array}\right]$ ；

式中，N为频谱长度， M频谱平移次数；m对应矩阵行，n对应矩阵列，t对应采样时刻点；

式中：

$k(f_m,t_n)=\frac{I_0[\frac{\mathrm{\α}}{f_m}\sqrt{1-(1-\frac{2n}{N-1})}}{I_0\left(\frac{\mathrm{\α}}{f_m}\right)}$ ；

式中，式中N为频谱长度，M频谱平移次数，α为常数，f_m 为特征频率；

令Kaiser窗矩阵k_(M×N)的傅里叶频谱为K_(M×N)，其中：

；

式中α为常数，ω_n代表为频域变量，N为频谱长度，f_m为特征频率；

将Kaiser矩阵与信号频谱矩阵相乘：

P_(M,N)=H_(M,N)×K_(M,N)；

式中，M为矩阵行数，N为矩阵列数目；

对所有特征频率f_m，对矩阵P_(M,N)各行做傅里叶反变换，得到离散K-S变换结果矩阵KS_(M,N) ：

${\mathrm{KS}}_{M\×N}=\left[\begin{array}{ccccc}\mathrm{ks}(f_1,t_1)& ...& \mathrm{ks}(f_1,t_n)& ...& \mathrm{ks}(f_1,t_N)\\ .& & .& & .\\ .& ...& .& ...& .\\ .& & .& & .\\ \mathrm{ks}(f_m,t_1)& ...& \mathrm{ks}(f_m,t_n)& ...& \mathrm{ks}(f_m,t_N)\\ .& & .& & .\\ .& ...& .& ...& .\\ .& & .& & .\\ \mathrm{ks}(f_M,t_1)& ...& \mathrm{ks}(f_M,t_n)& ...& \mathrm{ks}(f_M,t_N)\end{array}\right]$ ；

其中：

$\mathrm{ks}(m,n)=\frac{2}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}P(m,i)e^{-j\frac{2\mathrm{\πi}}{N}}$ ；

式中，m对应矩阵行元素，n对应矩阵列元素，M为矩阵行大小，N为矩阵列大小，f为对应频率，t对应采样时刻点；

由K-S复数矩阵分离K-S幅值矩阵A_M×N和相位矩阵

，分离公式为：

$A_{M\×N}=\sqrt{\mathrm{Re}{\left({\mathrm{KS}}_{M\×N}\right)}^2+\mathrm{Im}{\left({\mathrm{KS}}_{M\×N}\right)}^2}$ ；

；

式中，A_M×N和的行向量分别表示信号某一采样时刻的幅值和相位随频率变化的分布，它们的列向量分别表示信号某一频率处的幅值和相位随时间变化的分布，N为频谱长度，M频谱平移次数。

与现有技术相比，本发明提出的快速K-S变换，能用于电能质量信号时频参数分析，以FFT→Kaiser窗→IFFT为算法思路，采用Kaiser取代传统S变换的Gauss窗，且能实现随频率变化的Kaiser窗函数宽度、高度的自适应调节，能提高传统S变换时频分析的能量聚集性和频率自适应调节能力；通过特征频率判断，能减少计算量；且为电能质量控制、电力系统保护提供理论指导，为复杂电能质量扰动信号时频参数提取提供重要的理论基础；有利于时频分析方法的发展，为复杂信号处理提供新的研究思路；适用范围广泛，快速K-S变换与电能质量扰动信号时参数提取方法在振动信号分析、故障诊断、模式识别、噪声测量、生物医学检测领域具有广阔的应用前景。

附图说明

图1为本发明一实施例不同β下Kaiser窗时频特性；

（a）Kaiser窗时域特性；(b) Kaiser窗幅频特性；

图2为本发明一实施例基于快速K-S变换的电能质量扰动信号特征提取方法流程图。

图3为本发明一实施例中Kaiser窗与Gauss窗的性能比较图；

(a) Kaiser窗频率自适应调节幅频特性  ； (b) Gauss窗频率自适应调节幅频特性；

(c) f=100Hz的窗函数幅频特性对比 ；(d) f =3000Hz的窗函数幅频特性对比；

图4为本发明一实施例电能质量扰动信号仿真分析正常电压信号快速K-S变换检测结果图；

（a） 电压信号波形；（b）基波幅值曲线；（c）高频幅值和曲线；（d）基波相位曲线；（e）频域特性曲线；

图5为本发明一实施例仅含幅值下降的输入快速K-S变换电压骤降检测结果图；

（a）电压信号波形；（b） 基波幅值曲线；（c） 高频幅值和曲线；（d）基波相位曲线；（e） 频域特性曲线；

图6为本发明一实施例K-S变换和S变换时频能量聚集度对比图；

（a）S变换时域能量分布；(b) K-S变换时域能量分布；

具体实施方式

以下结合附图和实施例对本发明做进一步说明。

实施例1

本实施例包括以下步骤:

在S201步骤中，被测信号采样，将电能质量扰动信号输入电力系统分析仪内进行采样，得被测信号x(n)；

在S202步骤中，快速K-S（凯撒-S）变换：对步骤（1）采样好的被测信号进行K-S变换，得到特征频率f_m对应的扰动信号K-S变换谱，分析出特征频率为f_m对应的扰动信号的幅值、相位和起止、突变时刻；

（3）重复步骤（2），直到计算出所有特征频率f_m对应的扰动信号K-S变换谱，完成扰动信号的K-S变换KS(τ,f)，分析出所有扰动信号的幅值、相位和起止、突变时刻；

通过以上步骤完成电能质量扰动信号时频参数的提取。

所述步骤（2）中，K-S变换的方法包括以下步骤：

在步骤S2011中，傅里叶变换运算，对采样好的被测信号x(n)长度N点进行傅里叶变换运算，得到信号频谱y(n);

在步骤S2012中，特征频率提取，从步骤S2011中计算出的信号频谱y(n)中提取扰动信号特征频率f_m；

在步骤S2013中，自适应Kaiser窗，对步骤S2012提取的扰动信号特征频率f_m ,计算自适应优化Kaiser窗傅里叶变换K(f_m,ω_n)，ω_n代表为频域变量；

在步骤S2014中，频谱平移和频谱相乘，平移信号频谱y(n)到y(n+m)，计算y(m+n)K(f_m,ω_n)， y(n)为信号频谱，m为平移大小，K(f_m,ω_n)为特征频率f_m对应的Kaiser窗傅里叶变换；

在步骤S2015中，傅里叶反变换运算，计算y(m+n)K(f_m,ω_n)傅里叶反变换，得到扰动信号特征频率f_m对应的K-S变换谱ks(f_m,t_n)，t_n对应各采样时刻点；

在步骤S2016中，复数时频矩阵分析，根据K-S变换谱，进行K-S矩阵元素分析，由复数时频矩阵分离得到时频相位矩阵和时频幅值矩阵分析，通过对时频相位矩阵和时频幅值矩阵分析，得到特征频率为f_m的对应的扰动信号的幅值、相位和起止、突变时刻。

所述步骤S2013中，Kaiser窗为一组窗度和高度能调的窗函数，且能自由选择主瓣宽度和旁瓣高度之间的比重，其时域离散表达式为：

$k\left(n\right)=\frac{I_0\left[\mathrm{\β}\sqrt{1-(1-\frac{2n}{N-1})}\right]}{I_0\left(\mathrm{\β}\right)}$ ；

式中，β为形状参数，n为时域变量，N为窗函数长度 ，I₀为第一类修正零阶Bessel函数，其幂级数展开式：

$I_0\left(\mathrm{\β}\right)=1+\underset{m=1}{\overset{\∞}{\mathrm{\Σ}}}{\left[\frac{{(\mathrm{\β}/2)}^m}{m!}\right]}^2$ ；

式中，β为形状参数，用以调节主瓣宽度和旁瓣电平，m取正整数，Kaiser窗频域表达式为：

$K\left(\mathrm{\ω}\right)=\frac{N}{I_0\left(\mathrm{\β}\right)}\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{(\mathrm{N\ω}/2)}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{(\mathrm{N\ω}/2)}^2}}$ ；

式中，sinh为双曲正弦函数，N为窗函数长度，β为形状参数，ω为频域变量；

基于Kaser窗，定义自适应优化Kaiser窗形状参数为:

$\mathrm{\β}=\frac{\mathrm{\α}}{f}$ ；

式中，α为常数，f为被测信号频率，由上式可见，Kaiser窗口的高度和宽度随被测信号频率变化而变化，其幅值与频率大小成正比，宽度随频率的增大而变小，即能实现窗口宽度和高度随信号频率的自动调节。

参照图1，图中给出了β=[0,8,20,35]时Kaiser时域和频域示意图，由图1可见：

当β为0时，Kaiser退化为矩形窗，Kaiser窗旁瓣渐进衰减速度和主瓣宽度随着β的增而增加，时域窗口宽度，旁瓣峰值随β的增加而减小。

当β为0时，Kaiser窗频域主瓣窄，具有较高的频率分辨率；

当β为35时，其旁瓣峰值电平和衰减速度分别为-312.3dB和23dB/oct，时域窗口窄，可实现较高的时域分辨率。

由于Heisenberg测不准原理，信号时宽和带宽是一对矛盾的物理量，即不可能同时得到最优的时间分辨率和频率分辨率。在实际时频分析中，如何根据被测信号的性质，权衡时频分辨率的关系、使两指标达到平衡、满足两者分析要求至关重要。对于电能质量扰动信号，低频段伴有闪变、基波频率波动等扰动，因此在信号频率低时，采用时域窗口宽、频域窗口窄的高频率分辨率窗函数，有利于区别进行闪变、基波频率波动分析等扰动。在信号频率高时，采用时域窗口窄、频域窗口宽的高时间分辨率窗函数，便于对高频扰动进行起止、峰值、突变时刻定位。

Kaiser窗与Gauss窗的性能比较：

针对采样频率f_s=8kHz、采样点数N=1024的离散采样信号，分析对比了相同调节机制下Kaiser窗和Gauss窗的频率特性，图3为Kaiser窗和Gauss窗幅频特性的仿真试验结果。由图3可见：

1）通过Kaiser窗调节因子β、Gauss调节因子v，均可实现窗函数特性随频率的自适应调节；

2）在低频段（如f<500Hz），信号变化相对缓慢，与Gauss窗相比，Kaiser窗旁瓣衰减速率快，旁瓣电平低，更有利于消除各次谐波的混叠和频谱泄漏、闪变扰动，适合对基波和低次谐波频域信号（频率、幅值、相位）的准确提取；

3）随着频率升高，Kaiser和Gauss窗频域窗口宽度逐渐增加，频率变大后（如f>1000Hz）,高次谐波电网信号在很短的时域范围内幅值、频率变化大，此时，与Gauss窗相比，Kaiser窗主瓣较宽，时域窗口较窄，时间分辨率高，更适合对高频扰动信号起止、突变时刻的准确定位。

以上表明：相较于Gauss窗，Kaiser窗的时频适应性更强，自适应调节效果更明显，在各频段的实用性均优于Gauss窗。

所述K-S变换的处理过程为：

将信号x(t)的连续K-S变换KS(τ,f)定义如下：

$\mathrm{KS}(\mathrm{\&tau;},f)={\&Integral;}_{-\&infin;}^{\&infin;}x\left(t\right)k\left(\frac{\mathrm{\&alpha;}(\mathrm{\&tau;}-t)}{f}\right)e^{-i2\mathrm{\&pi;ft}}\mathrm{dt}$ ；

式中，k()为Kaiser窗函数，τ,t为时移因子，f为频率，α为常数，上式为卷积方程，利用卷积定理，将时域卷积运算简化为频域乘积运算；

令长度为N的连续采样离散采样信号x(n)：

x(n)=[x₁…x_n…x_N]；

式中，n对应时域序列，N为采样长度；

采样信号x的傅里叶频谱表示为：

y(n)=[y₁…y_n…y_N]；

式中，n对应频域序列，N为频谱长度；

式中，y_n满足：

$y_n=\frac{2}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_i{e}^{-j\frac{2\mathrm{\&pi;i}}{N}}$ ；

式中，N为采样长度， j代表虚数单位，x_i为采样序列，i为正整数；

特征频率f_m对应的谱线y_n满足以下公式：

y_n-1＜y_n＞y_n+1；

|y_n|＞ε；

式中，ε为设定阈值，n为对应频域序列，谱线y_n对应频率f_m为特征频率；

将信号频谱y(n)平移首尾相连构建频谱矩阵H_M×N：

$H_{M\&CenterDot;N}=\left[\begin{array}{ccccc}{y}_2& y_1& ...& y_N& y_1\\ y_1& y_4& ...& y_1& y_2\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ y_M& y_{M+1}& ...& y_{M-2}& y_{M-1}\\ y_{M+1}& y_{M+2}& ...& y_{M-2}& y_M\end{array}\right]$ ；

式中，N为频谱长度，M频谱平移次数，y对应频谱各谱线；

构建不同频率下Kaiser窗函数，记为Kaiser窗函数矩阵，矩阵各行对应各频率下自适应优化Kaiser窗函数：

$k_{M\&times;N}=\left[\begin{array}{ccccc}k(f_1,t_1)& ...& k(f_1,t_n)& ...& k(f_1,t_N)\\ .& & .& & .\\ .& ...& .& ...& .\\ .& & .& & .\\ k(f_m,t_1)& ...& k(f_m,t_n)& ...& k(f_m,t_N)\\ .& & .& & .\\ .& ...& .& ...& .\\ .& & .& & .\\ k(f_M,t_1)& ...& k(f_M,t_n)& ...& k(f_M,t_N)\end{array}\right]$ ；

式中，N为频谱长度， M频谱平移次数；m对应矩阵行，n对应矩阵列，t对应采样时刻点；

式中：

$k(f_m,t_n)=\frac{I_0[\frac{\mathrm{\&alpha;}}{f_m}\sqrt{1-(1-\frac{2n}{N-1})}}{I_0\left(\frac{\mathrm{\&alpha;}}{f_m}\right)}$ ；

式中，式中N为频谱长度，M频谱平移次数，α为常数，f_m 为特征频率；

令Kaiser窗矩阵k_(M×N)的傅里叶频谱为K_(M×N)，其中：

；

式中α为常数，ω_n代表为频域变量，N为频谱长度，f_m为特征频率；

将Kaiser矩阵与信号频谱矩阵相乘：

P_（M，N）=H_（M，N)×K_(M,N)；

式中，M为矩阵行数，N为矩阵列数目；

对所有特征频率f_m，对矩阵P_(M,N)各行做傅里叶反变换，得到离散K-S变换结果矩阵KS_(M,N) ：

${\mathrm{KS}}_{M\&times;N}=\left[\begin{array}{ccccc}\mathrm{ks}(f_1,t_1)& ...& \mathrm{ks}(f_1,t_n)& ...& \mathrm{ks}(f_1,t_N)\\ .& & .& & .\\ .& ...& .& ...& .\\ .& & .& & .\\ \mathrm{ks}(f_m,t_1)& ...& \mathrm{ks}(f_m,t_n)& ...& \mathrm{ks}(f_m,t_N)\\ .& & .& & .\\ .& ...& .& ...& .\\ .& & .& & .\\ \mathrm{ks}(f_M,t_1)& ...& \mathrm{ks}(f_M,t_n)& ...& \mathrm{ks}(f_M,t_N)\end{array}\right]$ ；

其中：

$\mathrm{ks}(m,n)=\frac{2}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P(m,i)e^{-j\frac{2\mathrm{\&pi;i}}{N}}$ ；

式中，m对应矩阵行元素，n对应矩阵列元素，M为矩阵行大小，N为矩阵列大小，f为对应频率，t对应采样时刻点；

由K-S复数矩阵分离K-S幅值矩阵A_M×N和相位矩阵

，分离公式为：

$A_{M\&times;N}=\sqrt{\mathrm{Re}{\left({\mathrm{KS}}_{M\&times;N}\right)}^2+\mathrm{Im}{\left({\mathrm{KS}}_{M\&times;N}\right)}^2}$ ；

；

式中，A_M×N和

的行向量分别表示信号某一采样时刻的幅值和相位随频率变化的分布，它们的列向量分别表示信号某一频率处的幅值和相位随时间变化的分布，N为频谱长度，M频谱平移次数。

K-S变换和S变换时频能量聚集度对比分析：

相较于Gauss窗函数，Kaiser窗函数频率自适应调节能力强。对谐波扰动信号分析进行S变换和K-S变换的分析，其幅值时频能量分布结果如图6所示。

由图(6)可见，利用K-S变换进行扰动信号时频分析，其幅值矩阵时频能量聚集性较高，谐波轮廓线清晰，利于基波、低次谐波扰动频率、幅值、相位和高次谐波扰动起止、突变时刻的求取。显然，基于Kaiser窗的K-S变换大大改善了S变换时结果曲线模糊、频谱混叠现象的缺点，在不增加计算量的情况下，提高了时频参数提取的准确性。

快速K-S变换运算时间测试：

快速K-S变换运算复杂度与快速S变换相同，按照FFT→Kaiser自适应优化窗→IFFT算法思路对离散信号x(n)进行K-S变换需要进行FFT和IFFT运算，且每个频率点的复数乘加运算次数和FFT运算量相同，运算量大。已知N点的FFT计算量为O(Nlog₂N)，则N点的S变换或K-S变换计算量为O(N²log₂N)，执行一次完整的N点K-S变换需要的时间为

T=2×N×cycles_FFT/f_主频；

式中cycles_FFT为DSP进行一次FFT的时钟周期， f_主频为DSP内核主频。

以德州仪器(TI)的6400系列DSP(主频720MHz)实现为例，执行一次长度分为512、1024的定点16位FFT和IFFT，分别需要13214、26286个时钟周期，由上式可知完成一次K-S变换分别需32.23 ms、74.76ms，计算时间远远超过工频电压周期(20ms)，不利于信号的实时快速分析。

若特征频率个数为r（扰动信号的实际发生个数），则快速K-S变换的计算量可降至O(rNlog₂N),执行一次N点快速K-S变换需要的时间为

T=2×r×cycles_FFT/f_主频；

式中，r为特征频率数，cycles_FFT为DSP进行一次FFT的时钟周期，f_主频为DSP内核主频。

若r=20，对于长度分别为512、1024的定点16位快速K-S变换，计算量仅为K-S变换的3.94%、1.97%，计算时间分别为1.269ms、1.474ms，远小于工频电网基波额定周期(20ms)。表1为S变换、K-S变换和快速K-S变换实际测量的运算时间对比，处理器选用TMS320VC6745。

表1   算法计算时间对比

由表1可知，快速K-S变换时间大大低于S变换、K-S变换运算时间，各采样点数下均能够有效满足电能质量扰动信号时频参数提取的实时性要求。

电能质量扰动信号仿真分析：

利用MATLAB仿真软件，对正常电压信号、电压骤降扰动信号进行快速K-S变换仿真分析。

1. 正常电压信号时频分析：

设供电系统的理想电压A₀为1p.u.、频率f₀为50Hz，正弦输入电压u(t)=A₀sin(2πf₀t)，其快速K-S变换检测结果如图4所示。

由图4可见，正弦电压信号的基波幅值和高频幅值和为常数，即骤降幅值为0，没有发生幅值突变；而基波相位偏差近似为0，即没有发生相位跳变；信号频谱只在基波处存在尖峰，没有谐波分量。上述参数表明信号未发生电压骤降、不含谐波分量，是正常的正弦电压。

2. 电压骤降扰动信号时频分析：

供电系统电压A₀为1p.u.、频率f₀为50Hz的电压骤降扰动信号数学模型为

；

式中，t₁为电压骤降的起始时刻；t₂为电压骤降的结束时刻；电压骤降持续时间t=t₂-t₁，且T<t<6T（T为基波周期）；骤降幅值a=10~90%；相位跳变为φ=0~90°。

图5为仅含幅值下降的输入快速K-S变换电压骤降检测结果，仿真时设定在0.044~0.125s之间发生幅值60%的电压骤降。由图5(b)和图5(c)，信号在0.04375s时发生了瞬时幅值下降，在0.125s时幅值开始恢复，骤降发生期间的幅值为原幅值的60.001%。此外，由图5(d)可知，被检扰动信号的基波电压相位偏差近似为0，即没有发生相位跳变。

Claims (4)

1.基于快速K-S变换电能质量扰动信号时频参数提取方法，包括以下步骤:

（1）被测信号采样：将电能质量扰动信号输入电力系统分析仪内进行采样，得被测信号x(n)；

（2）快速凯撒-S变换：对步骤（1）采样好的被测信号进行K-S变换，得到特征频率f_m对应的扰动信号K-S变换谱，分析出特征频率为f_m对应的扰动信号的幅值、相位和起止、突变时刻；

（3）重复步骤（2），直到计算出所有特征频率f_m对应的扰动信号K-S变换谱，完成扰动信号的K-S变换KS(τ,f)，分析出所有扰动信号的幅值、相位和起止、突变时刻；

通过以上步骤完成电能质量扰动信号时频参数的提取。

2.根据权利要求1所述的基于快速K-S变换电能质量扰动信号时频参数提取方法，其特征在于，所述步骤（2）中，K-S变换的方法包括以下步骤：

（1）傅里叶变换运算：对采样好的被测信号x(n)长度N点进行傅里叶变换运算，得到信号频谱y(n);

（2）特征频率提取:从步骤（1）中计算出的信号频谱y(n)中提取扰动信号特征频率f_m；

（3）自适应Kaiser窗：对步骤（2）提取的扰动信号特征频率f_m ,计算自适应优化Kaiser窗傅里叶变换K(f_m,ω_n)，ω_n代表为频域变量；

（4）频谱平移和频谱相乘：平移信号频谱y(n)到y(n+m)，计算y(m+n)K(f_m,ω_n)， y(n)为信号频谱，m为平移大小，K(f_m,ω_n)为特征频率f_m对应的Kaiser窗傅里叶变换；

（5）傅里叶反变换运算：计算y(m+n)K(f_m,ω_n)傅里叶反变换，得到扰动信号特征频率f_m对应的K-S变换谱ks(f_m,t_n)，t_n对应各采样时刻点；

（6）复数时频矩阵分析：根据K-S变换谱，进行K-S矩阵元素分析，由复数时频矩阵分离得到时频相位矩阵和时频幅值矩阵分析，通过对时频相位矩阵和时频幅值矩阵分析，得到特征频率为f_m的对应的扰动信号的幅值、相位和起止、突变时刻。

3.根据权利要求2所述的基于快速K-S变换电能质量扰动信号时频参数提取方法，其特征在于，所述步骤（3）中，Kaiser窗为一组窗度和高度能调的窗函数，且能自由选择主瓣宽度和旁瓣高度之间的比重，其时域离散表达式为：

$k\left(n\right)=\frac{I_0[\mathrm{\&beta;}\sqrt{1-(1-\frac{2n}{N-1})}}{I_0\left(\mathrm{\&beta;}\right)}$ ；

式中，β为形状参数，n为时域变量，N为窗函数长度 ，I₀为第一类修正零阶Bessel函数，其幂级数展开式：

$I_0\left(\mathrm{\&beta;}\right)=1+\underset{m=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}[\frac{{(\mathrm{\&beta;}/2)}^m}{m!}{]}^2$ ；

式中，β为形状参数，用以调节主瓣宽度和旁瓣电平，m取正整数，Kaiser窗频域表达式为：

$K\left(\mathrm{\&omega;}\right)=\frac{N}{I_0\left(\mathrm{\&beta;}\right)}\frac{\sinh \left(\sqrt{{\mathrm{\&beta;}}^2-{(\mathrm{N\&omega;}/2)}^2}\right)}{\sqrt{{\mathrm{\&beta;}}^2-{(\mathrm{N\&omega;}/2)}^2}}$ ；

式中，sinh为双曲正弦函数，N为窗函数长度，β为形状参数，ω为频域变量；

基于Kaser窗，定义自适应优化Kaiser窗形状参数为:

$\mathrm{\&beta;}=\frac{\mathrm{\&alpha;}}{f}$ ；

式中，α为常数，f为被测信号频率，由上式可见，Kaiser窗口的高度和宽度随被测信号频率变化而变化，其幅值与频率大小成正比，宽度随频率的增大而变小，即能实现窗口宽度和高度随信号频率的自动调节。

4.根据权利要求1或2所述的基于快速K-S变换电能质量扰动信号时频参数提取方法，其特征在于，所述K-S变换的处理过程为：

将信号x(t)的连续K-S变换KS(τ,f)定义如下：

$\mathrm{KS}(\mathrm{\&tau;},f)={\&Integral;}_{-\&infin;}^{\&infin;}x\left(t\right)k\left(\frac{\mathrm{\&alpha;}(\mathrm{\&tau;}-t)}{f}\right)e^{-i2\mathrm{\&pi;ft}}\mathrm{dt}$ ；

式中，k()为Kaiser窗函数，τ,t为时移因子，f为频率，α为常数，上式为卷积方程，利用卷积定理，将时域卷积运算简化为频域乘积运算；

令长度为N的连续采样离散采样信号x(n)：

$x\left(n\right)=[x_1...x_n...x_N]$ ；

式中，n对应时域序列，N为采样长度；

采样信号x的傅里叶频谱表示为：

$y\left(n\right)=[y_1...y_n...y_N]$ ；

式中，n对应频域序列，N为频谱长度；

式中，y_n满足：

$y_n=\frac{2}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_i{e}^{-j\frac{2\mathrm{\&pi;i}}{N}}$ ；

式中，N为采样长度， j代表虚数单位，x_i为采样序列，i为正整数；

特征频率f_m对应的谱线y_n满足以下公式：

$y_{n-1}<y_n<y_{n+1}$ ；

$\left|y_n\right|>\mathrm{\&epsiv;}$ ；

式中，ε为设定阈值，n为对应频域序列，谱线y_n对应频率f_m为特征频率；

将信号频谱y(n)平移首尾相连构建频谱矩阵H_M×N：

$H_{M\&times;N}=\left[\begin{array}{ccccc}{y}_1& y_1& ....& y_N& y_1\\ y_1& y_1& ...& y_1& y_1\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ y_M& y_{M+1}& ...& y_{M-1}& y_{M-1}\\ y_{M+1}& y_{M+1}& ...& y_{M-1}& y_M\end{array}\right]$ ；

式中，N为频谱长度，M频谱平移次数，y对应频谱各谱线；

构建不同频率下Kaiser窗函数，记为Kaiser窗函数矩阵，矩阵各行对应各频率下自适应优化Kaiser窗函数：

$k_{M\&times;N}=\left[\begin{array}{ccccc}k(f_1,t_1)& ...& k(f_1,t_n)& ...& k(f_1,t_N)\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ k(f_m,t_1)& & k(f_m,t_n)& & k(f_m,t_N)\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ k(f_M,t_1)& ...& k(f_M,t_n)& ...& k(f_M,t_N)\end{array}\right]$ ；

式中，N为频谱长度， M频谱平移次数；m对应矩阵行，n对应矩阵列，t对应采样时刻点；

式中：

$k(f_m,t_n)=\frac{I_0[\frac{\mathrm{\&alpha;}}{f_m}\sqrt{1-(1-\frac{2n}{N-1})}}{I_0\left(\frac{\mathrm{\&alpha;}}{f_m}\right)}$ ；

式中，式中N为频谱长度，M频谱平移次数，α为常数，f_m 为特征频率；

令Kaiser窗矩阵k_(M×N)的傅里叶频谱为K_(M×N)，其中：

；

式中α为常数，ω_n代表为频域变量，N为频谱长度，f_m为特征频率；

将Kaiser矩阵与信号频谱矩阵相乘：

$P_{(M,N)}=H_{(M,N)}\&times;K_{(M,N)}$ ；

式中，M为矩阵行数，N为矩阵列数目；

对所有特征频率f_m，对矩阵P_(M,N)各行做傅里叶反变换，得到离散K-S变换结果矩阵KS_(M,N) ：

${\mathrm{KS}}_{M\&times;N}=\left[\begin{array}{ccccc}\mathrm{ks}(f_1,t_1)& ...& \mathrm{ks}(f_1,t_n)& ...& \mathrm{ks}(f_1,t_N)\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ \mathrm{ks}(f_m,t_1)& \mathrm{ks}(f_m,t_n)& & \mathrm{ks}(f_m,t_N)\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ \mathrm{ks}(f_M,t_1)& ...& \mathrm{ks}(f_M,t_n)& ...& \mathrm{ks}(f_M,t_N)\end{array}\right]$ ；

其中：

$\mathrm{ks}(m,n)=\frac{2}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P(m,i)e^{j\frac{2\mathrm{\&pi;i}}{N}}$ ；

式中，m对应矩阵行元素，n对应矩阵列元素，M为矩阵行大小，N为矩阵列大小，f为对应频率，t对应采样时刻点；

由K-S复数矩阵分离K-S幅值矩阵A_M×N和相位矩阵，分离公式为：

$A_{M\&times;N}=\sqrt{\mathrm{Re}{\left({\mathrm{KS}}_{M\&times;N}\right)}^2+\mathrm{Im}{\left({\mathrm{KS}}_{M\&times;N}\right)}^2}$ ；

；

式中，A_M×N和

的行向量分别表示信号某一采样时刻的幅值和相位随频率变化的分布，它们的列向量分别表示信号某一频率处的幅值和相位随时间变化的分布，N为频谱长度，M频谱平移次数。

Images