多重分形去趨勢波動(dòng)分析及改進(jìn)決策樹在電能質(zhì)量分析中的應(yīng)用

張淑清，張 赟，劉海濤，胡 皓，李 華， 姚玉永，劉 勇，王 濤

(1.燕山大學(xué) 電氣工程學(xué)院，河北 秦皇島 066004；2.東北大學(xué) 信息科學(xué)與工程學(xué)院，遼寧 沈陽 110004；3.國網(wǎng)冀北電力有限公司唐山供電公司，河北 唐山 063000)

1 引 言

經(jīng)濟(jì)社會(huì)的不斷發(fā)展，各種新型電氣設(shè)備的使用，產(chǎn)生了大量的非線性負(fù)荷和沖擊性負(fù)荷。這些非線性負(fù)載的使用不可避免地造成電網(wǎng)的幅值和頻率發(fā)生偏差，造成電能質(zhì)量(power quality)的下降，不僅影響日常生活還會(huì)造成巨大的經(jīng)濟(jì)損失甚至電力事故。電能質(zhì)量信號(hào)具有非線性和非平穩(wěn)性的特點(diǎn)。目前電能質(zhì)量擾動(dòng)信號(hào)的分類識(shí)別方法主要為小波變換、S變換、EMD分解、HHT變換等時(shí)頻域分析法[1～4]。

多重分形去趨勢波動(dòng)分析(multifractal detrended fluctuation analysis,MFDFA)是在DFA基礎(chǔ)上提出的一種非線性時(shí)間序列的分析方法，可以有效消除干擾趨勢，估計(jì)多重分形譜[5]。文獻(xiàn)[6]采用MFDFA方法對(duì)風(fēng)電廠時(shí)間序列進(jìn)行分析，實(shí)驗(yàn)表明風(fēng)速變化與多重分形參數(shù)有一定的關(guān)聯(lián)性，有效判斷風(fēng)速的變化趨勢。文獻(xiàn)[7]針對(duì)非線性的脈沖信號(hào)，提出多重分形消除趨勢波動(dòng)分析和PSO-SVM的方法，有效識(shí)別了4種脈沖信號(hào)。文獻(xiàn)[8]提出了基于多重分形降趨波動(dòng)分析法的諧振接地系統(tǒng)故障選線新判據(jù)，仿真結(jié)果表明該判據(jù)抗噪能力強(qiáng)，計(jì)算速度快,克服了消弧線圈和故障條件等因素對(duì)故障選線的影響。文獻(xiàn)[9]針對(duì)齒輪箱故障診斷提出了一種基于多重分形的趨勢波動(dòng)分析和改進(jìn)的K均值集群的方法。文獻(xiàn)[10]針對(duì)EEG信號(hào)使用MFDFA方法提取了4個(gè)特征集，用作支持向量機(jī)和k-最近鄰分類器的輸入，實(shí)驗(yàn)證明所提出的MFDFA輔助特征提取方法獲得較高的分類精度。目前，MFDFA方法廣泛應(yīng)用于股票交易數(shù)據(jù)[11]、機(jī)械振動(dòng)信號(hào)[12]、電力系統(tǒng)故障選線[13]、人體腦電信號(hào)[14]等。

決策樹是由一種自上而下的遞歸算法，廣泛應(yīng)用于數(shù)據(jù)的預(yù)測，挖掘與分類。其核心算法最早是由Quinlan提出的ID3算法。C4.5是在ID3算法的基礎(chǔ)上進(jìn)行的改進(jìn)，以信息增益率代替ID3算法中的信息增益作為屬性選擇判斷條件[15]。由于計(jì)算信息增益率過程中大部分為對(duì)數(shù)運(yùn)算，C4.5運(yùn)算時(shí)間長，效率低。基于麥克勞林公式對(duì)C4.5算法改進(jìn)的決策樹模型，將計(jì)算過程中的對(duì)數(shù)運(yùn)算簡化，提高了運(yùn)算效率[16]。

本文提出一種基于MFDFA和麥克勞林公式改進(jìn)決策樹的電能質(zhì)量分析方法。首先通過分析證明了電能信號(hào)具有多重分形特性，然后計(jì)算信號(hào)的多重分形譜參數(shù)，從中選取3個(gè)參數(shù)(hqmax、αmin、α0)和信號(hào)能量E共4個(gè)參數(shù)作為信號(hào)的特征矩陣，最后，通過改進(jìn)決策樹實(shí)現(xiàn)對(duì)不同類型擾動(dòng)信號(hào)的識(shí)別。

通過實(shí)驗(yàn)及對(duì)比分析驗(yàn)證了MFDFA方法的優(yōu)越性：首先向信號(hào)中添加不同信噪比的白噪聲，證明了該方法具有較好的抗噪性；然后與總體平均經(jīng)驗(yàn)?zāi)Ｊ椒纸?ensemble empirical mode decomposition,EEMD)、希爾伯特-黃變換(Hilbert-Huang transform,HHT)和雙樹復(fù)小波變換(dual-tree complex wavelet transform,DTCWT)三種常見特征提取方法[17]進(jìn)行對(duì)比實(shí)驗(yàn)，用改進(jìn)的決策樹分別對(duì)上述提取的特征向量進(jìn)行分類識(shí)別，結(jié)果證明，用本文方法提取的特征向量進(jìn)行分類其精度優(yōu)于其它方法。

2 電能質(zhì)量信號(hào)的多重分形特征

2.1 多重分形去趨勢波動(dòng)分析

多重分形去趨勢波動(dòng)分析可以有效消除干擾趨勢，估計(jì)多重分形譜，多重分形去趨勢波動(dòng)分析思路是：

(1) 設(shè)非線性時(shí)間序列信號(hào)xk，構(gòu)造信號(hào)的離散時(shí)間序列Yi：

(1)

(2) 將時(shí)間序列Yi劃分為互不重疊的N個(gè)小區(qū)間，每個(gè)小區(qū)間包含s個(gè)數(shù)據(jù)?？紤]到Y(jié)i的長度N不一定整除s，造成剩余一段數(shù)據(jù)剩余，因此再次從數(shù)據(jù)的反方向以相同的長度s將時(shí)間序列分成互不重疊的N段數(shù)據(jù)，這樣一共得到2N段數(shù)據(jù)。

(3) 使用最小二乘法擬合法計(jì)算每段數(shù)據(jù)的局部趨勢，確定每一段數(shù)據(jù)的均方誤差F2(v,s)：

當(dāng)v=1,2,…,N時(shí)，

(2)

當(dāng)v=N+1,…,2N時(shí)，

(3)

式中：yv(i)是第v段數(shù)據(jù)的擬合多項(xiàng)式，多項(xiàng)式y(tǒng)v(i)的階數(shù)m越大趨勢消除的效果越好，相應(yīng)的計(jì)算時(shí)間也會(huì)相應(yīng)增加。當(dāng)m=1，2，3，…時(shí)，分別記為mfdfa1，mfdfa2，mfdfa3，…，通常m+2≤s≤N/4。

(4) 設(shè)波動(dòng)函數(shù)Fq(s)階數(shù)為q，則：

(4)

一般來說q可以取非零的實(shí)數(shù)，當(dāng)q=2時(shí)，MFDFA退化成DFA算法。

(5) 如果離散時(shí)間序列信號(hào)xk存在自相似性，則q階波動(dòng)函數(shù)Fq(s)和時(shí)間尺度s之間滿足冪律關(guān)系：

Fq(s)∝sh(q)

(5)

式中：h(q)為廣義Hurst指數(shù)，定義q～h(q)為廣義Hurst指數(shù)譜，其中：

(6)

h(q)的值隨著q值的改變而發(fā)生變化時(shí)，即q～h(q)表現(xiàn)為非線性函數(shù)，則時(shí)間序列展現(xiàn)出多重分形特征。

對(duì)于離散時(shí)間序列的廣義Hurst指數(shù)h(q)，多重分形理論中的質(zhì)量指數(shù)τ(q)和分形維數(shù)D(q)之間滿足：

τ(q)=qh(q)-1

(7)

(8)

通過Legendre變換，得到奇異指數(shù)α和多重分形譜f(α)：

α=τ(q)′=h(q)+qh(q)′

(9)

f(α)=qα-τ(q)=q[α-h(q)]+1

(10)

多重分形譜寬度Δα=αmax-αmin，反應(yīng)了整個(gè)結(jié)構(gòu)概率分布的不均勻程度，αmax是最大概率子集，αmin是最小概率子集，可以用來描述信號(hào)的多重分形強(qiáng)度。分形維數(shù)差ΔDq=Dqmax-Dqmin，ΔDq反映了最大波動(dòng)與最小波動(dòng)的頻率比，Dqmax是最大波動(dòng)點(diǎn)集的分形維數(shù)，Dqmin是最小波動(dòng)點(diǎn)集的分形維數(shù)，ΔDq值的大小表示電能擾動(dòng)信號(hào)中波動(dòng)最平穩(wěn)的子集個(gè)數(shù)與波動(dòng)最劇烈的子集個(gè)數(shù)的比例。

2.2 電能質(zhì)量信號(hào)的多重分形特征的發(fā)現(xiàn)

本文首先對(duì)電壓暫降(sag)、電壓暫升(swell)、暫態(tài)振蕩(transient oscillation)、諧波(harmonic)、閃變(flicker)、脈沖(transient pulse)等6種電能質(zhì)量信號(hào)進(jìn)行多重分形特征分析，證明電能質(zhì)量信號(hào)具有多重分形特征。根據(jù)國家標(biāo)準(zhǔn)技術(shù)委員會(huì)及相關(guān)行業(yè)標(biāo)委會(huì)制定的電能質(zhì)量國家標(biāo)準(zhǔn)[18]，6種電能質(zhì)量信號(hào)的擾動(dòng)模型及具體參數(shù)設(shè)置如表1所示，信號(hào)擾動(dòng)發(fā)生的幅值、起始時(shí)間和持續(xù)時(shí)間由Matlab隨機(jī)產(chǎn)生。

表1 電能質(zhì)量擾動(dòng)模型Tab.1 Power quality disturbance model

根據(jù)MFDFA算法的原理可以發(fā)現(xiàn)，廣義Hurst指數(shù)h(q)，質(zhì)量指數(shù)τ(q)，奇異指數(shù)α和多重分形譜f(α)四個(gè)參數(shù)可以描述信號(hào)的多重分形特性。

根據(jù)式(6)計(jì)算廣義Hurst指數(shù)h(q)，當(dāng)h(q)的值隨著q值的改變而發(fā)生變化時(shí)，即q～h(q)曲線表現(xiàn)為一個(gè)非線性函數(shù)，則時(shí)間序列展現(xiàn)出多重分形特征。圖1給出了的廣義Hurst指數(shù)，根據(jù)廣義Hurst理論，當(dāng)q>0.5時(shí)，時(shí)間序列時(shí)長程相關(guān)；當(dāng)q<0.5時(shí)，時(shí)間序列時(shí)短程相關(guān)。從圖中可以看出，每一種電能質(zhì)量信號(hào)的h(q)均大于0.5，因此可以確定電能質(zhì)量信號(hào)是一個(gè)長程相關(guān)的時(shí)間序列。

根據(jù)MFDFA算法的原理可以發(fā)現(xiàn)，廣義Hurst指數(shù)h(q)，質(zhì)量指數(shù)τ(q)，奇異指數(shù)α和多重分形譜f(α)四個(gè)參數(shù)可以描述信號(hào)的多重分形特性。

根據(jù)式(6)計(jì)算廣義Hurst指數(shù)h(q)，當(dāng)h(q)的值隨著q值的改變而發(fā)生變化時(shí)，即q～h(q)曲線表現(xiàn)為一個(gè)非線性函數(shù)，則時(shí)間序列展現(xiàn)出多重分形特征。圖1給出了的廣義Hurst指數(shù)，根據(jù)廣義Hurst理論，當(dāng)q>0.5時(shí)，時(shí)間序列時(shí)長程相關(guān)；當(dāng)q<0.5時(shí)，時(shí)間序列時(shí)短程相關(guān)。從圖中可以看出，每一種電能質(zhì)量信號(hào)的h(q)均大于0.5，因此可以確定電能質(zhì)量信號(hào)是一個(gè)長程相關(guān)的時(shí)間序列。

圖1 廣義Hurst指數(shù)h(q)Fig.1 Generalized Hurst index h(q)

根據(jù)式(7)計(jì)算信號(hào)的標(biāo)度指數(shù)τ(q)，判斷多重分形特性。當(dāng)τ(q)為一條直線時(shí)，表明信號(hào)函數(shù)是單分形的，若τ(q)表現(xiàn)出非線性，則表明信號(hào)函數(shù)是多重分形的。對(duì)電能質(zhì)量信號(hào)進(jìn)行分析，得到τ(q)隨q的變化曲線，如圖2所示。從圖中可看出τ(q)～q之間存在明顯的非線性關(guān)系，這表明了電能質(zhì)量信號(hào)具有多重分形性質(zhì)。

圖2 標(biāo)度指數(shù)τ(q)Fig.2 Scale index τ(q)

根據(jù)式(9)和式(10)可以得到信號(hào)的奇異指數(shù)α和多重分形譜f(α)。如果時(shí)間序列是多重分形特性的，α～f(α)曲線是一個(gè)單峰鐘形。圖3給出了6種電能質(zhì)量擾動(dòng)信號(hào)的多重分形譜圖，圖中曲線反映了奇異指數(shù)和分形維數(shù)的關(guān)系。從圖3中，可以明顯看出電能信號(hào)是具有多重分形特征的。

圖3 多重分形譜圖Fig.3 Multifractal spectrum

3 信號(hào)的多重分形特征提取

信號(hào)經(jīng)多重分形分析，計(jì)算出4個(gè)分形參數(shù)Hurst指數(shù)h(q)、尺度指數(shù)τ(q)，奇異指數(shù)α和多重分形譜f(α)，如圖1至圖3所示。從圖中可以看出不同類型的信號(hào)間的h(q)、f(α)和α三個(gè)參數(shù)有著明顯區(qū)別，可以從中選取能表征信號(hào)特征的參數(shù)。下面就以電壓暫降信號(hào)(sag)為例說明基于多重分形參數(shù)的特征提取過程。

圖4和圖5是某一電壓暫降信號(hào)(sag)的多重分形譜f(α)和Hurst指數(shù)h(q)。圖4中α～f(α)曲線中3個(gè)點(diǎn)是起始點(diǎn)，終止點(diǎn)和最高點(diǎn)，它們對(duì)應(yīng)的奇異指數(shù)分別為αmin、αmax、α0。圖5中q～h(q)曲線中3個(gè)點(diǎn)分別是q=-5，q=0，q=5時(shí)的廣義Hurst指數(shù)hqmax、hq0、hqmin。

圖4 sag信號(hào)的多重分形譜Fig.4 Multifractal spectrum of sag signal

圖5 sag信號(hào)的廣義Hurst指數(shù)Fig.5 Generalized Hurst index of sag signal

選取廣義Hurst指數(shù)的最大值hqmax，多重分形譜的奇異指數(shù)的最小值αmin，多重分形譜最高點(diǎn)處對(duì)應(yīng)的奇異指數(shù)α0；另外考慮到電能質(zhì)量信號(hào)發(fā)生擾動(dòng)時(shí)，其幅值大小會(huì)發(fā)生改變，信號(hào)所包含的能量值也會(huì)受到影響，因此可以用信號(hào)的能量值E作為其中一個(gè)特征量。其計(jì)算公式如下。

(11)

4 基于改進(jìn)決策樹的分類識(shí)別方法

C4.5以信息增益率作為屬性選擇判斷條件，同時(shí)增加了同步剪枝的操作，增強(qiáng)了決策樹模型的有效性。但計(jì)算信息增益率過程中大部分為對(duì)數(shù)運(yùn)算，為提高運(yùn)算效率，本文使用麥克勞林公式改進(jìn)的決策樹作為分類器。

4.1 C4.5決策樹算法

對(duì)于一個(gè)隨機(jī)變量S=[S1，S2，…，Si]，其中i=1,2,…,m，假設(shè)X是有m個(gè)不同屬性的樣本集，容量為n，每一類為Ci(i=1，2，…，m)，設(shè)Xi是每一類中的一個(gè)樣本，定義樣本分類所需的信息熵(也稱期望信息)為：

(12)

對(duì)樣本集X按照屬性A劃分為v個(gè)不同的類型，劃分后的新的信息熵為：

(13)

原信息熵與劃分后得到的新信息熵，二者間的差即為信息增益，

Gain(A)=Info(X)-InfoA(X)

(14)

信息增益進(jìn)行規(guī)范化操作所需的參考分裂信息定義為：

(15)

信息增益率的計(jì)算方式為：

(16)

4.2 麥克勞林公式改進(jìn)的C4.5決策樹算法

從C4.5算法原理中不難發(fā)現(xiàn)，求取信息增益率的運(yùn)算包含大量的對(duì)數(shù)計(jì)算，為減少計(jì)算時(shí)間，使用麥克勞林公式改進(jìn)信息增益率的計(jì)算方式，轉(zhuǎn)化成非對(duì)數(shù)的運(yùn)算。式(12)中的信息熵轉(zhuǎn)化為：

(17)

將ln(ni/n)用麥克勞林公式展開，得到：

(18)

式(17)由此可以簡化為：

(19)

此時(shí)，式(13)中的新信息熵和式(15)中的分裂信息的計(jì)算式更新為：

(20)

(21)

改進(jìn)后的信息增益率為：

(22)

式中：Gain′(A)=I′(X)-E′(A)，根據(jù)計(jì)算出信息增益率，生成決策樹模型，利用if-then規(guī)則實(shí)現(xiàn)數(shù)據(jù)的有序劃分。

5 實(shí)驗(yàn)及性能分析

本文采用的數(shù)據(jù)包括電壓暫降、電壓暫升、暫態(tài)振蕩、諧波、閃變、脈沖等6種電能質(zhì)量擾動(dòng)信號(hào)。每種類型均通過MATLAB仿真實(shí)現(xiàn)，每種信號(hào)隨機(jī)仿真100組，采樣頻率 5 000 Hz，信號(hào)長度為8 000點(diǎn)。

5.1 抗噪性能分析

將MFDFA方法分析得到的電能質(zhì)量擾動(dòng)信號(hào)的特征向量矩陣，作為輸入數(shù)據(jù)由決策樹分類器進(jìn)行擾動(dòng)識(shí)別，在隨機(jī)仿真的100組電能質(zhì)量擾動(dòng)信號(hào)中，選取前70組作為訓(xùn)練樣本，后30組作為測試樣本。記電壓暫降為標(biāo)簽1，電壓暫升為標(biāo)簽2，暫態(tài)振蕩為標(biāo)簽3，諧波為標(biāo)簽4，閃變?yōu)闃?biāo)簽5，脈沖為標(biāo)簽6。

為了說明MFDFA方法的抗噪性能，用麥克勞林公式改進(jìn)的C4.5決策樹算法分別對(duì)測試集信號(hào)無噪聲的條件下和加入噪聲條件下的信號(hào)進(jìn)行分類實(shí)驗(yàn)，得到圖6和圖7。圖6是對(duì)無噪聲信號(hào)進(jìn)行分類的結(jié)果，圖7是對(duì)加入不同程度高斯白噪聲的信號(hào)進(jìn)行分類的結(jié)果。

圖6 無噪聲測試集信號(hào)的MFDFA分類結(jié)果Fig.6 Classification results of the MFDFA to the test set signals without noise

圖7 加入噪聲的測試集信號(hào)的MFDFA分類結(jié)果Fig.7 Classification results of the MFDFA to the test set signals added with noises

通過圖6和圖7的效果表明，利用MFDFA方法對(duì)無噪聲測試集信號(hào)分類準(zhǔn)確度達(dá)到100%、對(duì)加入噪聲后的測試集信號(hào)分類準(zhǔn)確度仍然較高(加入10 dB噪聲的分類準(zhǔn)確度為98.89%，加入35 dB噪聲的分類準(zhǔn)確度為96.10%，加入50 dB噪聲的分類準(zhǔn)確度為95.00%)。因此，MFDFA方法對(duì)無噪聲和加噪聲這兩種情形分類準(zhǔn)確性都比較高，證明了本文方法受噪聲影響小，具有較好的抗噪性。

5.2 特征提取實(shí)驗(yàn)對(duì)比分析

使用DTCWT，HHT，EEMD三種常見的時(shí)頻分析方法對(duì)上述電能擾動(dòng)信號(hào)進(jìn)行特征提取。以諧波信號(hào)(har)為例說明，圖8是諧波信號(hào)的五層復(fù)小波分解結(jié)果圖，圖9是對(duì)諧波信號(hào)HHT變換后得到的高中低頻分量，圖10是經(jīng)過EEMD變換后的各層IMF分量。

圖8 諧波信號(hào)DTCWT分析Fig.8 Harmonic signal DTCWT analysis

圖9 諧波信號(hào)的HHT分析Fig.9 Harmonic signal HHT analysis

圖10 諧波信號(hào)的EEMD分解圖Fig.10 Harmonic signal EEMD analysis

采用決策樹對(duì)上述3種方法的分解結(jié)果進(jìn)行分類，其中，決策樹參數(shù)設(shè)置均與上文一致。表1為4種方法的分類準(zhǔn)確度，可以看出由MFDFA方法提取出的特征量其分類精度高于其他3種方法。

從表2和圖11中可以看出，4種特征提取方法中，MFDFA方法的識(shí)別準(zhǔn)確率是最高的，其它3種方法由高到低依次是DTCWT、HHT和EEMD。

圖11 不同方法的誤差曲線Fig.11 Error curve for different methods

表2 不同算法的分類結(jié)果Tab.2 Classification results of different methods

表3給出了4種方法在信噪比為0、10、35、50 dB時(shí)的分類準(zhǔn)確度，可以看出MFDFA方法提取的特征量識(shí)別率、抗噪性均優(yōu)于其他3種方法。

表3 不同信噪比下的分類準(zhǔn)確度Tab.3 Classification accuracy under different SND

6 結(jié) 論

本文首先證明了電能質(zhì)量信號(hào)具有多重分形特征。據(jù)此提出利用多重分形去趨勢波動(dòng)分析方法提取3個(gè)分形參數(shù)，與信號(hào)能量共同作為特征向量，使用改進(jìn)決策樹作為分類器識(shí)別了6種常見的電能質(zhì)量擾動(dòng)信號(hào)。通過仿真實(shí)驗(yàn)驗(yàn)證，該方法在無噪聲的情況下可以精確識(shí)別6種擾動(dòng)信號(hào)，在有噪聲的情況下該方法的分類準(zhǔn)確度可達(dá)95%,表明了該方法具有較好的抗噪性能。通過與DTCWT、HHT和EEMD三種方法進(jìn)行對(duì)比實(shí)驗(yàn)，結(jié)果表明本文所提方法表現(xiàn)出更好的識(shí)別結(jié)果，是一種電能質(zhì)量擾動(dòng)分析有效的新方法。