楊 揚，劉道杰，李 煒

(1.山東石油化工學(xué)院 機(jī)械與控制工程學(xué)院，山東 東營 257061；2.中國石油化工股份有限公司 勝利油田分公司，山東 東營 257061)

0 引言

據(jù)統(tǒng)計，我國停電事故占電力系統(tǒng)故障的95%以上[1-2]，各類故障發(fā)生時，故障選線相對困難。近年來，電網(wǎng)信息管理系統(tǒng)大規(guī)模推廣應(yīng)用，可實現(xiàn)用電信息采集、配電自動化等功能。與此同時，隨著系統(tǒng)長時間運行，產(chǎn)生了運檢、調(diào)度、營銷等方面的大量多源數(shù)據(jù)，這些海量數(shù)據(jù)為研究電網(wǎng)故障提供了有利條件。伴隨電網(wǎng)數(shù)據(jù)多維度、強(qiáng)耦合、高隨機(jī)性的特征，計算機(jī)終端數(shù)據(jù)異常龐大，計算工作呈現(xiàn)指數(shù)增長趨勢。導(dǎo)致在大電網(wǎng)系統(tǒng)中出現(xiàn)網(wǎng)損增大、電壓質(zhì)量變差等問題，影響到系統(tǒng)安全平穩(wěn)運行，嚴(yán)重時會引起電壓崩潰造成大面積的停電事故。因此，如何有效利用、處理這些海量數(shù)據(jù)，保障電網(wǎng)安全穩(wěn)定運行，是當(dāng)前電力大數(shù)據(jù)時代下亟待解決的問題[3]。

近年來，國家電網(wǎng)持續(xù)加大對電力大數(shù)據(jù)的研究與應(yīng)用，出現(xiàn)了大量電力調(diào)度自動化MIS(管理信息系統(tǒng))，伴隨產(chǎn)生大量的動態(tài)數(shù)據(jù)，推動了對大數(shù)據(jù)時代電網(wǎng)故障檢測和分析的研究。由于高斯噪聲、測量誤差等因素造成的數(shù)據(jù)具有隨機(jī)性，采用隨機(jī)矩陣?yán)碚摽蓪﹄娋W(wǎng)進(jìn)行有效的分析。如文獻(xiàn)[4]通過對數(shù)據(jù)流中的變異點進(jìn)行檢測，借助樣本協(xié)方差矩陣最大特征值對低信噪比環(huán)境電網(wǎng)異常狀態(tài)進(jìn)行檢測。文獻(xiàn)[5]利用隨機(jī)矩陣?yán)碚搶﹄娋W(wǎng)負(fù)荷曲線中的異常點進(jìn)行了檢測和修正。文獻(xiàn)[6]利用數(shù)據(jù)源和失效模式的映射，通過處理大量數(shù)據(jù)粗糙集信息熵來簡化信息。但是隨著電力系統(tǒng)規(guī)模的不斷擴(kuò)大，需要假設(shè)和簡化構(gòu)建模型和參數(shù)，計算量大且精確度不高，難以實現(xiàn)電網(wǎng)實時檢測分析，同時特征數(shù)據(jù)具有不同類型和結(jié)構(gòu)，這給數(shù)據(jù)存儲、運維策略帶來很多問題。應(yīng)用隨機(jī)矩陣原理對配電網(wǎng)絡(luò)故障選線中多源異常數(shù)據(jù)的處理，在不需要假設(shè)、簡化的情況下，可以減少分析過程中的偶然性、不確定性，該方法能夠準(zhǔn)確、實時地獲得分析結(jié)果，從而實現(xiàn)了配電網(wǎng)的故障選線。在大數(shù)據(jù)的基礎(chǔ)上，通過對電力系統(tǒng)的歷史數(shù)據(jù)的分析，可以使電網(wǎng)的故障選線更加有效。

本文結(jié)合單環(huán)定理、特征值譜分析等隨機(jī)矩陣?yán)碚搶﹄娋W(wǎng)故障進(jìn)行檢測，提出了一種結(jié)合平均譜半徑(MSR)理論和M-P定律雙重判據(jù)定位方法。將數(shù)據(jù)的統(tǒng)計學(xué)特征與故障選線策略相對應(yīng)，為大數(shù)據(jù)構(gòu)架解決電網(wǎng)故障檢測提供了新思路。本文通過IEEE39節(jié)點標(biāo)準(zhǔn)系統(tǒng)和某油田實際電網(wǎng)兩個算例進(jìn)行驗證，結(jié)果表明了該方法的準(zhǔn)確性和有效性。

1 高維隨機(jī)矩陣相關(guān)理論

1.1 M-P定律

1976年Marcenko和Pastur提出并證明了M-P法則，用于描述任意矩陣的特征值，其矩陣經(jīng)驗頻譜分布滿足某種規(guī)則[7]。

在矩陣分析中，一般采用經(jīng)驗頻譜分布函數(shù)，一個m階的方陣U，它的經(jīng)驗譜分布函數(shù)為：

(1)

對高維隨機(jī)矩陣Xmxn，將其各列單獨提取出來，形成N個維數(shù)為M的向量，從這M個向量中任意選一個向量：

xa=[x1a，x2a，x3a，…，x(m-1)axma]

(2)

則Xmxn的樣本協(xié)方差矩陣為：

(3)

一般地，高維隨機(jī)矩陣Xmxn的樣本協(xié)方差矩陣還可以做如下定義：

(4)

隨機(jī)矩陣Xmxn的協(xié)方差矩陣Sn，其特征值λi收斂于函數(shù)F(λi)，滿足F(λi)極限譜密度函數(shù)分布：

(5)

通過M-P定律，可以很好地估算出隨機(jī)矩陣中各個要素的標(biāo)準(zhǔn)偏差，并與特征值的頻譜密度作比較，使其直觀地反映出這些偏差程度。

1.2 單環(huán)定理

在高維隨機(jī)理論中，通過對矩陣的經(jīng)驗特征值分配概率來量化結(jié)果。在不受干擾的情況下，它的特征值在復(fù)平面上向一個圓環(huán)聚斂。相反在某一特定因素作用下，隨機(jī)矩陣的特征值分布將不再聚斂。步驟如下：

隨機(jī)矩陣Xmxn做歸一化處理后成為，得到標(biāo)準(zhǔn)的非Hermitian隨機(jī)矩陣：

(6)

式中，μ(xi)，σ(yi)分別為隨機(jī)矩陣Xmxn的第i列均值和標(biāo)準(zhǔn)差2]滿足：

(7)

非Hermitian隨機(jī)矩陣Ymxn可以按照公式(8)求其奇異值等價矩陣Yu，i；

(8)

式中，U為Haar酉矩陣，并且有：

(9)

奇異值等價矩陣Yu，i是一個高維隨機(jī)矩陣，其滿足公式(10)性質(zhì)：

Yu，i=UnΛnVn

(10)

其中：Λn是一個復(fù)對角陣；UnVn是兩個互不關(guān)聯(lián)的隨機(jī)酉矩陣，并且嚴(yán)格服從Haar分布。

假設(shè)Yu，i有W個奇異值等價矩陣，那么求所有等價矩陣積Z：

(11)

(12)

式中，i=1，2，…，W；σ(zi)為Z的標(biāo)準(zhǔn)差。

(13)

單環(huán)定理能將采樣矩陣的特征降維為奇異值的頻譜分布，并具有成熟的可視化方法，從而實現(xiàn)對測量數(shù)據(jù)的異常情況地直接比對。圖1采用Matlab中plot函數(shù)繪制示意圖，圖中細(xì)實線代表r=1的外環(huán)，細(xì)虛線代表長度為(1-δ)L/2的內(nèi)環(huán)，※號代表的是λ2，◎代表的是λ1。電網(wǎng)在非異常狀態(tài)下，對差分、規(guī)范化的采樣矩陣Xmxn進(jìn)行歸一化，所得的矩陣單元近似為獨立高斯分布。

圖1 單環(huán)定理示意圖

1.3 平均譜半徑(MSR)理論

采用隨機(jī)矩陣法來分析特定事件時，用平均譜半徑來判定(MSR)[8]：

(14)

式中，rMSR為矩陣的譜平均半徑，λi為矩陣的M個特征值。圖2為IEEE39節(jié)點9號發(fā)生故障時MSR計算結(jié)果。由于采用滑動時窗的原因，從第200個數(shù)據(jù)開始顯示，正常狀態(tài)下的MSR分布在內(nèi)環(huán)和外環(huán)之間，當(dāng)故障發(fā)生時會跌落至內(nèi)環(huán)以下。如圖2所示，設(shè)置5 s發(fā)生單相接地故障，MSR小于內(nèi)環(huán)半徑，表明系統(tǒng)異常工作。MSR的波動性和失效時間都與模擬值一致，失效持續(xù)時間與Tw相關(guān)，MSR實現(xiàn)了從宏觀角度來判定系統(tǒng)工作狀態(tài)。

1.4 樣本協(xié)方差矩陣?yán)碚?/h3>

電網(wǎng)中選取n個節(jié)點量測數(shù)據(jù)作為空間樣本，其中每個節(jié)點共有m個狀態(tài)變量，構(gòu)成N個變量，即N=n*m。在任意采樣時刻Ti，所有節(jié)點量測數(shù)據(jù)構(gòu)成一個列向量：

x(ti)=[x1i，x2i，…，xni]H

(15)

在一個時間段內(nèi)，N個向量在時序上矩陣分布為公式(16)：

X=[x(t1)，x(t2)，…，x(ti)，…]

(16)

為了可以對電力網(wǎng)絡(luò)數(shù)據(jù)進(jìn)行實時分析，采用了滑動時間窗[10-12](圖3)加以輔助，實時窗口長度為Mw，寬度為Nw。當(dāng)數(shù)據(jù)都是正常值的時候，生成的矩陣是滿足某種分布規(guī)律，而遇到異常數(shù)據(jù)出現(xiàn)時，時窗生成的矩陣的奇異性就會變化。采樣時向后平移一個采樣點，因為前Tw寬度數(shù)據(jù)沒有，所以只能檢測之后發(fā)生的故障，移動時窗每次僅對Tw寬度數(shù)據(jù)處理，可以減少歷史處理對當(dāng)前數(shù)據(jù)狀態(tài)的影響(具體分析見3.1)。ti時刻的實時滑動時間窗包含ti時刻的當(dāng)前數(shù)據(jù)和(Nw-1)個相鄰的歷史數(shù)據(jù)，即

圖3 滑動時間窗口

(17)

在狀態(tài)評估時間窗內(nèi)，對系統(tǒng)的M個可測狀態(tài)參量測量N次，所有測量數(shù)據(jù)組成的具有電網(wǎng)時空特性的M行N列數(shù)據(jù)：

(18)

2 電網(wǎng)故障檢測步驟與閾值設(shè)定

本算法故障檢測分兩部分，一部分是確定故障時間的(如圖4)，第二部分為判別故障饋線及故障類型(如圖5)。

圖4 確定故障時間流程圖

圖5 確定故障饋線與類型流程圖

2.1 故障檢測步驟

2.1.1 確定故障時間

從數(shù)據(jù)源采集數(shù)據(jù)后，采用滑動時間窗口處理數(shù)據(jù)，融合歷史數(shù)據(jù)和當(dāng)前數(shù)據(jù)作為量測時間序列D，對矩陣進(jìn)行差分、加躁、歸一化預(yù)處理生成Haar矩陣U，差分的目的是使樣本數(shù)據(jù)平穩(wěn)化，若一次差分后仍未符合要求，則采取二次差分，一般二次差分即可滿足平穩(wěn)要求。

2.1.2 確定故障類型

通過1)得到故障時間后，選取故障時間附近時窗寬度的三相電壓數(shù)據(jù)作為數(shù)據(jù)源同樣采用滑動時窗的方法得到時間序列矩陣，并對其預(yù)處理，得到協(xié)方差矩陣，將所有節(jié)點的ABC三相各分為一組，通過計算頻譜偏差比較超過閾值的大小確定故障饋線，通過匯集所有節(jié)點故障相后綜合分析得到故障類型。若某節(jié)點僅有一相超過閾值，則判定為單相故障；若某節(jié)點有兩相超過閾值，則判定為相間故障；同理，若三相均超過閾值，則判定為三相故障[9]。

2.2 閾值設(shè)定

電網(wǎng)出現(xiàn)故障時，矩陣譜的變化對矩陣元素變化的影響最大，特別是最大特征值和最小特征值，因此定義t時刻的頻譜偏差ds(t)為：

(19)

通過對各節(jié)點位置分析矩陣的頻差進(jìn)行對比，對電網(wǎng)故障進(jìn)行準(zhǔn)確定位[13-15]。在有大量節(jié)點的情況下，相鄰的結(jié)點往往具有近似的頻譜偏差，很難找出最大頻譜偏差的結(jié)點。為了準(zhǔn)確地識別故障節(jié)點，本文對節(jié)點頻譜偏差進(jìn)行了以下改進(jìn)：

(20)

式中，dts節(jié)點t改善了原指標(biāo)數(shù)據(jù)輸出，dmax為各結(jié)點最大偏差值，ddif1為最大值與次最大值的差值。由上式可知，對于任意一個節(jié)點t，均有dmax-dt≧ddif1，進(jìn)一步可得dts≤2/ddif1。同時，相對原偏差較大的節(jié)點和原偏差值幅度較為接近，最大與第二大差值為較小的實數(shù)，經(jīng)上述公式處理后，偏差可被判別。

3 算例分析

3.1 算例1—IEEE39模型分析

算例1基于圖6新英格蘭IEEE39節(jié)點模型，仿真總時長T=10 s，系統(tǒng)運行頻率為60 Hz，數(shù)據(jù)采樣頻率也為60 Hz，選擇每個節(jié)點及線路測量儀(117個)三相電壓電流功率等k個采樣參數(shù)構(gòu)成數(shù)據(jù)源矩陣U，設(shè)置時窗寬度T_windows=200。

圖6 IEEE39節(jié)點系統(tǒng)

3.1.1 異常狀態(tài)判定

以IEEE39節(jié)點系統(tǒng)節(jié)點負(fù)荷增加作為異常狀態(tài)的判定，區(qū)別于正常狀態(tài)和故障狀態(tài)的電網(wǎng)第三種狀態(tài)異常狀態(tài)。本案例以節(jié)點負(fù)荷變化驗證算法判斷的準(zhǔn)確性。設(shè)置節(jié)點20設(shè)置5 s時從7.5 MW變化到30 MW，持續(xù)時間0.25 s，研究判據(jù)對異常時刻的準(zhǔn)確性[3]。設(shè)置的負(fù)荷變化節(jié)點及變化幅度如表1所示。

表1 負(fù)荷變化表

由圖7(a)可知譜平均半徑大于內(nèi)環(huán)半徑附近，未跌落至內(nèi)環(huán)以內(nèi)，區(qū)別于故障狀態(tài)。從頻譜偏差的圖7(b)中看出在5 s時出現(xiàn)大幅升高，說明出現(xiàn)非正常狀態(tài)。

圖7 頻譜偏差圖

通過比較圓環(huán)分布圖8(a)和(b)可知，正常狀態(tài)時奇異值分布基本在圓環(huán)里面，異常狀態(tài)奇異值有部分分布在環(huán)外[3]。

圖8 狀態(tài)圖

3.1.2 故障時間確定

算例1用MSR和頻譜偏差分析運行整體態(tài)勢，采用滑動時間窗口選取含有故障時間所有節(jié)點的三相電壓作為數(shù)據(jù)源矩陣進(jìn)行具體分析[7]，并用M-P定律和單環(huán)定理來雙重判斷故障時刻，驗證判定方法是否在電網(wǎng)系統(tǒng)運行中有效性。

1)相同節(jié)點發(fā)生不同故障：

設(shè)置節(jié)點4處在t=5 s時分別發(fā)生單相、相間、相間接地、三相故障，故障持續(xù)0.25 s。由于算例中滑動時間窗口采樣點數(shù)T=200，故系統(tǒng)運行狀態(tài)監(jiān)控測量的開始時間為200/60 s，0到(200/60)s無法檢測出故障不顯示數(shù)據(jù)，故譜平均半徑圖像開始的顯示時刻為第200/60時刻。

從隨機(jī)矩陣法和MSR的定義中可以看出，在MSR小于內(nèi)環(huán)時，系統(tǒng)的工作是異常的。在圖9、圖10中，MSR和譜偏離度的波動性和異常時間都與提前設(shè)定值一致。其中，異常狀態(tài)持續(xù)時間和Tw相關(guān)，還可以推斷出異常的恢復(fù)時間。MSR值越低，譜偏離度越高，表明其與電網(wǎng)的運行相關(guān)性越大。

圖10 節(jié)點9發(fā)生單相故障時的頻譜偏差

如圖10所示，正常時刻的特征值基本分布在圓環(huán)內(nèi)，其采樣時間窗口內(nèi)數(shù)據(jù)的譜平均半徑為0.840 4，在內(nèi)環(huán)和外環(huán)中間。5 s單相故障時刻的特征值基本分布在內(nèi)環(huán)外側(cè)，當(dāng)5.017 s時檢測到故障發(fā)生，其單相、相間、三相故障的譜平均半徑(取穩(wěn)定時段下的數(shù)值)分別為0.275 06、0.129 18、0.101 94均小于內(nèi)環(huán)半徑。當(dāng)采用譜平均半徑和頻譜偏差檢測到故障時間后，為進(jìn)一步驗證，作出異常數(shù)據(jù)和正常數(shù)據(jù)的單環(huán)定理圖進(jìn)行對比，利用M-P定律和單環(huán)定理來雙重驗證準(zhǔn)確性。

圖11為故障單環(huán)圖，從特征點的分布集中在圓環(huán)內(nèi)可以識別出異常數(shù)據(jù)。圖12 M-P律圖對比看出正常數(shù)據(jù)和異常數(shù)據(jù)的區(qū)別，異常數(shù)據(jù)的核密度曲線與極限譜密度函數(shù)分布偏離了理想狀態(tài)。

圖11 故障單環(huán)圖

圖12 故障M-P律圖

2)不同節(jié)點發(fā)生相同故障：

圖13表示不同節(jié)點發(fā)生不同故障的MSR歸納圖，節(jié)點9和節(jié)點18在5 s分別設(shè)置持續(xù)時間為0.25 s的單相、相間、相間接地、三相故障，故障程度越高，MSR跌幅越大。MSR數(shù)值由大到小的順序依次為單相接地〉相間故障〉相間接地〉三相接地，故障程度由大到小的順序依次為三相故障〉相間接地〉相間故障〉單相故障，所以該方法可以實現(xiàn)不同故障判別。

圖13 不同節(jié)點發(fā)生不同故障的MSR圖

如果出現(xiàn)短路故障，則多個測點的數(shù)據(jù)都會發(fā)生同樣趨勢變化，矩陣的相關(guān)性會增強(qiáng)，不能滿足獨立的高斯分布。同一故障在多個結(jié)點出現(xiàn)時，所產(chǎn)生的關(guān)聯(lián)程度也不一樣，因而特征值向中心聚集的程度也不一樣，異常的程度與系統(tǒng)內(nèi)部的關(guān)聯(lián)度呈正相關(guān)，與特征值向中心聚集的程度成正比。

3.1.3 故障類型識別

利用IEEE39模型對單相故障、相間故障、三相故障的仿真分析各節(jié)點三相電流，分別分為A相、B相、C相電流采用改進(jìn)后的頻譜偏差來測試故障類型的準(zhǔn)確性。由閾值設(shè)定公式頻譜偏差為5.2，得節(jié)點4發(fā)生各種故障時的三相電壓頻譜偏差。當(dāng)檢測到故障后，進(jìn)行故障定位和故障類型判別，設(shè)定故障為節(jié)點4分別設(shè)定單相故相間故障、三相故障，設(shè)置故障時間5 s，持續(xù)時間0.25 s，采集滑動時窗長度的包含故障時間所有節(jié)點的三相電流并計算頻譜偏差dRa、dRb、dRc，如圖14(a)、(b)、(c)中，橫坐標(biāo)為節(jié)點編號，縱坐標(biāo)為三相電流頻譜偏差。

3.2 算例2—A油田電網(wǎng)模型分析

以A油田某年1～6月期間的全部線路電壓監(jiān)測數(shù)據(jù)為例，每隔1小時即可進(jìn)行數(shù)據(jù)采集。應(yīng)用高維隨機(jī)矩陣?yán)碚摚葟臅r間維度上計算出頻譜偏差和譜平均半徑，初步確定故障時間，在經(jīng)過M-P定律和單環(huán)定理雙重判定故障時間，選取故障時間附近的時窗寬度數(shù)據(jù)再構(gòu)成數(shù)據(jù)源，對頻譜的偏離進(jìn)行了空間維度的分析，找到異常數(shù)據(jù)的來源，確定故障饋線。

3.2.1 確定故障時刻

首先采用滑動時間窗口生成時序矩陣，再經(jīng)過差分、加躁、歸一化處理得到協(xié)方差矩陣，最后求取特征值、譜平均半徑和譜偏離度，如圖15、16所示。在圈中對應(yīng)時刻，頻譜偏差和MSR都是異常的，偏離正常指標(biāo)。

圖15 譜平均半徑

結(jié)果與算例1理想模型對比差異較大，可以看出實際生產(chǎn)中與理論算法的差異性，使得在應(yīng)用到實際環(huán)境中出現(xiàn)不適應(yīng)現(xiàn)象。但是故障時刻相較于正常狀態(tài)區(qū)別較大，相對明顯的MSR下跌仍然有較好的判斷能力，故標(biāo)準(zhǔn)算例及實際系統(tǒng)均能較好地看出本算法的可靠性，能精確的確定故障時間。圖16中表示1 080 h和3 800 h為故障時刻，進(jìn)一步使用M-P和單環(huán)定理雙重判據(jù)，每個隨機(jī)矩陣的特征值特性如圖17(a)和(b)所示。

圖16 頻譜偏差

圖17 譜密度函數(shù)與M-P定律比較

對比圖17發(fā)現(xiàn)，當(dāng)特征值增加時，正常狀態(tài)的核密度估算曲線的概率密度會先增加后降低，曲線的重合度較故障狀態(tài)更大。故障狀態(tài)的核密度估計曲線則與M-P定律重合度小，核密度估計曲線隨著特征值的增大而呈現(xiàn)出先升后降的趨勢，與正常狀態(tài)時有很大差異。由此看出，故障狀態(tài)時電壓、電流數(shù)據(jù)與正常狀態(tài)時電壓、電流數(shù)據(jù)存在很大差別，矩陣不再具有獨立的分布特征，從而判定在3 080 h系統(tǒng)發(fā)生故障。

3.2.2 故障位置判定

以3 080 h異常時間為例，對頻譜的偏差進(jìn)行了空間維度的分析，并找到了異常數(shù)據(jù)的來源。在圖18中給出了3 080 h的頻譜偏差，其中Ⅰ線、Ⅵ線頻譜偏差程度與其他線有較大偏差，故判定故障饋線。

圖18 電流數(shù)據(jù)的頻譜偏差(空間維度)

4 結(jié)束語

根據(jù)隨機(jī)矩陣相關(guān)理論進(jìn)行電網(wǎng)異常分析、線路故障時刻和故障類型判定及故障選線。通過應(yīng)用兩個算例數(shù)據(jù)源建立隨機(jī)矩陣的方法，并對其進(jìn)行了分析驗證算法。在技術(shù)上，它避免了傳統(tǒng)的人工拉路方式，不能適應(yīng)配網(wǎng)故障后的自動選線。通過對故障選線的研究，對不同故障線路所產(chǎn)生的隨機(jī)矩陣的特性分布進(jìn)行了分析，利用M-P律和核密度曲線重合度和單環(huán)原理對故障線路進(jìn)行快速、準(zhǔn)確的選擇。

1)采集三相電壓電流數(shù)據(jù)源構(gòu)建高維隨機(jī)矩陣，采用滑動時間窗口生成時間序列并對其進(jìn)行差分、加噪、歸一化生成協(xié)方差矩陣，求特征值和特征向量。利用譜偏離度、譜平均半徑等指標(biāo)對電網(wǎng)線路進(jìn)行分析。

2)提出一種MSR和M-P定律雙重判據(jù)定位方法，通過IEEE39節(jié)點標(biāo)準(zhǔn)系統(tǒng)及某油田電網(wǎng)某年1～6月份的全部線路電壓、電流數(shù)據(jù)兩個算例，通過比較結(jié)果表明所提方法的準(zhǔn)確性和有效性。

本文應(yīng)用隨機(jī)矩陣法對電網(wǎng)異常和故障狀態(tài)進(jìn)行研究，根據(jù)隨機(jī)矩陣極限譜分布的原理，給出了頻譜偏離程度和譜平均半徑的定量計算方法，相比其他算法，該算法的優(yōu)勢在于譜偏離度的大小之前小了很多，從而大大擴(kuò)大了故障檢測的范圍，提高了故障識別的效率。