利用矩陣束算法提取水力機(jī)組參數(shù)的振蕩特性

曾保安，曾 云；張立翔（昆明理工大學(xué)建筑工程學(xué)院工程力學(xué)系，昆明 650500）

前言

水力發(fā)電機(jī)組的穩(wěn)定性是衡量機(jī)組性能的指標(biāo)之一，機(jī)組的穩(wěn)定性將關(guān)系到整個(gè)電站及其水工廠房的安全性、國民經(jīng)濟(jì)利益和整個(gè)電力系 統(tǒng)能否安全運(yùn)行。因此，機(jī)組穩(wěn)定性受到學(xué)者們廣泛關(guān)注，而表征穩(wěn)定性的參數(shù)分別是振動(dòng)、擺度和壓力脈動(dòng)，其中振動(dòng)是機(jī)組穩(wěn)定運(yùn)行最重要的指標(biāo)。

水力發(fā)電機(jī)組的振動(dòng)按其形式可分為水力振動(dòng)、機(jī)械振動(dòng)和電磁振動(dòng)。這三種振動(dòng)耦合作用于水力機(jī)組，以致造成水力機(jī)組振動(dòng)機(jī)理很復(fù)雜。雖然國內(nèi)外很多研究者對其做出了巨大的努力，但對其振動(dòng)機(jī)理的理論研究仍然欠缺。一些研究者試圖從振動(dòng)信號中提取機(jī)組的振蕩信息，找出誘發(fā)機(jī)組振蕩的誘因，以便達(dá)到改善機(jī)組穩(wěn)定性的效果。目前對振動(dòng)信號的提取方法有傅里葉變換法、短時(shí)傅里葉變換法、HHT變換法、小波分析法、prony算法、矩陣束算法等。其中傅里葉變換法不能反應(yīng)振蕩的阻尼特性和瞬時(shí)頻率；HHT變換法會出現(xiàn)模態(tài)混疊的現(xiàn)象；小波分析方法是目前信號分析比較完備的方法，但是存在波基選取困難的局限；文獻(xiàn)[1]、[2]用prony方法對電力系統(tǒng)暫態(tài)或者大規(guī)模系統(tǒng)干擾進(jìn)行穩(wěn)定性研究，發(fā)現(xiàn)該方法對噪音很敏感，抗噪能力很差。文獻(xiàn)[3]、[4]提出一種簡化分析的方法，不依賴于對象模型可應(yīng)用于從現(xiàn)場實(shí)測數(shù)據(jù)提取機(jī)組參數(shù)振蕩特性，可對機(jī)組任一參數(shù)的振蕩特性及其影響因素進(jìn)行詳細(xì)的分析。文獻(xiàn)[5]、[6]研究的是矩陣束算法，此種方法最關(guān)鍵的步驟在于極點(diǎn)求解時(shí)進(jìn)行了去噪處理，采用奇異值分解(SVD)和降低秩的方法把數(shù)據(jù)中隱含的噪音過濾掉，以免噪音產(chǎn)生虛假極點(diǎn)，其抗噪能力有很大提高。文獻(xiàn)[8～12]將幾種不同形式的矩陣束算法做了計(jì)算精度和計(jì)算復(fù)雜程度的比較性研究。驗(yàn)證了矩陣束算法具有極強(qiáng)的抗噪能力和廣泛使用性。

基于以上的研究結(jié)果，本文將矩陣束算法做了一些改進(jìn)，并用改進(jìn)的矩陣束算法提取同一工況下水力發(fā)電機(jī)組角速度、功角、有功、進(jìn)口處水頭、流量和水輪機(jī)出力等主要參數(shù)的模態(tài)信息，進(jìn)而建立運(yùn)動(dòng)模型，初始工況、擾動(dòng)強(qiáng)度和擾動(dòng)方向是參數(shù)振蕩特性的主要影響因素，本文將對這些主要因素如何定量影響振幅、頻率和衰減因子的變化進(jìn)行全面分析，總結(jié)出模態(tài)信息之間存在的耦聯(lián)關(guān)系，這些研究和探討為控制和設(shè)計(jì)提供依據(jù)。

1 矩陣束算法

1.1 矩陣束基本算法

（1）模態(tài)數(shù)的確定。通過采樣數(shù)據(jù)y(1)，y(2),…，y(N)構(gòu)造一個(gè)(N-L+1)×L階 Hankel矩陣：

通過奇異值分解，Y=UΣVT，其中σi是其對角矩陣Σ對角線上的第i個(gè)奇異值。由于現(xiàn)場實(shí)測數(shù)據(jù)和模型存在噪聲，會產(chǎn)生虛假極點(diǎn)。可用歸一化奇異值(σi/σmax)≥β（β為閾值）來辨識出真實(shí)極點(diǎn)與虛假極點(diǎn)。把發(fā)生突變的歸一化奇異值作為初始閾值，逐漸增大或者減小閾值，選取信噪比達(dá)到最大時(shí)所對應(yīng)的最大下標(biāo)i作為最大模態(tài)數(shù)M，那么前M個(gè)極點(diǎn)為真實(shí)極點(diǎn)，其對應(yīng)的歸一化奇異值作為閾值。

（2）構(gòu)造2個(gè)Hankel矩陣。由左奇異矩陣V的前M個(gè)主奇異向量構(gòu)成濾噪矩陣V′=[v1,v2,…,vM]，構(gòu)造出兩個(gè)Hankel矩陣：

其中，V1、V2分別是V′刪掉第一行、最后一行而得到；Σ′的前M個(gè)奇異值與Σ相同，其余的均為0；

（3）極點(diǎn)與留數(shù)的求解。通過對矩陣束Y1和Y2定義，得到以下關(guān)系：

由（3）可知，當(dāng) λ≠zi時(shí)，對角陣Z0-λI的第i行不為零，即Y2-λY1的秩為M，當(dāng) λ=zi時(shí)，Z0-λI的第i行等于0；矩陣Y2-λY1的秩降為M-1，由相關(guān)理論得出，極點(diǎn)zi為矩陣束{Y2,Y1}的廣義特征值；即：

其中，Y1+是Y1的偽逆矩陣。

留數(shù)的求解可通過最小二乘法計(jì)算：

極點(diǎn)與留數(shù)求解完畢，利用(10)式求解振幅Ai、相位θi、頻率fi和衰減因子αi。

逼近函數(shù)為：

（4）擬合精度的衡量。擬合曲線與原始數(shù)據(jù)曲線之間的重合度用信噪比(SNR)來衡量；真實(shí)數(shù)據(jù)為y(i)，逼近數(shù)據(jù)為x(i)。

一般認(rèn)為SNR大于20就是一種可以接受的逼近結(jié)果，該值越大，逼近效果越好。

1.2 算法改進(jìn)

（1）重構(gòu)數(shù)據(jù)。式(9)右邊的數(shù)據(jù)y(1)，y(2)，…，y(n)應(yīng)使用消噪后的模型數(shù)據(jù)y′(1)，y′(2)，…，y′(n)，通過式子Y′=UΣ′VT求出矩陣Y′，把Y′反對角線上的數(shù)據(jù)取平均值作為模型數(shù)據(jù)。

（2）目前矩陣束算法所研究的振蕩模型通常是某個(gè)量圍繞平衡點(diǎn)的振蕩，其相應(yīng)的相位角只分布在一、四象限，可用（10）式直接求出。但本文研究的振蕩模型屬于階躍式的振蕩，其相位角應(yīng)根據(jù)留數(shù)的實(shí)部與虛部對應(yīng)分布到四個(gè)象限上。

（3）影響矩陣束算法精度的因素

1）數(shù)據(jù)長度至少應(yīng)取到包含全部振蕩信息。

2）采樣頻率保證大于2fmax，fmax是低頻率振蕩的最大頻率。

3）閾值是消除噪音合理與否的關(guān)鍵；一般選取突變的歸一化奇異值作為初始閾值，逐漸增大或者減小，信噪比最大時(shí)所對應(yīng)的歸一化奇異值作為閾值。

4）Hankel矩陣列數(shù)L也會影響矩陣束算法的精度，一般取L/3-L/2。

2 水力機(jī)組模型

2.1 水力系統(tǒng)模型

水力系統(tǒng)采用一管多機(jī)的形式，如圖1所示，機(jī)組2和機(jī)組3假設(shè)為穩(wěn)定運(yùn)行狀態(tài)。鋼管與支管可用等效管i等效。

圖1 水力系統(tǒng)模型

其中，x=[x1x2x3x4x5]，x4=q為流量，x5=y為導(dǎo)葉開度。

水力系統(tǒng)數(shù)學(xué)模型：

fp(i)為第i路支管與鋼管的等效管道損失系數(shù)，Z(i)為第i路支管與鋼管的等效管道涌浪阻抗，T(i)為第i路支管與鋼管的等效管道彈性時(shí)間，y為主接力器位移，fpT為隧道水頭損失系數(shù)，TwT為隧道水流慣性時(shí)間系數(shù)，Ty為主接力器時(shí)間常數(shù)，qi為i路支管機(jī)組的流量，qj為除第i路支管外機(jī)組的流量。

2.2 水輪機(jī)力矩計(jì)算模型

水輪機(jī)力矩模型：

式中，mt為水輪機(jī)力矩，At為水輪機(jī)增益系數(shù)，qnl為水輪機(jī)空載流量，ht為水輪機(jī)進(jìn)口處的水頭，q為水輪機(jī)流量。

2.3 發(fā)電機(jī)模型

發(fā)電機(jī)采用單機(jī)無窮大三階實(shí)用模型：

mt為輸入機(jī)械力矩；Ef為勵(lì)磁電動(dòng)勢；Td′0為d軸開路暫態(tài)時(shí)間常數(shù)；Tj為機(jī)組慣性時(shí)間常數(shù)ωB=314rad/s為角速度基值；ω為機(jī)組角速度；Δω=ω－1；δ為功角；Us為無窮大系統(tǒng)電壓；XdΣ=Xd+XT+XL；Xd為d軸電抗；XT為變壓器電抗；XL為線路等效電抗；XqΣ=Xq+XT+XL；Xq為q軸電抗；D為阻尼系數(shù)；E′q為q軸瞬變電動(dòng)勢;X′dΣ=X′d+XT+XL；X′d為d軸次暫態(tài)電抗。

3 模型計(jì)算

水力系統(tǒng)模型、發(fā)電機(jī)模型、水輪機(jī)出力模型、典型的并聯(lián)PID調(diào)速器和PI勵(lì)磁控制器構(gòu)成完整的水力機(jī)組系統(tǒng)進(jìn)行模型計(jì)算。

3.1 發(fā)電機(jī)功角計(jì)算模型

初始工況p=0.5p.u.,目標(biāo)工況pc=1p.u.，在正常調(diào)節(jié)時(shí)，其建模的計(jì)算過程如下：

（1）按頻率10Hz進(jìn)行采樣，得到功角的一組仿真數(shù)據(jù)，并利用這些數(shù)據(jù)構(gòu)成一個(gè)Hankel矩陣。通過奇異值分解，選取突變的歸一化奇異值0.000048作為初始閾值，逐漸增大或減小閾值，以信噪比達(dá)到最大時(shí)所對應(yīng)的下標(biāo)i作為最大模態(tài)數(shù)M=7，對應(yīng)的閾值β= 0 .0006。

（2）構(gòu)造兩個(gè)Hankel矩陣和求解模型數(shù)據(jù)。

（3）求出極點(diǎn)與留數(shù)。求出的極點(diǎn)有3個(gè)實(shí)數(shù)和2對共軛復(fù)數(shù)，通過(10)式可知，實(shí)數(shù)極點(diǎn)對應(yīng)的頻率為0；復(fù)數(shù)極點(diǎn)以共軛成對形式出現(xiàn)，由特征值分析理論可知，每一對復(fù)數(shù)極點(diǎn)對應(yīng)一個(gè)振蕩模態(tài)。

（4）通過(11)式把振蕩頻率為0與不為0的運(yùn)動(dòng)模型分別表示出來。

功角振蕩頻率為0的運(yùn)動(dòng)模型：

功角振蕩頻率不為0的運(yùn)動(dòng)模型：

功角運(yùn)動(dòng)模型為δ1(t+δ2)（t）即

采用上式計(jì)算并與仿真數(shù)據(jù)進(jìn)行比較，結(jié)果如圖2所示，功角運(yùn)動(dòng)模型曲線δ(t)與仿真曲線基本上完全重合。

圖2 p=0.5p.u., pc=1p.u.時(shí)功角的響應(yīng)曲線

3.2 水輪機(jī)出力計(jì)算模型

采用與功角相同的工況和計(jì)算步驟，提取水輪機(jī)出力的運(yùn)動(dòng)模型。

水輪機(jī)出力振蕩頻率為0的運(yùn)動(dòng)模型：

水輪機(jī)出力振蕩頻率不為0的運(yùn)動(dòng)模型：

水輪機(jī)出力的運(yùn)動(dòng)模型為：

采用上式計(jì)算并與仿真數(shù)據(jù)進(jìn)行比較，結(jié)果如圖3所示。水輪機(jī)出力運(yùn)動(dòng)模型曲線mt(t)與仿真曲線基本上完全重合。

圖3 p=0.5p.u., pc=1p.u.時(shí)水輪機(jī)出力的響應(yīng)曲線

采用相同的方法研究機(jī)組主要參數(shù)的運(yùn)動(dòng)模型時(shí)，發(fā)現(xiàn)模態(tài)信息之間存在耦聯(lián)關(guān)系，即：

（1）進(jìn)口處的水頭、流量和水輪機(jī)出力三個(gè)水力參數(shù)的振蕩頻率相近；功角、有功和角速度三個(gè)電氣參數(shù)的振蕩頻率相近；衰減因子α2相近；各參數(shù)的衰減因子r1相近。

（2）水力參數(shù)的擾動(dòng)頻率為 0.48Hz左右；電氣參數(shù)的振動(dòng)頻率為：第一基頻0.51Hz左右，第二基頻1.1Hz左右；水力參數(shù)的頻率相近，電氣參數(shù)的頻率相近，說明參數(shù)的振蕩模態(tài)只與參數(shù)類型有關(guān)。因此，在調(diào)節(jié)過程中，參數(shù)產(chǎn)生振蕩是機(jī)組結(jié)構(gòu)本身所固有的，其頻率為固有頻率，其振蕩模態(tài)為固有振蕩模態(tài)。

通過采用相同的方法提取各參數(shù)的模態(tài)信息，并建立運(yùn)動(dòng)模型，且與文獻(xiàn)[1]建立的運(yùn)動(dòng)模型基本一致。因此，可將水力機(jī)組參數(shù)運(yùn)動(dòng)模型歸結(jié)為兩種基本形態(tài)：

1）周期衰減運(yùn)動(dòng)模型

其中，Ai、αi、fi、θi分別為第i個(gè)振蕩模態(tài)的振幅，衰減因子、頻率、相位角。m是周期衰減運(yùn)動(dòng)階數(shù)。

2）過阻尼運(yùn)動(dòng)模型

xz為運(yùn)動(dòng)終值（系統(tǒng)變量的平衡值），Ci、ri分別是第i個(gè)過阻尼運(yùn)動(dòng)模態(tài)的振幅、衰減因子，n是過阻尼運(yùn)動(dòng)階數(shù)。

4 模型關(guān)聯(lián)仿真

4.1 初始工況對各參數(shù)模態(tài)信息的影響

初始工況p對功角、水輪機(jī)出力模態(tài)信息的影響分別如表1、表2，其中，目標(biāo)工況pc=1p.u.不變。

表1 初始工況對功角模態(tài)信息的影響

表2 初始工況對水輪機(jī)出力模態(tài)信息的影響

采用相同的方法提取其他各電氣參數(shù)、水力參數(shù)的模態(tài)信息隨初始工況變化分別與表1、表2中的變化規(guī)律一致。即：

（1）各參數(shù)的過阻尼運(yùn)動(dòng)的階數(shù)、周期衰減運(yùn)動(dòng)的階數(shù)不隨初始工況而變化，各參數(shù)的信噪比隨著初始工況減小而減小。

（2）各主要參數(shù)的振蕩頻率受初始工況影響小；電氣參數(shù)之間的頻率是相近的，水力參數(shù)之間頻率也相近。

（3）各參數(shù)的各階振幅隨初始工況減小而增大。

（4）各參數(shù)的衰減因子r1和α1隨初始工況減小而減小；各參數(shù)的衰減因子r2隨初始工況減小而增大；電氣參數(shù)的衰減因子α2隨初始工況減小而增大。

4.2 擾動(dòng)強(qiáng)度對各參數(shù)模態(tài)信息的影響

初始工況p=0.7p.u.，擾動(dòng)強(qiáng)度Δp對功角、水輪機(jī)出力模態(tài)信息的影響分別如表3、表4。

表3 擾動(dòng)強(qiáng)度對功角模態(tài)信息的影響

表4 擾動(dòng)強(qiáng)度對水輪機(jī)出力模態(tài)信息的影響

采用相同的方法提取其他各電氣參數(shù)、水力參數(shù)的模態(tài)信息隨擾動(dòng)強(qiáng)度的變化分別與表3、表4中的變化規(guī)律一致。即：

（1）各參數(shù)的過阻尼運(yùn)動(dòng)的階數(shù)、周期衰減運(yùn)動(dòng)的階數(shù)不隨擾動(dòng)強(qiáng)度而變化，各參數(shù)的信噪比隨著擾動(dòng)強(qiáng)度增大而減小。

（2）各參數(shù)的振蕩頻率受擾動(dòng)強(qiáng)度影響小，電氣參數(shù)的振蕩頻率相近，水力參數(shù)的振蕩頻率也相近，且電氣參數(shù)的振蕩頻率比水力參數(shù)的振蕩頻率多一階。

（3）各參數(shù)的振幅隨擾動(dòng)強(qiáng)度增大而增大，不受擾動(dòng)方向的影響。

（4）各參數(shù)的衰減因子r1相近；電氣參數(shù)的衰減因子α2相近；各參數(shù)的各階衰減因子隨正向擾動(dòng)強(qiáng)度增大而減小；各階衰減因子隨負(fù)向擾動(dòng)強(qiáng)度增大而增大。

（5）改變初始工況，Δp與表3、表4一致，可以得出各參數(shù)的模態(tài)信息隨擾動(dòng)強(qiáng)度的變化不受初始工況影響。

5 結(jié)論

（1）改進(jìn)的MP算法用于提取水力機(jī)組各參數(shù)模態(tài)信息有效性得以驗(yàn)證，體現(xiàn)出良好的抗噪能力；該方法可通過實(shí)測數(shù)據(jù)提取各參數(shù)模態(tài)信息，不依賴于對象系統(tǒng)模型，因此，可用于電網(wǎng)和機(jī)組的在線監(jiān)測和故障診斷。

（2）通過模態(tài)信息建立的兩種基本運(yùn)動(dòng)形態(tài)，能直觀、有效地對水力機(jī)組的振蕩特性進(jìn)行描述；可通過振動(dòng)頻率推斷出誘發(fā)振動(dòng)的誘因。

（3）探討了各種影響因素下各參數(shù)模態(tài)信息的變化趨勢和耦聯(lián)關(guān)系，為控制設(shè)計(jì)提供依據(jù)。

[1]Hauer J F, Demeure C J, Scharf L L. Initial results in prony analysis of power system response signals [J].IEEE Trans on Power Systems, 1990, 5(1): 80-89.

[2]Grund C E, Paserba J J, J F Hauer. Comparison of prony and eigenanalysis for power system control design[J]. IEEE Transactions on Power Systems,1993,8 (3): 964-971.

[3]曾 云, 張立翔, 王煜. 水輪發(fā)電機(jī)組振蕩特性的簡化分析方法[J]. 電機(jī)與控制學(xué)報(bào), 2009, 13(增刊1): 25-29.

[4]曾云, 沈祖詒, 曹林寧. 低頻振蕩下發(fā)電機(jī)響應(yīng)特性的量化分析[J]. 電力系統(tǒng)及其自動(dòng)化學(xué)報(bào),2008, 20(6): 83-87.

[5]Tapan K.Sarkar , odilon Pereira. 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 .

[6]Tapan Kumar Sarkar, Jinhwan Koh. Application of the matrix pencil method for estimating the SEM(Singularity Expansion Method) poles of source-free transient responses from multiple look directions[J]. IEEE Transactions on Antennas and Propagation, 2000, 48(4): 612-618.

[7]Kunder P. 電力系統(tǒng)穩(wěn)定與控制[M]. 北京:中國電力出版社, 2002.

[8]Muhammad Faisal Khan, Muhammad Tufail.Comparative analysis of various matrix pencil methods for direction of arrival estimation[C].Image Analysis and Signal Processing (IASP),International Conference on,Islamabad,Pakistan,9-11 April 2010.

[9]朱瑞可, 李興源. 矩陣束算法在同步電機(jī)參數(shù)辨識中的應(yīng)用. 電力系統(tǒng)及其自動(dòng)化學(xué)報(bào), 2012, 36(6): 52-56.

[10]Fran?ois Sarrazin, Ala Sharaiha. Comparison between matrix pencil and prony methods applied on noisy antenna responses[C]. Loug hborough Antennas & Propagation Conference, Loughborough, UK, 14-15 November 2011.

[11]Yanhui Liu, Zaiping Nie. Reducing the number of elements in a linear antenna array by the matrix pencil method[J]. IEEE Transactions on Antennas and Propagation,2008, 56(9):2955-2962.

[12]Nuri Yilmazer, Jinhwan Koh. Utilization of a unitary transform for efficient computation in the matrix pencil method to find the direction of arrival[J]. IEEE Transactions on Antennas and Propagation, 2006, 54(1): 175-181.