CN105606893A - 基于空间平滑修正music的电力间谐波检测方法 - Google Patents

Abstract

本发明公开了一种基于空间平滑修正MUSIC的电力间谐波检测方法，对谱分析中的频率分辨率高，无频谱泄露和栅栏效应，谐波间谐波的频率、幅值和相位测量精度高，幅值相对较小的间谐波亦可分辨；本发明首次将空间平滑修正MUSIC方法应用在电力间谐波检测领域，并能对相干电力信号解相干，具有稳健性，适于相干信号源的检测分析；本发明对谐波间谐波信号的幅值和相位检测方面采用最小二乘法，相比混沌检测理论和神经网络模型效率更高，方法简单易于实现，幅值误差小，可在嵌入式系统实现。

Description

基于空间平滑修正MUSIC的电力间谐波检测方法

技术领域

本发明涉及通信领域，尤其是涉及一种基于空间平滑修正MUSIC的电力间谐波检测方法。

背景技术

现代社会由于各种非线性、冲击性和波动性负荷以及高功耗装置的广泛使用，使得电力系统中产生多种形式的谐波和间谐波干扰，严重污染了电网，引起电能质量下降，影响工业生产和人们正常生活。整数倍基波频率的正弦波定义为谐波，非整数倍基波频率的正弦波定义为间谐波，其会增加能量损耗、引发与供电系统的谐振，降低负荷功率因数、还会产生过零点漂移、电压闪变等一系列问题。电能质量遭到严重损坏并影响到电网的安全经济运行。因此必须对其进行治理，而治理的前提是对谐波和间谐波频率和幅值等参数的准确检测。

间谐波的精确检测要比整数次谐波检测困难，能对谐波进行精确检测的方法几乎都不能适应间谐波的检测。由于间谐波幅值一般都比较小，而且时变性和随机性强，要实现对间谐波的精确检测难度很大。已有的间谐波检测方法有快速傅里叶变换法(FFT)、小波变换法(WT)、瞬时无功功率理论等。快速傅里叶变换是当前应用最广泛的一种谐波测量方法，但是栅栏效应和频谱泄漏这个缺点它不能有效克服，这在很大程度上满足不了电力系统高精度的要求。小波变换具有良好的时频局部化特性及自适应性，但小波变换是线性变换，其存在边界效应，频率分辨率也比较低。基于瞬时无功功率理论的谐波检测方法虽具有较好的实时性，但该方法只适用于波形对称、电流不含零序分量且无畸变的三相电压。

已有文献表明，赵娥等人在《中国科技论文在线》发表的“一种基于加权空间平滑的新MUSIC算法”，该文提出利用协方差矩阵的自相关和互相关信息将其子阵阵元数阶子矩阵加权平均，提高了解相干能力，具有较高的频率分辨率和较低的信噪比门限，但无法估计出信号源数即无法检测电力系统谐波间谐波个数，且不能对间谐波的幅值和相位进行测量，因而无法直接应用在电力谐波间谐波检测中。专利申请号为200910154680.1的“基于MUSIC谱估计和HBF神经网络的电力系统间谐波检测方法”适用于间谐波短数据的离线检测和在线检测，此方法频率分辨率高，而且频率和幅值相对准确，但权值的调整学习次数不确定，性能指标有可能长时间达不到误差标准，计算量和复杂度较高，制约实时检测间谐波。

现有技术的不足之处是，间谐波检测的频率分辨率低，增加分析数据长度会使分辨率提高，但不能解决根本问题，而又造成计算量的增加，由于间谐波的频率是基波频率的非整数倍，理想的同步采样难以实现，因此造成测量误差较大；此外，Prony法中的最小二乘估计不能检测电力系统间谐波的个数和频率，无法直接检测电力系统间谐波。

发明内容

发明目的：针对上述现有技术中存在的缺陷，本发明旨在提供一种基于空间平滑修正MUSIC的电力间谐波检测方法。

技术方案：一种基于空间平滑修正MUSIC的电力间谐波检测方法，包括如下步骤：

(1)将模拟的电流电压信号通过采样频率为f_s的模数转换器转化为数字信号，得到的采样数据长度为n，由平滑秩序列法估计出信号源数K；

(2)应用空间平滑修正MUSIC算法对得到的信号源数进行加权空间平滑修正，得到修正的空间平滑矩阵R^fb，根据空间平滑MUSIC谱估计得到信号中各个谐波和间谐波频率；

(3)由空间平滑修正MUSIC谱估计得到信号中各个谐波和间谐波频率后，通过扩展Prony方法中的最小二乘估计计算出这些谐波间谐波分量的幅值和相位。

进一步的，步骤(1)具体包括如下子步骤：

(1.1)设电力系统谐波间谐波波形为：

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

以固定采样频率f_s对电力系统谐波间谐波波形进行采样：

$y\left(n\right)=\underset{i=0}{\overset{H}{\Σ}}{b}_h{z}_h^n+N\left(t\right),n=0,1,...,N-1---(1-2)$

式中，z_h＝exp(j2πf_h△t)，y(n)为采样信号，H为谐波次数，n为采样点数，f_h为第h次谐波的频率，b_h为第h次谐波的幅值；

(1.2)用y(n)生成M×M维矩阵R₀，定义一个(M-k)×M维矩阵I_M-k,j，

I_M-k,j＝[0…0I0…0](1-3)

式中，k是正整数，I为单位矩阵，前j列和后k-j列为0矢量；

将R₀分成交叉重叠矩阵即

$R_0^{\left(k\right)}=\frac{1}{k+1}\underset{i=0}{\overset{k}{\Σ}}{I}_{M-k,i}{R}_0{I}_{M-k,i}^{T}k=1,2,...,M-1---(1-4)$

相关群的数目为：

$Q=\underset{i=1}{\overset{L}{\Σ}}{g}_i---(1-5)$

式中gi(i＝1，2，…，L)，即每个群有i个相干源，L是最大的相关源数；

总的信号源数为 $K=\underset{q=1}{\overset{Q}{\Σ}}{f}_q---\left(1-6\right)$

式中，fq表示第q组相关群的信号源数，即谐波间谐波个数。

进一步的，步骤(2)具体包括如下子步骤：

(2.1)应用空间平滑修正MUSIC算法：窄带情况下的一均匀线阵，第s个阵元接收的数据为

$x_s\left(t\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_s\left({\mathrm{\θ}}_i\right)s_i\left(t\right)+n_s\left(t\right),s=1,2,...,M---(2-1)$

式中，τ_si＝(s-1)dsinθ_i/c，M为阵元数，N为信号源数，c为信号传播速度，d为均匀线阵的间距；

将具有M个阵元的均匀线阵分成p个相互交错的子阵，每个子阵有m个阵元数，M＝p+m-1；

选参考子阵为第一个子阵，则第k个子阵的数据模型为

x_k(t)＝[x_kx_k+1…x_k+m-1]＝AD^(k-1)s(t)+n_k(t)(2-2)

其中

$D=\left[\begin{array}{cccc}{e}^{{\mathrm{j\β}}_1}& 0& \mathrm{...}& 0\\ 0& e^{{\mathrm{j\β}}_2}& \mathrm{...}& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& \mathrm{...}& e^{{\mathrm{j\β}}_N}\end{array}\right]---\left(2-3\right)$

于是该子阵数据协方差矩阵为

R_k＝AD^(k-1)R_s(D^(k-1))^HA^H+σ²I(2-4)

取修正的前向平滑协方差矩阵为：

$R^f=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{R}_i={\mathrm{AR}}_s^f{A}^H+{\mathrm{\σ}}^{2}I---\left(2-5\right)$

其中， $R_s^f=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{D}^{(i-1)}{R}_s{\left(D^{(i-1)}\right)}^H,$

采用后向空间平滑的方法划分子阵，第i个子阵的数据矢量：

$x_i^b\left(t\right)={\left[\begin{array}{cccc}{x}_{M-i+1}& x_{M-i}& ...& x_{M-m-i+2}\end{array}\right]}^{*}---(2-6)$

后向第p-k+1个子阵和前向第k个子阵有如下关系：

$x_{p-k+1}^b\left(t\right)={\mathrm{Jx}}_k^{*}\left(t\right)={\mathrm{JA}}^{*}{D}^{-(k-1)}{s}^{*}\left(t\right)+{\mathrm{Jn}}_k^{*}\left(t\right)---(2-7)$

其中，J为m维的交换矩阵。

因此后向平滑第p-k+1个子阵的数据协方差矩阵为

$R_{p-k+1}^b={\mathrm{JA}}^{*}{D}^{-(k-1)}{R}_s^{*}{\left(D^{-(k-1)}\right)}^H{\left(A^{H}J\right)}^{*}+{\mathrm{\σ}}^{2}I---(2-8)$

且JA^*＝AD^-(m-1)，则：

$R_{p-k+1}^b={\mathrm{AD}}^{-(m+k-2)}{R}_s^{*}{D}^{(m+k-2)}{A}^H+{\mathrm{\σ}}^{2}I---(2-9)$

则修正的后向空间平滑数据矩阵为

$R^b=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{R}_{p-i+1}^b={\mathrm{AR}}_s^b{A}^H+{\mathrm{\σ}}^{2}I---(2-10)$

其中， $R_s^b=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{D}^{(m+i-2)}{R}_s{D}^{-(m+i-2)}$

修正的空间平滑MUSIC算法即：

$R_W^{fb}=\underset{i=1}{\overset{p}{\Σ}}\underset{j=1}{\overset{p}{\Σ}}{R}_{ij}{W}_{ij}^f+\underset{i=1}{\overset{p}{\Σ}}\underset{j=1}{\overset{p}{\Σ}}J{\left(R_{ij}\right)}^{*}{\mathrm{JW}}_{ij}^b=R_W^f+R_W^b---(2-11)$

式中，W^b为后向平滑权矩阵，W^f为前向平滑权矩阵，W^f、W^b均为p×p的加权矩阵，为修正的后向加权矩阵，为修正的前向加权矩阵；

当W^b＝0.5I，W^f＝0.5I，且i＝j时，

$R_W^{fb}=\frac{1}{2}(R^f+R^b)=R^{fb}---(2-12)$

对进行特征值分解：

$R_W^{fb}=U_S{\mathrm{\Σ}}_S{U}_S^H+U_N{\mathrm{\Σ}}_N{U}_N^H---(2-13)$

式中，U_S＝[e₁…e_N]是大特征值对应的信号子空间，U_N＝[e_N+1…e_M]是小特征值对应的噪声子空间，Σ_S为大特征值组成的对角阵,Σ_N为小特征值组成的对角阵；

(2.2)定义多项式：

$f\left(z\right)=e_i^H{\left[\begin{array}{cccc}1& z& ...& z^{M-1}\end{array}\right]}^T,i=N+1,...,M---(2-14)$

当z＝exp(jw)时，p(exp(jw))为空间频率为w的导向矢量，p(exp(jw))＝p_m为信号的导向矢量，信号子空间正交于噪声子空间，信号子空间的导向矢量也正交于噪声子空间，即a^H(θ)U_N＝0，则：

$f\left(z\right)=z^{M-1}{p}^T\left(Z^{-1}\right)U_N{U}_N^{H}p\left(z\right)---(2-15)$

其中，p(z)＝[1z…z^M-1]^T，当z＝exp(jw)时，f(z)的根处在单位圆上，f(z)的阶数为2(M-1)，则f(z)有(M-1)对根，且每对根互相共轭，其中有N个根z₁，…，z_N均分布在单位圆上，且

z_i＝exp(jw_i)i＝1≤i≤N(2-16)

则频率估计为

$f_i=\frac{f_{s}angle\left(z_i\right)}{2\mathrm{\π}}---(2-17).$

进一步的，步骤(3)具体包括：

空间平滑MUSIC谱估计得到信号中各个谐波和间谐波频率fi以及个数K后，利用Prony方法求取各频率分量的幅值和相位：

$Zb=\hat{y}---(3-1)$

式中，

$Z=\left[\begin{array}{cccc}1& 1& \mathrm{...}& 1\\ z_1& z_2& \mathrm{...}& z_K\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ z_1^{N-1}& z_2^{N-1}& \mathrm{...}& z_K^{N-1}\end{array}\right]---(3-2)$

b＝[b₁,b₂,…,b_p]^T

(3-3)

$\hat{y}={\left[\begin{array}{cccc}\hat{y}\left(0\right)& \hat{y}\left(1\right)& ...& \hat{y}(N-1)\end{array}\right]}^T---(3-4)$

由于z_i各不相同，矩阵Z是满列秩，方程(3-1)的最小二乘解为

$b={\[Z^{H}Z\]}^{-1}{Z}^H\hat{y}---(3-5)$

求出b_i后根据式(3-5)来计算信号中各频率分量的幅值A_i和相位

进一步的，步骤(1.2)中所述R₀即：

$R_0=\left[\begin{array}{cccc}y\left(0\right)& y\left(1\right)& \mathrm{...}& y(M-1)\\ y\left(1\right)& y\left(2\right)& \mathrm{...}& y\left(M\right)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ y(M-1)& y\left(M\right)& \mathrm{...}& y(2M-2)\end{array}\right]---(1-7)$

式中，M为阵元数，2M为总采样点数。

有益效果：

(1)本发明对谱分析中的频率分辨率高，无频谱泄露和栅栏效应，谐波间谐波的频率、幅值和相位测量精度高，幅值相对较小的间谐波亦可分辨；

(2)本发明首次将空间平滑修正MUSIC方法应用在电力间谐波检测领域，并能对相干电力信号解相干，具有稳健性，适于相干信号源的检测分析；

(3)本发明对谐波间谐波信号的幅值和相位检测方面采用最小二乘法，相比混沌检测理论和神经网络模型效率更高，方法简单易于实现，幅值误差小，可在嵌入式系统实现。

附图说明

图1是空间平滑修正MUSIC在电力间谐波检测中的应用的框图；

图2是前向空间平滑算法原理图；

图3是后向空间平滑算法原理图。

具体实施方式

下面结合附图对本发明作进一步的说明，如图1所示，一种基于空间平滑修正MUSIC的电力间谐波检测方法，包括如下步骤：

(1)将模拟的电流电压信号通过采样频率为f_s的模数转换器转化为数字信号，得到的采样数据长度为n，由平滑秩序列法估计出信号源数K；

步骤(1)具体包括如下子步骤：

(1.1)设电力系统谐波间谐波波形为：

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

以固定采样频率f_s对电力系统谐波间谐波波形进行采样：

$y\left(n\right)=\underset{i=0}{\overset{H}{\Σ}}{b}_h{z}_h^n+N\left(t\right),n=0,1,...,N-1---(1-2)$

式中，z_h＝exp(j2πf_h△t)，y(n)为采样信号，H为谐波次数，n为采样点数，f_h为第h次谐波的频率，b_h为第h次谐波的幅值；

(1.2)用y(n)生成M×M维矩阵R₀，

$R_0=\left[\begin{array}{cccc}y\left(0\right)& y\left(1\right)& \mathrm{...}& y(M-1)\\ y\left(1\right)& y\left(2\right)& \mathrm{...}& y\left(M\right)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ y(M-1)& y\left(M\right)& \mathrm{...}& y(2M-2)\end{array}\right]---(1-7)$

式中，M为阵元数，2M为总采样点数。

定义一个(M-k)×M维矩阵I_M-k,j，

I_M-k,j＝[0…0I0…0](1-3)

式中，k是正整数，I为单位矩阵，前j列和后k-j列为0矢量；

将R₀分成交叉重叠矩阵

$R_0=\left[\begin{array}{cccc}y\left(0\right)& y\left(1\right)& \mathrm{...}& y(M-1)\\ y\left(1\right)& y\left(2\right)& \mathrm{...}& y\left(M\right)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ y(M-1)& y\left(M\right)& \mathrm{...}& y(2M-2)\end{array}\right]---(1-7)$

式中，M为阵元数，2M为总采样点数。

$R_0^{\left(k\right)}=\frac{1}{k+1}\underset{i=0}{\overset{k}{\Σ}}{I}_{M-k,i}{R}_0{I}_{M-k,i}^{T}k=1,2,...,M-1---(1-4)$

上式中的交叉重叠矩阵序列属于前向空间平滑矩阵。假设信号源有L组相关源的群组成，分别表示为g_i(i＝1,2,…,L)，如i＝1则表明该群是单个独立信号；i＝3说明该群有三个相干源，L是最大的相关源数。若g₃＝4则说明有四个相关群，每群有三个相干源，则相关群的数目为：

$Q=\underset{i=1}{\overset{L}{\Σ}}{g}_i---(1-5)$

式中gi(i＝1，2，…，L)，即每个群有i个相干源，L是最大的相关源数；

总的信号源数为 $K=\underset{q=1}{\overset{Q}{\Σ}}{f}_q---\left(1-6\right)$

式中，fq表示第q组相关群的信号源数，即谐波间谐波个数。

(2)应用空间平滑修正MUSIC算法对得到的信号源数进行加权空间平滑修正，得到修正的空间平滑矩阵R^fb，根据空间平滑MUSIC谱估计得到信号中各个谐波和间谐波频率；

MUSIC谱估计有求根和频谱搜索两种方式，下面以求根MUSIC为例介绍空间平滑MUSIC谱估计方法。当信号源相干时，阵元接收的数据协方差矩阵的秩降为1，信号源数大于信号子空间的维数，这会使得有些相干源的噪声子空间和导向矢量不完全正交，从而无法正确估计信号源信息。空间平滑MUSIC通过一系列处理使得信号协方差矩阵的秩得到有效恢复，对修正后的空间平滑矩阵R^fb进行特征值分解，得到噪声和信号子空间，由这两个子空间的正交性构造出多项式，再通过求解此多项式的根得到谐波间谐波的检测频率值。

步骤(2)具体包括如下子步骤：

(2.1)应用空间平滑修正MUSIC算法：窄带情况下的一均匀线阵，第s个阵元接收的数据为

$x_s\left(t\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_s\left({\mathrm{\θ}}_i\right)s_i\left(t\right)+n_s\left(t\right),s=1,2,...,M---(2-1)$

式中，τ_si＝(s-1)dsinθ_i/c，M为阵元数，N为信号源数，c为信号传播速度，d为均匀线阵的间距；

前、后向空间平滑MUSIC算法恢复满秩协方差矩阵需要N个子阵，空间平滑修正算法只需要0.5N个子阵恢复，因此修正平滑阵列孔径损失小。当有N个信源、每个子阵的阵元数为N+1时，要分辨这些信源，修正平滑共需M＝1.5N个阵元，而单向平滑至少需要M＝2N个阵元。

如图2所示，将具有M个阵元的均匀线阵分成p个相互交错的子阵，每个子阵有m个阵元数，M＝p+m-1；

选参考子阵为第一个子阵，则第k个子阵的数据模型为

x_k(t)＝[x_kx_k+1…x_k+m-1]＝AD^(k-1)s(t)+n_k(t)(2-2)

其中

$D=\left[\begin{array}{cccc}{e}^{{\mathrm{j\β}}_1}& 0& \mathrm{...}& 0\\ 0& e^{{\mathrm{j\β}}_2}& \mathrm{...}& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& \mathrm{...}& e^{{\mathrm{j\β}}_N}\end{array}\right]---\left(2-3\right)$

于是该子阵数据协方差矩阵为

R_k＝AD^(k-1)R_s(D^(k-1))^HA^H+σ²I(2-4)

取修正的前向平滑协方差矩阵为：

$R^f=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{R}_i={\mathrm{AR}}_s^f{A}^H+{\mathrm{\σ}}^{2}I---\left(2-5\right)$

其中， $R_s^f=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{D}^{(i-1)}{R}_s{\left(D^{(i-1)}\right)}^H,$

如图3所示，采用后向空间平滑的方法划分子阵，第i个子阵的数据矢量：

$x_i^b\left(t\right)={\left[\begin{array}{cccc}{x}_{M-i+1}& x_{M-i}& ...& x_{M-m-i+2}\end{array}\right]}^{*}---(2-6)$

后向第p-k+1个子阵和前向第k个子阵有如下关系：

$x_{p-k+1}^b\left(t\right)={\mathrm{Jx}}_k^{*}\left(t\right)={\mathrm{JA}}^{*}{D}^{-(k-1)}{s}^{*}\left(t\right)+{\mathrm{Jn}}_k^{*}\left(t\right)---(2-7)$

其中，J为m维的交换矩阵。

因此后向平滑第p-k+1个子阵的数据协方差矩阵为

$R_{p-k+1}^b={\mathrm{JA}}^{*}{D}^{-(k-1)}{R}_s^{*}{\left(D^{-(k-1)}\right)}^H{\left(A^{H}J\right)}^{*}+{\mathrm{\σ}}^{2}I---(2-8)$

且JA^*＝AD^-(m-1)，则：

$R_{p-k+1}^b={\mathrm{AD}}^{-(m+k-2)}{R}_s^{*}{D}^{(m+k-2)}{A}^H+{\mathrm{\σ}}^{2}I---(2-9)$

则修正的后向空间平滑数据矩阵为

$R^b=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{R}_{p-i+1}^b={\mathrm{AR}}_s^b{A}^H+{\mathrm{\σ}}^{2}I---(2-10)$

其中， $R_s^b=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{D}^{(m+i-2)}{R}_s{D}^{-(m+i-2)}$

修正的空间平滑MUSIC算法即：

$R_W^{fb}=\underset{i=1}{\overset{p}{\Σ}}\underset{j=1}{\overset{p}{\Σ}}{R}_{ij}{W}_{ij}^f+\underset{i=1}{\overset{p}{\Σ}}\underset{j=1}{\overset{p}{\Σ}}J{\left(R_{ij}\right)}^{*}{\mathrm{JW}}_{ij}^b=R_W^f+R_W^b---(2-11)$

式中，W^b为后向平滑权矩阵，W^f为前向平滑权矩阵，W^f、W^b均为p×p的加权矩阵，为修正的后向加权矩阵，为修正的前向加权矩阵；上式的目的是实现对相干信号源的解相干，其实质是对数据协方差矩阵的各子阵进行加权求和。

当W^b＝0.5I，W^f＝0.5I，且i＝j时，

$R_W^{fb}=\frac{1}{2}(R^f+R^b)=R^{fb}---(2-12)$

对进行特征值分解：

$R_W^{fb}=U_S{\mathrm{\Σ}}_S{U}_S^H+U_N{\mathrm{\Σ}}_N{U}_N^H---(2-13)$

式中，U_S＝[e₁…e_N]是大特征值对应的信号子空间，U_N＝[e_N+1…e_M]是小特征值对应的噪声子空间，Σ_S为大特征值组成的对角阵,Σ_N为小特征值组成的对角阵；

(2.2)定义多项式：

$f\left(z\right)=e_i^H{\left[\begin{array}{cccc}1& z& ...& z^{M-1}\end{array}\right]}^T,i=N+1,...,M---(2-14)$

当z＝exp(jw)时，p(exp(jw))为空间频率为w的导向矢量，p(exp(jw))＝p_m为信号的导向矢量，信号子空间正交于噪声子空间，信号子空间的导向矢量也正交于噪声子空间，即a^H(θ)U_N＝0，则：

$f\left(z\right)=z^{M-1}{p}^T\left(Z^{-1}\right)U_N{U}_N^{H}p\left(z\right)---(2-15)$

其中，p(z)＝[1z…z^M-1]^T，当z＝exp(jw)时，f(z)的根处在单位圆上，f(z)的阶数为2(M-1)，则f(z)有(M-1)对根，且每对根互相共轭，其中有N个根z₁，…，z_N均分布在单位圆上，且

z_i＝exp(jw_i)i＝1≤i≤N(2-16)

在实际情况下，噪声存在于信号中，数据矩阵存在误差，这时只需求取接近于单位圆上的N个根即可，则频率估计为

$f_i=\frac{f_{s}angle\left(z_i\right)}{2\mathrm{\π}}---(2-17).$

(3)由空间平滑修正MUSIC谱估计得到信号中各个谐波和间谐波频率后，通过扩展Prony方法中的最小二乘估计计算出这些谐波间谐波分量的幅值和相位。

步骤(3)具体包括：

式(1-2)所示的信号可表示为

$\hat{y}\left(n\right)=\underset{i=0}{\overset{H}{\Σ}}{b}_i{z}_i^n---(2-18)$

式中，z_i＝exp(j2πf_iΔt)，用代替y(n)。

由以上步骤得到信号中各个谐波和间谐波频率fi以及个数K后，z_i就成为已知量，把指数模型(31)简化为参数b_i的线性方程，用矩阵形式表示，即：

$Zb=\hat{y}---(3-1)$

式中，

$Z=\left[\begin{array}{cccc}1& 1& \mathrm{...}& 1\\ z_1& z_2& \mathrm{...}& z_K\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ z_1^{N-1}& z_2^{N-1}& \mathrm{...}& z_K^{N-1}\end{array}\right]---(3-2)$

b＝[b₁,b₂,…,b_p]^T

(3-3)

$\hat{y}={\left[\begin{array}{cccc}\hat{y}\left(0\right)& \hat{y}\left(1\right)& ...& \hat{y}(N-1)\end{array}\right]}^T---(3-4)$

由于z_i各不相同，矩阵Z是满列秩，方程(3-1)的最小二乘解为

$b={\[Z^{H}Z\]}^{-1}{Z}^H\hat{y}---(3-5)$

求出b_i后根据式(3-5)来计算信号中各频率分量的幅值A_i和相位

Claims (5)

1.一种基于空间平滑修正MUSIC的电力间谐波检测方法，其特征在于，包括如下步骤：

(1)将模拟的电流电压信号通过采样频率为f_s的模数转换器转化为数字信号，得到的采样数据长度为n，由平滑秩序列法估计出信号源数K；

(2)应用空间平滑修正MUSIC算法对得到的信号源数进行加权空间平滑修正，得到修正的空间平滑矩阵R^fb，根据空间平滑MUSIC谱估计得到信号中各个谐波和间谐波频率；

(3)由空间平滑修正MUSIC谱估计得到信号中各个谐波和间谐波频率后，通过扩展Prony方法中的最小二乘估计计算出这些谐波间谐波分量的幅值和相位。

2.根据权利要求1所述的一种基于空间平滑修正MUSIC的电力间谐波检测方法，其特征在于，所述步骤(1)具体包括如下子步骤：

(1.1)设电力系统谐波间谐波波形为：

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

以固定采样频率f_s对电力系统谐波间谐波波形进行采样：

$y\left(n\right)=\underset{i=0}{\overset{H}{\Σ}}{b}_h{z}_h^n+N\left(t\right),n=0,1,...,N-1---(1-2)$

式中，z_h＝exp(j2πf_h△t)，y(n)为采样信号，H为谐波次数，n为采样点数，f_h为第h次谐波的频率，b_h为第h次谐波的幅值；

(1.2)用y(n)生成M×M维矩阵R₀，定义一个(M-k)×M维矩阵I_M-k,j，

I_M-k,j＝[0…0I0…0](1-3)

式中，k是正整数，I为单位矩阵，前j列和后k-j列为0矢量；

将R₀分成交叉重叠矩阵即

$R_0^{\left(k\right)}=\frac{1}{k+1}\underset{i=0}{\overset{k}{\Σ}}{I}_{M-k,i}{R}_0{I}_{M-k,i}^T,k=1,2,...,M-1---(1-4)$

相关群的数目为：

$Q=\underset{i=1}{\overset{L}{\Σ}}{g}_i---(1-5)$

式中gi(i＝1，2，…，L)，即每个群有i个相干源，L是最大的相关源数；

总的信号源数为 $K=\underset{q=1}{\overset{Q}{\Σ}}{f}_q---(1-6)$

式中，fq表示第q组相关群的信号源数，即谐波间谐波个数。

3.根据权利要求1所述的一种基于空间平滑修正MUSIC的电力间谐波检测方法，其特征在于，所述步骤(2)具体包括如下子步骤：

(2.1)应用空间平滑修正MUSIC算法：窄带情况下的一均匀线阵，第s个阵元接收的数据为

$x_s\left(t\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_s\left({\mathrm{\θ}}_i\right)s_i\left(t\right)+n_s\left(t\right),s=1,2,...,M---(2-1)$

式中，τ_si＝(s-1)dsinθ_i/c，M为阵元数，N为信号源数，c为信号传播速度，d为均匀线阵的间距；

将具有M个阵元的均匀线阵分成p个相互交错的子阵，每个子阵有m个阵元数，M＝p+m-1；

选参考子阵为第一个子阵，则第k个子阵的数据模型为

x_k(t)＝[x_kx_k+1…x_k+m-1]＝AD^(k-1)s(t)+n_k(t)(2-2)

其中

$D=\left[\begin{array}{cccc}{e}^{{\mathrm{j\β}}_1}& 0& ...& 0\\ 0& e^{{\mathrm{j\β}}_2}& ...& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& ...& e^{{\mathrm{j\β}}_N}\end{array}\right]---(2-3)$

于是该子阵数据协方差矩阵为

R_k＝AD^(k-1)R_s(D^(k-1))^HA^H+σ²I(2-4)

取修正的前向平滑协方差矩阵为：

$R^f=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{R}_i={\mathrm{AR}}_s^f{A}^H+{\mathrm{\σ}}^{2}I---(2-5)$

其中， $R_s^f=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{D}^{(i-1)}{R}_s{\left(D^{(i-1)}\right)}^H,$

采用后向空间平滑的方法划分子阵，第i个子阵的数据矢量：

$x_i^b\left(t\right)={\left[\begin{array}{cccc}{x}_{M-i+1}& x_{M-i}& ...& x_{M-m-i+2}\end{array}\right]}^{*}---(2-6)$

后向第p-k+1个子阵和前向第k个子阵有如下关系：

$x_{p-k+1}^b\left(t\right)={\mathrm{Jx}}_k^{*}\left(t\right)={\mathrm{JA}}^{*}{D}^{-(k-1)}{s}^{*}\left(t\right)+{\mathrm{Jn}}_k^{*}\left(t\right)---(2-7)$

其中，J为m维的交换矩阵。

因此后向平滑第p-k+1个子阵的数据协方差矩阵为

$R_{p-k+1}^b={\mathrm{JA}}^{*}{D}^{-(k-1)}{R}_s^{*}{\left(D^{-(k-1)}\right)}^H{\left(A^{H}J\right)}^{*}+{\mathrm{\σ}}^{2}I---(2-8)$

且JA^*＝AD^-(m-1)，则：

$R_{p-k+1}^b={\mathrm{AD}}^{-(m+k-2)}{R}_s^{*}{D}^{(m+k-2)}{A}^H+{\mathrm{\σ}}^{2}I---(2-9)$

则修正的后向空间平滑数据矩阵为

$R^b=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{R}_{p-i+1}^b={\mathrm{AR}}_s^b{A}^H+{\mathrm{\σ}}^{2}I---(2-10)$

其中， $R_s^b=\frac{1}{p}\underset{i=1}{\overset{p}{\Σ}}{D}^{(m+i-2)}{R}_s{D}^{-(m+i-2)}$

修正的空间平滑MUSIC算法即：

$R_W^{fb}=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{R}_{ij}{W}_{ij}^f+\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\underset{j=1}{\overset{p}{\mathrm{\Σ}}}J{\left(R_{ij}\right)}^{*}{\mathrm{JW}}_{ij}^b=R_W^f+R_W^b---(2-11)$

式中，W^b为后向平滑权矩阵，W^f为前向平滑权矩阵，W^f、W^b均为p×p的加权矩阵，为修正的后向加权矩阵，为修正的前向加权矩阵；

当W^b＝0.5I，W^f＝0.5I，且i＝j时，

$R_W^{fb}=\frac{1}{2}(R^f+R^b)=R^{fb}---(2-12)$

对进行特征值分解：

$R_W^{fb}=U_S{\mathrm{\Σ}}_S{U}_S^H+U_N{\mathrm{\Σ}}_N{U}_N^H---(2-13)$

式中，U_S＝[e₁…e_N]是大特征值对应的信号子空间，U_N＝[e_N+1…e_M]是小特征值对应的噪声子空间，∑_S为大特征值组成的对角阵,∑_N为小特征值组成的对角阵；

(2.2)定义多项式：

$f\left(z\right)=e_i^H{\left[\begin{array}{cccc}1& z& ...& z^{M-1}\end{array}\right]}^T,i=N+1,...,M---(2-14)$

当z＝exp(jw)时，p(exp(jw))为空间频率为w的导向矢量，p(exp(jw))＝p_m为信号的导向矢量，信号子空间正交于噪声子空间，信号子空间的导向矢量也正交于噪声子空间，即a^H(θ)U_N＝0，则：

$f\left(z\right)=z^{M-1}{p}^T\left(Z^{-1}\right)U_N{U}_N^{H}p\left(z\right)---(2-15)$

其中，p(z)＝[1z…z^M-1]^T，当z＝exp(jw)时，f(z)的根处在单位圆上，f(z)的阶数为2(M-1)，则f(z)有(M-1)对根，且每对根互相共轭，其中有N个根z₁，…，z_N均分布在单位圆上，且

z_i＝exp(jw_i)i＝1≤i≤N(2-16)

则频率估计为

$f_i=\frac{f_{s}angle\left(z_i\right)}{2\mathrm{\π}}---(2-17).$

4.根据权利要求1所述的一种基于空间平滑修正MUSIC的电力间谐波检测方法，其特征在于，所述步骤(3)具体包括：

空间平滑MUSIC谱估计得到信号中各个谐波和间谐波频率fi以及个数K后，利用Prony方法求取各频率分量的幅值和相位：

$Zb=\hat{y}---(3-1)$

式中，

b＝[b₁,b₂,…,b_p]^T(3-3)

$\hat{y}={\left[\begin{array}{cccc}\hat{y}\left(0\right)& \hat{y}\left(1\right)& ...& \hat{y}(N-1)\end{array}\right]}^T---(3-4)$

由于z_i各不相同，矩阵Z是满列秩，方程(3-1)的最小二乘解为

$b={\[Z^{H}Z\]}^{-1}{Z}^H\hat{y}---(3-5)$

求出b_i后根据式(3-5)来计算信号中各频率分量的幅值A_i和相位

5.根据权利要求1所述的一种基于空间平滑修正MUSIC的电力间谐波检测方法，其特征在于，步骤(1.2)中所述R₀即：

$R_0=\left[\begin{array}{cccc}y\left(0\right)& y\left(1\right)& ...& y(M-1)\\ y\left(1\right)& y\left(2\right)& ...& y\left(M\right)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ y(M-1)& y\left(M\right)& ...& y(2M-2)\end{array}\right]---(1-7)$

式中，M为阵元数，2M为总采样点数。