一種快速矩陣束相量提取方法的研究

楊 洋， 肖湘寧， 陳鵬偉， 房 釗

(新能源電力系統(tǒng)國家重點實驗室， 華北電力大學， 北京 102206)

一種快速矩陣束相量提取方法的研究

楊 洋， 肖湘寧， 陳鵬偉， 房 釗

(新能源電力系統(tǒng)國家重點實驗室， 華北電力大學， 北京 102206)

電力系統(tǒng)的監(jiān)控、保護及模態(tài)識別通常依賴于相量參數(shù)的準確程度，因而快速準確地提取各類相量十分重要。本文基于矩陣束算法，利用前后兩個耦合時段的輸入信號采樣陣構成一個方陣，該方陣的特征值不僅表示前后兩個時刻工頻、諧波及其他各類分量的復幅值之比，也表示著各分量的極值點信息。通過求解該方陣的特征根，在已知初始相量前提下，即可得到所需各分量的相量值。為減少求解計算量，本文利用方陣的Krylov子空間來構造相量的最小多項式，快速求解代表特定分量的特征值。理論分析和仿真結果表明，本文提出的快速矩陣束相量提取方法，與傳統(tǒng)矩陣束法相比，計算量顯著減小且實現(xiàn)方式簡單，不受非目標分量及頻率偏移的影響，在較短數(shù)據(jù)窗內仍具有較高精度。該方法適用于工頻信息的提取，也適用于各次諧波信息。

相量； 矩陣束； Krylov子空間； 最小多項式

1 引言

電力系統(tǒng)工頻以及其他各類諧波相量參數(shù)，在電力系統(tǒng)事件錄波、狀態(tài)估計、系統(tǒng)保護、系統(tǒng)控制及系統(tǒng)仿真等各個方面有著廣泛的應用。

隨著對相量或其中某個特征量測量的深入研究，國內外學者已提出了很多相量測量算法，如三點算法[1]、離散傅里葉變換(Discrete Fourier Transform，DFT)算法及其改進算法[2-5]、最小二乘法(Least Square，LS)曲線擬合、Prony算法、矩陣束法[6,7](Matrix Pencil，MP)、卡爾曼濾波算法[8,9]、基于泰勒展開模型的算法[10]、小波變換[11]、定頻率相量提取算法[12]等。

三點算法原理簡單，運算量小，但其測量精度低，易受諧波和噪聲干擾。廣泛運用的DFT算法由于未考慮到非整次諧波和衰減非周期分量的影響，會產(chǎn)生頻譜泄露和欄柵效應，其改進算法雖然在一定程度上改善了性能，但都帶來了迭代計算、數(shù)據(jù)窗限制等問題。LS算法是將輸入的暫態(tài)電氣量與一個預設的含有非周期分量及某些諧波分量的函數(shù)模型按最小平方誤差原理進行擬合[13]，其雖然可以消除任意暫態(tài)分量，但提取精度依賴于預設信號模型的完整程度，而預設模型的復雜性也會增加計算時間?？柭鼮V波算法與LS算法一樣，需要預設模型，根據(jù)設定的狀態(tài)數(shù)目進行卡爾曼濾波，濾除各次諧波與噪聲。

Prony算法與矩陣束算法都是基于衰減指數(shù)和的信號模型，前者通過對特征多項式的求解來確定信號模型的極點信息[14]，后者能通過奇異值分解來確定模型階數(shù)，并構造兩個相關矩陣求解與極點信息等價的廣義特征值[6,7]。但是，上述兩個基于衰減指數(shù)和模型的方法，其計算精度依賴于較長的數(shù)據(jù)窗與較高的模型階數(shù)，而且矩陣求逆、奇異值分解或其他矩陣運算過程都較為復雜，不可避免地增加了計算量。小波變換、基于泰勒展開、定頻率相量提取算法等方法雖然在穩(wěn)定性方面得到改善，并且有效地提高了計算結果的精度，但同時也不可避免地增加了計算量、依賴較長的數(shù)據(jù)窗、引入迭代或存在其他缺陷。

上述方法都無法將精度、響應時間和運算速度有效統(tǒng)一起來，即現(xiàn)有單一方法無法完全滿足電力系統(tǒng)實際工況下相量提取的實時性與準確性要求。

本文在矩陣束算法的基礎上，提出了一種適用于電力系統(tǒng)電磁暫態(tài)仿真的快速矩陣束相量提取(Fast Matrix Pencil，F(xiàn)MP)算法，利用采樣時段耦合的兩個輸入信號采樣陣構成一個方陣，從該方陣特征值中獲得所需相量變化信息。為減少計算量，利用方陣的Krylov子空間來構造向量的最小多項式，實現(xiàn)對所需部分特征值的快速求解，同時得到極值點信息與相量變化信息，實現(xiàn)對工頻與其他各類諧波的相量提取。

2 算法原理

2.1 矩陣束法

電力系統(tǒng)中大部分被觀測信號都可以看作是一系列具有任意振幅、相位、頻率和衰減因子的指數(shù)函數(shù)的和，即由衰減余弦分量和衰減直流分量組成：

(1)

式中，信號y(t)分別有q1個衰減余弦分量和q2個衰減直流分量；Ai、Aj為對應分量幅值；θi為對應分量相位；αi,αj為衰減因子；fi為對應分量頻率。

令M=2q1+q2，則式(1)的離散時間函數(shù)形式為：

(2)

(3)

(4)

式中，Δt為采樣間隔；zi表示信號模型的極點信息；pi表示留數(shù)信息。

由采樣序列y(n)構造出2個(N-L)×L的Hankel矩陣Y1與Y2：

(5)

(6)

式中，L稱為束參數(shù)，合理地選取可以在一定程度上消除數(shù)據(jù)所含噪聲，通常取值范圍在N/3～N/2之間，從而Y1與Y2可表示為：

Y1=Z1PZ2

(7)

Y2=Z1PZ0Z2

(8)

式中

Z0=diag[z1,z2,…,zM]

P=diag[p1,p2,…,pM]

定義矩陣束

(9)

容易證明，Y2-λY1的秩為M，而當λ=zi時，Z0-λI的第i行等于零，矩陣的秩降低為M-1，可以看出，信號的極點即為矩陣對{Y2,Y1}的廣義特征值，故而對上述極點信息的求取轉化為求解式(10)特征值的過程，其表達式為：

(10)

求得zi與M后，留數(shù)pi可利用最小二乘法的線性方程式(11)計算得到：

(11)

根據(jù)式(12)～式(15)，可求解出信號中包括基波分量在內的各類分量的幅值、相位、角頻率以及衰減因子：

(12)

(13)

(14)

(15)

2.2 快速矩陣束法

本文所提出的快速矩陣束法是利用采樣時段耦合的兩個輸入信號采樣來構成一個方陣，并從該方陣特征值中獲得所需相量變化信息，實現(xiàn)相量的快速提取。

對于y(kΔt)，分別從前一時刻t1和后一時刻t2開始截取數(shù)據(jù)窗長為N的序列，按式(16)與式(17)構成(N-L+1)×LHankel采樣矩陣YE與YL：

(16)

(17)

式中，Δn表示兩個測量起始時刻之間相差的采樣次數(shù)，針對穩(wěn)態(tài)與暫態(tài)情況可取不小于1的不同值。

式(16)與式(17)可寫成：

YE=Z1PZ2

(18)

YL=Z1P′Z2

(19)

將YE的偽逆左乘YL，得到：

(20)

由式(20)可知：

(21)

YL重新表示為：

YL=Z1PZΔnZ2

(22)

定義一個新的矩陣束如下：

(23)

(24)

在相鄰時刻，極值點變化量可忽略，即近似認為極值點信息相同。因此在初始穩(wěn)定狀態(tài)、已知初始相量的情況下，采用上述方法時，除在第一次求解時需要構造兩個采樣矩陣，后續(xù)相量提取可利用前一時段的采樣陣。并且該方法無需求解模態(tài)數(shù)目M，即可以同時提取出較為準確的工頻、各類諧波分量前后時刻相量變化信息。

2.3 特征值快速求取

對信號中某個或若干分量的提取工作并不需要完全求解所有特征值，本文通過構造向量的最小多項式來快速求解矩陣的目標特征值，以減小計算量，提高快速矩陣束方法的計算速度。

2.3.1 最小多項式的構造

已有相關定理證明[15]，對于任意矩陣A，存在非零m次的最小多項式fm(x)，使得fm(A)r0=0，且fm(x)=0的每個根一定都是矩陣的特征值，其中r0為A的某一列向量。因此矩陣A的部分特征值可通過求解fm(x)來獲得。為求得向量r0的最小多項式，需找到最小的正整數(shù)m，滿足r0,Ar0,…，Amr0線性相關，然后按照式(25)求解得到r0的最小多項式，如式(26)所示：

[r0,Ar0,…,Amr0](x0,x1,…,xm-1,1)T=0

(25)

f(z)=x0+x1z+…+xm-1zm-1+zm=0

(26)

由于構造的矩陣YE與YL須包含待觀測信號所有的模態(tài)信息，矩陣階數(shù)可能較大，會導致向量組r0,Ar0,…,Amr0的線性相關性開始變得模糊，并且Krylov矩陣[r0,Ar0,…,Amr0]不可避免變得病態(tài)，即不容易確定向量r0的最小多項式次數(shù)和精確求解線性方程組式(25)。

通過歸一化Krylov[r0,Ar0,…,Amr0]矩陣的列，降低其病態(tài)程度，然后對其做QR分解，以R的秩作為r0的最小多項式次數(shù)，并求解：

(27)

從而得到r0的最小多項式：

(28)

2.3.2 最小多項式的求解

由于僅提取所需分量的對應特征值，且特征值所代表的極點信息變化范圍一般較小，所以根據(jù)式(24)，給定較為合適的迭代初始值，迭代求解出目標特征值，從而可在較小的迭代次數(shù)下獲得高精度的零點值。本文采用密勒法求解多項式零點。

3 仿真驗證

本文在Matlab中實現(xiàn)快速矩陣束算法，通過給定信號評估各類諧波、衰減分量、不同數(shù)據(jù)窗長、頻率偏移對工頻及諧波相量提取精度的影響，并與離散傅里葉算法、最小二乘算法、Prony算法及傳統(tǒng)矩陣束算法的提取結果進行比較。本文中，若無特殊說明，所有仿真的采樣頻率取10kHz，基波頻率為50Hz，數(shù)據(jù)窗長取一個周波20ms，測量誤差為0～40ms內連續(xù)測量的最大誤差。

3.1 各種諧波和衰減分量的影響

為便于比較不同分量對快速矩陣束相量提取方法的影響，在不同給定信號中重復給定特定分量，構造的待測信號分別為：

(1)含整次諧波信號

10cos(3ω0t+π/4)

(2)含整次、非整次諧波信號

10cos(3.35ω0t+π/4)

(3)含整次、非整次諧波及衰減直流分量信號

10cos(3.35ω0t+π/4)+10e-t/0.1

(4)含整次、非整次諧波及衰減余弦分量信號

10cos(3.35ω0t+π/4)+

10e-t/0.1cos(4.45ω0t+π/3)

分別以上述信號中的基波相量，3.35次諧波相量為提取目標，連續(xù)滑窗測量40ms，采用本文的快速矩陣束法對信號x1(t)、x2(t)提取結果分別如圖1和圖2所示。

圖1 含整次諧波信號滑窗測量工頻相量結果Fig.1 Fundamental phasor extraction results containing integer harmonic components

圖2 含非整次諧波信號滑窗測量3.35次諧波相量結果Fig.2 3.35 harmonic phasor extraction results of signal containing non-integer harmonic components

針對含直流衰減分量和余弦衰減分量的信號，分別以2次諧波與3.35次諧波分量為目標提取分量，滑窗測量幅值最大誤差分別為3.01e-9%和6.53e-10%，圖形與圖1和圖2類似，限于篇幅，不再重復列出。仿真結果表明，快速矩陣束相量提取結果不受非目標分量的影響。

分別采用離散傅里葉算法、最小二乘算法、Prony算法及傳統(tǒng)矩陣束算法對上述4種信號進行測量，以基波分量為提取目標，其中離散傅里葉算法為全周傅里葉算法，最小二乘法的預設模型為不大于5次的工頻、整次諧波與衰減分量的組合。提取結果如表1所示，其中Eamp表示幅值提取最大相對誤差，Eang表示相位提取最大絕對誤差。

表1 不同信號類型下各算法的提取性能對比

表1的結果表明快速矩陣束方法由于累計誤差的存在，在相角提取誤差上稍大于傳統(tǒng)矩陣束算法，但幅值提取僅通過極值點信息獲得，因此結果更為精確。其總體精度遠優(yōu)于傅氏算法及最小二乘算法。

3.2 頻率偏移的影響

當實際頻率偏離額定值時，各分量相量的幅值、相角的測量誤差都與初始相角有關。以0°為初始相角，分別采用DFT算法、最小二乘法、Prony算法及快速矩陣束算法對信號x3(t)連續(xù)滑窗測量，設定偏移頻率為±0.5Hz與±1Hz，即取實際頻率為49Hz、49.5Hz、50.5Hz、51Hz，以20ms數(shù)據(jù)窗連續(xù)測量40ms，采樣頻率設定為10kHz。幅值提取結果如圖3所示。

圖3 頻率偏移下各算法幅值提取性能對比Fig.3 Comparison of amplitude extraction performance among various algorithms with frequency offset

圖3的仿真結果表明，快速矩陣束算法同Prony算法一樣，在幅值提取上幾乎不受頻率偏移的影響，且遠優(yōu)于DFT算法及最小二乘算法。

3.3 數(shù)據(jù)窗長的影響

對含不同分量的待測信號而言，其最小窗長限制也不同。針對信號x3(t)，改變窗長進行工頻相量的滑窗測量，測量結果如表2所示，其中×表示幅值相對誤差大于30%，相位絕對誤差大于30°。

由表2可以看出，數(shù)據(jù)窗長為基波周期非整數(shù)倍時，測量結果依然精確；幅值測量誤差隨著數(shù)據(jù)窗長的增大而減小，算法抗噪性能逐漸增強；相角測量誤差由于累計作用，存在一定波動性，但始終保持較為準確的提取精度。

4 結論

本文提出了一種適用于電力系統(tǒng)電磁暫態(tài)仿真的快速矩陣束相量提取方法，它具有如下特點。

表2 不同數(shù)據(jù)窗長對快速矩陣束算法的影響

(1)與傳統(tǒng)矩陣束算法相比，幅值提取僅涉及到矩陣求逆運算，因此提取結果更為準確。

(2)采用本文構造的向量最小多項式來求解部分特征值，實現(xiàn)簡單，極大地減小了計算量。

(3)工頻、各類諧波相量的提取結果不受非目標分量及頻率偏移的影響。

(4)對數(shù)據(jù)窗長的限制小，提取結果在較短數(shù)據(jù)窗長內仍具有較高的精度，且數(shù)據(jù)窗越長，幅值提取結果越精確。

[ 1] 王暉，常鮮戎，鄭煥坤(Wang Hui, Chang Xianrong, Zheng Huankun)．基于改進瞬時對稱分量與三點算法的相量測量新算法(A method based on improved instantaneous symmetrical components and three-point algorithm for synchronized phasor measurement)[J]．電力系統(tǒng)保護與控制(Power System Protection and Control)，2011，39(19)：115-120．

[ 2] Yang J Z，Liu C W．A precise calculation of power system frequency and phasor[J]．IEEE Transactions on Power Delivery, 2000，15(2)：494-498．

[ 3] Maohai Wang，Yuanzhang Sun．A practical, precise method for frequency tracking and phasor estimation[J]．IEEE Transactions on Power Delivery，2004，19 (4)：1547-1552．

[ 4] T S Sidhu，Xudong Zhang，V Balamourougan．A new half-cycle phasor estimation algorithm[J]. IEEE Transactions on Power Delivery，2005，20(2)：1299-1305．

[ 5] David Hart，Damir Novosel，Yi Hu，et al．A new frequency tracking and phasor estimation algorithm for generator protection[J]．IEEE Transactions on Power Delivery，1997，12(3)：1064-1073．

[ 6] Sarkar T K, Pereira O. Using the matrix pencil method to estimate the parameters of a sum of complex exponentials[J]. IEEE Antennas and Propagation Magazine, 1995,37(1): 48-55.

[ 7] Hua Y, Sarkar T K. Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[J]. IEEE Transactions on Acoustics, Speech and Signal Processing, 1990,38(5): 814-824.

[ 8] Kennedy K, Lightbody G，Yacamini R．Power system harmonic analysis using the Kalman filter[A]. IEEE Power Engineering Society General Meeting [C]. 2003. 2：757-761．

[ 9] Routray A，Pradhan A K，Rao K P．A novel Kalman filter for frequency estimation of distorted signals in power systems[J]．IEEE Transactions on Instrumentation and Measurement,2002，51(3)：469-479．

[10] 麥瑞坤,何正友，薄志謙(Mai Ruikun, He Zhengyou, Bo Zhiqian). 基于泰勒展開模型的同步相量估計新算法(A new transient method for earth fault line selection in resonant grounded power system based on morphological filtering)[J].電力系統(tǒng)自動化(Automation of Electric Power Systems), 2008,32(12):22-26.

[11] Wilkinson W A，Cox M D．Discrete wavelet analysis of power system transients[J]．IEEE Transaction on Power Systems，1996，11(4)：2038-2044．

[12] 索南加樂，王斌，王莉, 等(Suonan Jiale, Wang Bin, Wang Li, et al.)．快速相量提取算法(A fast phasor calculation algorithm for power systems)[J]．中國電機工程學報(Proceedings of the CSEE)，2013，33(1)：123-129．

[13] Sachdev M S，Giray MM．A least error square technique for determining power system frequency [J]．IEEE Transactions on Power Apparatus and Systems，1985，PAS-104(2)：437-444．

[14] 丁屹峰，程浩忠，呂干云, 等(Ding Yifeng, Chen Haozhong, Lv Ganyun, et al.). 基于Prony算法的諧波和間諧波頻譜估計 (Spectrum estimation of harmonics and interharmonics based on Prony algorithm)[J].電工技術學報(Transactions of China Electrotechnical Society), 2005，20(10):94-98.

[15] 李宏(Li Hong)．矩陣特征值問題的若干數(shù)值算法研究(Research on some numerical algorithms for solving eigenvalue problem)[D]．廣州：華南理工大學(Guangzhou: South China University of Technology)，2011．

Research on fast matrix pencil method for phasor extraction

YANG Yang, XIAO Xiang-ning, CHEN Peng-wei, FANG Zhao

(State Key Laboratory for Alternate Electrical Power System with Renewable Energy Sources, North China Electric Power University, Beijing 102206, China)

Power system monitoring, protection and modal identification rely on the accuracy of phasor parameters, thus quickly and accurately extracting phasor is very important. This paper proposes a fast matrix pencil (FMP) method to estimate parameters of electrical signals based on matrix pencil (MP) method. A square matrix was formed by two time domain coupling matrices of input sampling signals. Its eigenvalues represent the changing information of fundamental and harmonic phasors, as well as the pole information of all components. Under the knowledge of initial phasors, the concerned phasors can be obtained by solving the square matrix eigenvalues. In order to reduce computation burden, this paper uses the Krylov subspace of square matrix to construct a phasor minimal polynomial, which can be used to quickly calculate concerned eigenvalues of specific components. Theoretical analysis and simulation results demonstrate that the FMP method significantly reduces the calculation burden and is simple to conduct as compared with the MP method. It is immune to the effect of non-target components and frequency offset, also could keep a high accuracy with short data window. The method can be applied to the extraction of fundamental phasors, as well as harmonic phasors.

phasor; matrix pencil method; Krylov subspace; minimal polynomial

2015-04-20

中央高?；究蒲袠I(yè)務費專項資金資助項目(2015XS22)

楊 洋(1989-)， 男， 河北籍， 博士研究生， 研究方向為電力系統(tǒng)仿真與分析； 肖湘寧(1953-)， 男， 湖南籍， 教授， 博士生導師， 從事電力系統(tǒng)電能質量等方面的教學與科研工作。

TM711

A

1003-3076(2016)02-0001-06