【HHT孤岛检测】基于希尔伯特-黄变换(Hilbert-Huang Transform，HHT)GRNN的孤岛检测算法

1.软件版本

MATLAB2013b

2.本算法理论知识

        该算法在总结现有的电网孤岛检测方法基础之上，将时频分析方法应用到了孤岛检测之中，给出了基于希尔伯特-黄变换的暂态分析新能源集中式孤岛检测方法。该方法对于电网中的公共耦合点(PCC)处电压进行经验模态分解(EMD)，获取固有模态函数(IMF)，实现暂态特征分量的提取和分析，并通过GRNN神经网络实现孤岛现象的检测。

        希尔伯特-黄变换(Hilbert-Huang  Transform，HHT)是一种可以精确识别的自适应信号处理方法，其主要用于非平稳信号的处理与分析。希尔伯特-黄变换的核心技术是EMD经验模态分解和HSA希尔伯特谱分析。

       结合GRNN广义回归神经网络，实现孤岛检测，整个算法流程图如下：

          最后将该方法通过对公共耦合点(PCC)处电压进行改进型的经验模态分解(EMD)，并得到固有模态函数(IMF)，实现暂态特征分量的提取和分析，在检测特征量的同时，利用电压相位变化量、功率波动量等多种判据综合判断孤岛现象的发生。整个算法的流程图如下图所示：

        从图的算法流程可知，基于希尔伯特-黄变换和GRNN神经网络的新型孤岛检测算法，该算法首先对公共耦合点的电压信号进行采集，然后将电压通过希尔伯特-黄变换提取器谐波特征数据，将特征信息通过GRNN广义回归神经网络进行判决是否产生孤岛，因此，该方法具有不依然功率、阈值等参数的优势，通过神经网络进行孤岛检测的方法可以避免直接通过观测判断而导致的孤岛误判等问题，并且可以在最大程度上消除孤岛检测忙去的影响。

3.部分核心代码

function YH=func_hilbert(imf);
N  = length(imf);  
k  = 0:N-1;  
% t  = Time1;  
% signal transform  
% 结论：sin信号Hilbert变换后为cos信号  
y  = imf;  
yh = hilbert(y);    % matlab函数得到信号是合成的复信号  
yi = imag(yh);      % 虚部为书上定义的Hilbert变换  
 
% 检验两次Hilbert变换的结果（理论上为原信号的负值）  
% 结论：两次Hilbert变换的结果为原信号的负值  
yih = hilbert(yi);  
yii = imag(yih);  
% 信号与其Hilbert变换的正交性  
% 结论：Hilbert变换后的信号与原信号正交  
% 谱分析  
% 结论：Hilbert变换后合成的复信号的谱没有大于奈氏频率的频谱，即其谱为单边的  
NFFT = 2^nextpow2(N);  
% f = fs*linspace(0,1,NFFT);  
Y  = fft(y, NFFT)/N;  
YH = fft(yh, NFFT)/N;  
  
 

function [imf,res] = func_emd(x,t);

add     = t(2)-t(1);
defstop = [0.05,0.5,0.05];
stop    = defstop;
 
S       = size(x);
if S(1) > 1
   x = x';
end
S       = size(t);
if S(1) > 1
   t = t';
end
S       = size(stop);
if S(1) > 1
   stop = stop';
   S    = size(stop);
end

if S(2) < 3
   stop(3)=defstop(3);
end

if S(2) < 2
   stop(2)=defstop(2);
end

sd  = stop(1);
sd2 = stop(2);
tol = stop(3);

MAXITERATIONS = 200;
NBSYM         = 2;
lx            = length(x);
sdt(lx)       = 0;
sdt           = sdt+sd;
sd2t(lx)      = 0;
sd2t          = sd2t+sd2;
ner           = lx;
nzr           = lx;
r             = x;
imf           = [];
k             = 1;
nbit          = 0;
NbIt          = 0;

while ner > 2  
  m    = r;
  mp   = m;
  sx   = sd+1;
  test = 0;
  
  [indmin,indmax,indzer] = func_minmaxzo(m);
  lm=length(indmin);
  lM=length(indmax);
  nem=lm + lM;
  nzm=indzer;
  
  j=1;
  
  % sifting loop
  while ( mean(sx > sd) > tol | any(sx > sd2) | (abs(nzm-nem)>1)) & (test == 0) & nbit<MAXITERATIONS

        [maxbl,minbl]=func_baoluo( m,t,indmax,indmin );
        midbl = (maxbl + minbl)/2;

        m = m - midbl;

        % eval（abs(midbl))./(abs(maxbl-minbl));
        s = mean(sx);

        if nem < 3
           test = 1;
        end

        mp   = m;
        nbit = nbit+1;
        NbIt = NbIt+1;
  end
  
  imf(k,:) = m;
  nbits(k) = nbit;
  k        = k+1;
  r        = r - m;
  [indmin,indmax] = func_minmaxzo(r);
  ner = length(indmin) + length(indmax);
  
  nbit=1;

  if (max(r) - min(r)) < (1e-10)*(max(x) - min(x))
    break
  end
  
end

imf(k,:) = r;

clc;
clear;
close all;
warning off;

 

load data.mat 

P_train = [Train(:,2)/1000,Train(:,3)/100,Train(:,4)];
T_trian = Train(:,1);

P_test = [Test(:,2)/1000,Test(:,3)/100,Test(:,4)];
T_test = Test(:,1);

net = newgrnn(P_train',T_trian',0.0605);
tic;
y   = round([net(P_test')]');

length(find(y==T_test))
toc;

4.操作步骤与仿真结论

EMD分解之后的IMF仿真

 IMF经过Hilbert变换之后得到的瞬时频率仿真结果

       根据采集得到的电压信号的幅度，基波频率以及希尔伯特-黄变换之后的THD值三个指标作为特征数据。

        从上述表格的测试结论可知，算法获得了98%检测率，这说明对100组未知样本测试，可以成功检测98组，且这100组样本的检测时间仅需0.038548s，

5.参考文献

[1]Moziza C J. Interconnection protection of IPP generators at commercial/industrial facilities[J]. IEEE Transactions on Industry Applications, 2001, 37(03):681-689．

[2]Katiraei F, Iravani  M  R,  Lehn  P  W.  Micro-grid  autonomous  operation  during  and  subsequent  to islanding process[J]. IEEE Transactions on Power Delivery, 2005, 20(01):248-257.

[3]IEEE Std.1547-2003, IEEE Standard for Interconnecting Distributed Resources With Electric Power System.A02-51