基于樣本協方差矩陣特征特性的電網多重擾動定位方法

李洪乾 韓 松 周忠強

（1. 貴州大學電氣工程學院 貴陽 550025 2. 貴州電網有限公司電力調度控制中心 貴陽 550002）

0 引言

隨著同步相量測量裝置、廣域測量系統（Wide Area Measurement System, WAMS）的廣泛應用以及5G等信息通信技術的快速演進，面向電網狀態的隨機波動性和復雜耦合性[1]，從數據驅動角度出發，借助大數據技術[2]，將電網量測系統采集的海量電力內外部多源數據進行知識提取和價值應用，或為新一代電力系統安全分析與控制帶來新的解決方案和技術路徑[3-4]。

在大數據分析方法中，隨機矩陣理論（Random Matrix Theory, RMT）以大維統計原理作為根基，從宏觀角度反映系統當前運行狀態[5]。一方面，在基于RMT的電力系統分析理論與方法研究中，平均譜半徑（Mean Spectral Radius, MSR）指標已初步實現對電網異常場景識別[6]、用戶竊電偵測[7]和電網擾動傳播分析[8]等應用。相較MSR指標，樣本協方差矩陣最大特征值（Maximum Eigenvalue of Sample Covariance Matrix，Max-ESCM）指標不僅具備相近功能，且更適用于低信噪比場景[9-10]，在大規模電網中計算效率更高[11]。而基于RMT中的樣本協方差矩陣最小特征向量（Minimum Eigenvector of Sample Covariance Matrix, Min-ESCM）指標[12]雖然在數學領域已進行了研究，但目前在電力系統領域還少有探索。另一方面，在基于RMT的電力系統領域應用探索過程中，對電網擾動定位已進行了初步研究。文獻[13]基于RMT并結合增廣矩陣法提出了一種故障時刻確定和故障區域定位方法，通過分別計算及比較各節點 MSR指標之差，實現了對故障節點的定位。在文獻[13]的基礎上，文獻[14]通過分別構建每條線路的增廣矩陣，并對其 MSR指標進行計算及越限判別，從而定位出故障線路。不同于文獻[13-14]，文獻[15]基于RMT通過不斷對電網矩陣進行分塊計算并在圖形界面中進行對比，實現了對電網擾動區域的定位。然而，隨著系統規模的增大，上述文獻所需運算次數將大幅增加，導致計算效率偏低。此外，這些文獻中，擾動節點是通過指標間比較的方式來定位的。因此，當電網發生多重擾動[16]時，存在難以對其進行有效定位、適應性差的問題。

鑒于此，為了提升 RMT在電網擾動定位應用中的計算效率和適應性，本文基于 Max-ESCM 及Min-ESCM提出了一種適用于多重擾動場景的電網擾動定位方法。借助DIgSILENT和Matlab R2014a軟件，通過IEEE 54機118節點系統算例，表明了該方法的有效性和高效性。

1 基于Spiked模型的電網擾動識別

1.1 數據源矩陣的構建與預處理

在廣域測量系統中，相量測量單元（Phasor Measurement Unit, PMU）能采集海量具有統一時間戳的電氣量測數據[6]。將篩選后的量測數據按照時空特性排列，可構成一個二維矩陣，即數據源矩陣Xs，如式（1）所示。

Xs不僅存在隨機噪聲的干擾，還受到負荷隨機波動造成的影響，其檢測模型可寫為

式中，Xp為未受噪聲污染的信號矩陣；ψ為負荷隨機波動率，波動范圍設置為±1%；η為噪聲矩陣；m為噪聲幅值。

得到Xs后，對其采用滑動窗口技術可得N×T維窗口數據矩陣X，并將X的行向量按式（3）進行標準化處理，得到標準非Hermitian矩陣Xn。

式中，xi=(xi,1,xi,2,…,xi,T)；μ(xi)、σ(xi)分別為xi的均值和標準差；μ(xni)、σ(xni)分別為非 Hermitian矩陣行向量xni的均值和標準差。

1.2 基于Mar?henko-Pastur定律的擾動識別

RMT中 M-P定律（Mar?henko-Pastur law）[17]對大維樣本協方差矩陣的漸進譜分布特性進行了描述。采用M-P定律可對經過預處理后的電網量測數據矩陣進行分析。當電網中無擾動事件發生時，各節點采集的量測數據整體相對平穩且具有一定隨機性，此時其 Max-ESCM指標也將滿足其統計規律，收斂于一定范圍。但若電網中出現負荷突變或線路故障等擾動情況，則采集的量測數據將發生突變而且系統的隨機性也會被破壞，導致Max-ESCM指標越過其正常統計邊界。因此，可以利用該定律識別電網中擾動事件的發生。對于M-P定律[17]的原理與方法描述如下。

設Xn={xi,j}1≤i≤N,1≤j≤T為一個N×T維的隨機矩陣，且每個元素滿足獨立同分布。當μ(x)=0，σ2(x)＜∞時，則Xn的N×N維樣本協方差矩陣S為

式中，上標H表示共軛轉置。此時，樣本協方差矩陣S也為實對稱矩陣，當維容比c=N/T∈(0,∞)不變時，S的經驗譜密度將收斂于M-P律，譜密度函數可表示為

式中，σ2為刻度參數，σ2=1；a和b分別為譜密度函數中特征值的理論下、上確界。

1.3 基于Spiked模型的Max-ESCM動態閾值

傳統M-P律中，閾值采用的是式（6）中的邊界值，但該邊界值a與b僅與矩陣的c有關，而忽略了樣本所受到的干擾，故該邊界值存在一定的局限性。因此，為了增強對擾動事件判定的適應性，本文采用前期研究中基于 Spiked模型的動態閾值[9]，對 Max-ESCM指標進行越限判別。該動態閾值在傳統閾值模型的基礎上，進一步考慮了樣本所受干擾的影響，增強了擾動識別的適應性，其具體表達式為

式中，γλ為Max-ESCM指標的動態閾值；ρ為當前時刻全局信噪比估算值；c為維容比；α為比例系數，0≤α≤1，可根據滑動窗寬度T進行調整[9]，一般取α=0.5。

2 基于“相變”現象的電網多重擾動定位

2.1 多重擾動

電力系統是復雜的非線性時變系統，在運行過程中會不可避免地受到一系列擾動，使其結構發生改變、潮流重新分配，甚至導致系統失穩[18]。雖然電網中發生多重擾動的概率較低，但若僅針對電網單一擾動進行分析是不全面的。因此，從多重擾動觀點出發研究電網擾動定位是有必要的，并且這也更能反映實際電網擾動事件的隨機性和突發性。

電網中多重擾動是由引起系統運行狀態變化的多個事件在時間、空間、類型上組合而成[16]。根據時間先后發生順序可分為一個節點發生擾動，另一個節點也在相隔較小的時間內發生擾動的同時多重擾動，以及一個節點擾動形成后，影響或傳播到其他節點也發生擾動的相繼多重擾動。

2.2 基于Min-ESCM的多重擾動定位

識別出電網有擾動事件發生后，對所有擾動源進行快速定位對電網安全運行具有重要意義。區別于傳統模型驅動的電網擾動定位方法，基于 Min-ESCM的電網多重擾動定位方法是一種數據驅動的方法，其具體原理如下。

對于任意一個樣本協方差矩陣S，其主對角線呈自相關，副對角線呈互相關。進行特征分解可得

式中，v為矩陣S的特征向量；λ為矩陣S的特征值。

根據定理[19]：設A=(aij)N×N是非負不可約矩陣，v=[v1v2…vn]T是一個正向量，令Xij=aijvj/vi(Xij≥0)，若則λ是A的模最大特征值，結合式（8）可得[20]

由式（10）可知，在矩陣S的vmin中，vmin的第i個元素vmini的變化主要與矩陣S的第i行有關聯。當矩陣S的第i行出現波動時，將導致vmin上第i個元素的值出現顯著變化，且與其他元素有明顯區別。因此，對于由電網量測數據處理得到的樣本協方差矩陣，當電網中出現擾動時，擾動節點k、l…的數據將出現波動，這時將導致Min-ESCM指標中第k、l…號元素的數值相較于其他元素出現顯著變化。故可依據該原理對電網多重擾動事件進行快速定位。

2.3 基于“相變”現象的Min-ESCM動態閾值

“相變”現象原指物理系統中的臨界現象，當物質達到某臨界點時，將會從一種狀態轉換到另外一種狀態。而對于大維實對稱矩陣的特征向量，當對應特征值達到一定范圍時，該特征向量中各元素將呈現出特定的分布規律，在數學領域上稱其出現了“相變”現象[12,21-23]。

假設c∈(0, 1]，有一N×1維單位向量ek，其中僅第k個元素的值為1。當矩陣S的特征值處于不同范圍時，對應特征向量將收斂于[12]

式中，γv為Min-ESCM指標中每一個元素的閾值，當Min-ESCM指標中第k、l…號元素的值大于γv時，說明系統受到的擾動是由第k、l…節點引起的；β為比例系數，0≤β≤1，可根據滑動窗寬度T進行調整[9]，一般取0.5。

3 所提方法流程圖與步驟

電網多重擾動定位方法流程如圖1所示。具體步驟如下：

1）由式（1），將PMU采集的量測數據構造為數據源矩陣Xs。如為模擬現場信號，可通過式（2）引入噪聲和隨機波動負荷。

2）根據式（3），對矩陣Xs進行預處理后，結合式（4）獲取其N×N維樣本協方差矩陣S。

3）計算S的最大特征值λmax并將其作為擾動識別指標Max-ESCM。

4）由式（7）計算當前時刻 Max-ESCM 基于Spiked模型的動態閾值γλ。

5）判斷λmax≥γλ是否成立，若成立，則判定電網有擾動事件發生，需進一步對其進行定位并執行步驟6），否則重復步驟2）～步驟5）。

圖1 所提電網多重擾動定位方法流程Fig.1 The flowchart of the proposed multiple disturbance positioning method for power system

6）計算擾動定位判別指標Min-ESCM即Vmin。

7）根據式（12）計算Vmin基于“相變”現象的動態閾值γv。并設置初始節點k=1。

8）判斷Vk≥γv是否成立，其中Vk為Vmin中第k個元素，若成立，記錄引起擾動事件的節點k，并進行步驟9），否則直接進行步驟9）。

9）判斷k＞N是否成立，若成立，輸出所有擾動節點編號，否則k=k+1并回到步驟8）。

4 算例分析

借助DIgSILENT軟件，對如圖2所示的IEEE 54機118節點系統[24]開展時域仿真獲取測試數據，其中仿真步長均為Δt=0.01s。依據第3節方法步驟在 Matlab R2014a軟件中編制算法程序，以驗證本文所提多重擾動定位方法的有效性。

4.1 同時多重擾動定位測試

為模擬電網中同時多重擾動場景，設置節點21與相隔較遠的節點 95在t1001采樣時刻同時出現負荷躍變，具體見表1。

表1 一類同時躍變的合成負荷Tab.1 The synthetic loads with simultaneous step-up change

圖2 IEEE118節點測試系統Fig.2 IEEE 118-bus test system

按照第3節步驟1），選取全網118個節點的電壓幅值作為數據源進行測試。為模擬PMU量測數據中噪聲干擾和隨機負荷波動，在該信號中引入高斯噪聲源[10]。這樣可由式（1）和式（2）構建一個118維數據源矩陣進行分析。由步驟 2），設滑動窗口寬度T為240，則Max-ESCM指標變化曲線從采樣時刻t240開始，此時維容比c=0.5[13]。將118×240維窗口矩陣按式（3）標準化后，根據式（4）獲取其樣本協方差矩陣。進一步地，由步驟3）和步驟4），可分別獲得當前時刻Max-ESCM指標及其對應動態閾值。

依次對不同時刻滑動窗口數據進行計算，得到t240～t2000采樣時刻的Max-ESCM指標及其對應動態閾值曲線，如圖3所示。由圖3可知，在t1000采樣時刻之前，雖然受到了隨機負荷的波動及噪聲的干擾，但Max-ESCM指標無明顯變化。然而在t1001采樣時刻，圖3中Max-ESCM指標由2.89近似階躍地增加至 15.09并越過當前時刻閾值，說明電網當前有擾動事件發生，這可能與所設節點 21和節點95有功負荷突變有關。因此，依據步驟5），需對引起該事件的擾動節點進行快速定位。

圖3 Max-ESCM指標計算結果Fig.3 The results from Max-ESCM

由步驟6）與步驟7），計算Min-ESCM及其對應閾值γv。然后，根據步驟8）與步驟9），對當前時刻Min-ESCM指標中的每個元素進行判斷，輸出所有大于其動態閾值的元素，即為引起該擾動事件的節點編號，如圖4所示。

圖4 Min-ESCM指標計算結果Fig.4 The results from Min-ESCM

對比圖4a與圖4b可知，在t1000采樣時刻，因電網中沒有擾動事件發生，所以Min-ESCM指標中各元素分布相對隨機均勻，且均未超過對應閾值。但在t1001采樣時刻，因電網中有擾動事件發生，使得系統隨機性被破壞，此時Min-ESCM指標中的第21號元素與第 95號元素的數值相較于其他元素出現了顯著變化，并且超過了當前閾值，故可判定引起該事件的擾動節點編號為 21與 95，這與實際設置的擾動情況一致。因此，這說明了本文所提方法對電網同時多重擾動定位是有效的。

4.2 相繼多重擾動定位測試

對于相繼多重擾動場景的模擬，分別設置了“滑動窗內”場景，即節點117與相鄰節點14在間隔一個滑動窗內相繼出現有功負荷躍變，和“滑動窗外”場景，即節點117與節點14在間隔一個滑動窗外相繼出現有功負荷躍變，見表2。

表2 一類相繼躍變的合成負荷Tab.2 The synthetic loads with successive step-up change

按照第 3節擾動識別的步驟，可分別得到上述兩種場景下Max-ESCM指標的變化曲線，如圖5所示。由圖5a與圖 5b可知，在擾動時刻Max-ESCM指標的值均超過了當前閾值，說明系統此時有擾動事件發生，需進一步對擾動節點進行定位。此外，在再次出現擾動事件的時刻，Max-ESCM曲線均又一次出現了波動，且波動點發生時刻與相繼擾動發生時刻幾乎一致。因此，可判斷此時存在相繼擾動事件。

圖5 Max-ESCM指標計算結果Fig.5 The results from Max-ESCM

根據第3節擾動節點定位的步驟，可分別得到上述兩種場景下Min-ESCM指標的變化曲線，如圖6所示。從圖6a與圖6b可知，在t1001采樣時刻，兩種場景下的Min-ESCM指標中均僅有第117號元素大于當前時刻動態閾值，這表明該方法能夠對電網首先發生的單一擾動事件進行有效定位。

圖6 Min-ESCM指標計算結果Fig.6 The results from Min-ESCM

其次，對比圖6a與圖6b可得，當相繼多重擾動發生時刻大于或小于時間窗時，Min-ESCM指標在同一時刻將出現不同現象。對于“滑動窗內”場景，由于 14號節點有功負荷在相隔 200采樣時刻之后發生了躍變，其相繼擾動發生時刻小于一個滑動窗寬度，此時 Min-ESCM指標中的14號元素與117號元素在t1200采樣時刻均超過了當前閾值，故可判斷出 14號節點此時有擾動事件發生。而對于“滑動窗外”場景，因為14號節點有功負荷相隔300采樣時刻之后才發生了躍變，其相繼擾動發生時刻大于一個滑動窗寬度，故在t1300采樣時刻的 Min-ESCM指標中僅有14號元素出現突變，從而定位出14號節點此時有擾動事件發生。可見該方法能夠有效地定位電網相繼多重擾動。

4.3 多重短路故障定位測試

設置節點21與相隔較遠的節點95在采樣時刻t1001發生了三相短路，故障持續時間為0.49s。采樣時刻為t1050后，系統恢復正常，共2 000個采樣時刻。根據第3節擾動識別的步驟，可得該場景下Max-ESCM指標的變化曲線，如圖7所示。

圖7 Max-ESCM指標計算結果Fig.7 The results from Max-ESCM

由圖7可見，在t1001采樣時刻Max-ESCM指標達到了22.67且越過了當前時刻閾值。依據第3節步驟 5），可判斷出系統當前時刻有擾動事件發生，需對引起該擾動事件的節點進行快速定位。由第 3節擾動定位的流程可得圖8。

圖8 Min-ESCM指標計算結果Fig.8 The results from Min-ESCM

觀察圖8a與圖8b可得，Min-ESCM指標中的第21號元素與第95號元素在t1001采樣時刻的數值明顯與其他元素不同（其他元素均小于0.1），而且超過了當前閾值。這與實際設置的擾動節點情況一致，可見所提方法對多重短路故障事件定位是有效的。

5 討論

為支撐本文所提方法的潛在應用，對 PMU信號的時間延遲問題，以及 PMU數據的時效性和空間位置的關聯性對電網多重擾動識別和定位的事前輔助可行性討論如下。

1）關于PMU信號的時間延遲問題。由于PMU發送數據延時抖動、通信協議、通信鏈路、傳送距離、通信鏈路的負載情況、通信通道帶寬等因素的影響，上傳到調度控制中心的 PMU實測數據的延時變化范圍較大，具有一定的分布特性，一般在幾十毫秒到幾百毫秒之間變化[25]。針對該問題，一方面，從時間維度來看，可考慮數字信號處理方法來處理，例如：由于不同PMU數據的等效采樣頻率可能不同，可以認為等效采樣頻率低的數據類型在采樣間隔內數值相等[26]；另一方面，從空間維度來看，可考慮文獻[11]所提分區并行計算策略，一個省級或地市級電網范圍內設置一個分區子站，一個大區或省級電網范圍內設置一個主站，以減少通信距離，降低時間延遲范圍。

2）PMU數據的時效性和空間位置的關聯性對電網多重擾動識別和定位的事前輔助是可行的。一方面，基于實時數據的閾值設定是電網多重擾動識別的一個關鍵因素。傳統上，基于RMT的擾動識別與定位方法的研究中，其判據一般均采用靜態閾值，即往往僅考慮理論模型數據或工程經驗。文獻[9,27]綜合 PMU歷史數據和實時數據的動態閾值或許更為合理有效。因此，可以認為PMU的時效性有助于動態閾值的快速生成[11]，繼而有助于電網態勢感知有效性的提高[28]。另一方面，基于空間位置的分區預警是電網多重擾動識別和定位的一個重要環節。文獻[29]指出，不論是單一擾動還是復雜擾動都可能包含傳導型關系。復雜擾動在發展過程中同時跨越多個行政區域的情況較少，體現出區域電網間較強的抗擾動能力。這與文獻[11]一致。因此，可以認為利用 PMU空間位置的關聯性，能夠有助于對此類擾動提前預警。此外，需要注意的是，相對傳統的離線靜態分區預警，結合 PMU數據的在線動態分區預警也是一個值得關注的領域。

6 結論

本文基于Max-ESCM及Min-ESCM提出了一種適用于多重擾動場景的電網擾動定位方法。借助DIgSILENT和Matlab軟件，通過IEEE 54機118節點系統算例，表明了該方法的有效性，同時得到以下結論：

1）相較傳統基于 RMT的電網擾動定位方法，本文所提方法無需進行多次運算，能夠在判定電網發生擾動事件的同時，輸出引起該擾動的所有節點。

2）基于實對稱矩陣特征向量的“相變”現象，利用 Max-ESCM指標及Min-ESCM指標實現了對電網不同場景下多重擾動的快速定位，有助于電網多重擾動的排查。