周期激勵下的風機塔架模態識別方法

胡嘉苗，杭曉晨，朱 銳，曹芝腑，姜 東,

(1.南京林業大學機械電子工程學院，南京210037；2.東南大學空天機械動力學研究所, 南京211189)

準確獲取風機塔架的模態參數對于及早發現風機異常、避免風機事故具有至關重要的作用[1]。通過結構模態辨識可得到結構的固有頻率和振型，其被廣泛應用于振動系統分析、振動故障診斷與預測、結構動力優化設計、結構損傷檢測等方面[2–3]。由于實際環境中模擬風力機的載荷困難，且由于塔筒細長、鋼壁較薄等，難以制作縮比試驗模型，以及載荷特殊輸入量難測量等問題，運行模態分析(Operational modal analysis，OMA)作為一種基于實測數據的分析工具，在研究風力機動力學特性方面具有廣闊的應用前景[4–5]。結構動力學領域最早應用于風機模態參數測試的OMA算法是James和Carne的自然激勵技術(Natural excitation technique，NExT)[6–7]，該技術不依賴于人工激勵，大大加快了數據采集速度。目前，環境激勵下模態參數識別方法主要有：峰值拾取法[8]、時間序列法[9]、隨機減量法[10]、NExT 法、隨機子空間法[11]等。

由于周期性的重力、氣動載荷作用以及轉子旋轉等因素，運行中的風機的穩定性分析應該在一個周期性的框架內進行[12–13]。在過去的十幾年間，子空間模態辨識算法得到了飛速發展，并且在多輸入-多數出（Multiple-input multiple-output，MIMO)系統辨識領域取得極大的成功[14–15]。然而現有的時域方法都是考慮在平穩的白噪聲激勵下的情況，對于激勵中含有周期成分的情況研究較少[16]。為了識別運行中的風機塔架的模態參數，選擇一種能夠識別具有大量測量通道的周期系統的算法是很重要的。Allen 等[17–18]確定了在運行中的風機的周期性振型，但精度較低。Yang等[19]識別出了風力渦輪機葉片的高分辨率模態振型，但這些振型是在轉子不工作時的靜態振型。Verhargen 等最早提出了MOESP 算法[20]，該算法直接從輸入輸出數據中估計出系統廣義能觀空間的基底從而完成模態參數辨識[21]。Verhargen[22]在1994年提出了PO-MOESP 算法(PO指past output)，是MOESP 算法的一種變形，其思想是通過處理系統的輸出矩陣方程，經過投影計算得出可觀測矩陣的列空間，從而獲得模態參數。Verdult 等[23]通過減小可觀測矩陣的維數，減小了計算系統矩陣的工作量。然而，對于受到周期激勵作用的結構，很難有效地辨識其高分辨率振型。鄧先來等[22]利用一種無相移數字濾波器對響應信號進行處理，消除了周期激勵對模態辨識的影響。董霄峰等[24]通過改進特征系統實現法(Eigensystem realization algorithm，ERA)與概率密度函數法(Probability density function，PDF)結合的運行模態識別方法，剔除了不同工況下葉輪和機組轉頻的諧波成分干擾。

對運行中的結構進行運行模態分析，需要準確測量在環境激勵作用下的各個測點的響應[25]。本文研究一種在環境激勵下基于周期性-過去輸出多變量輸出誤差狀態空間算法(Periodic PO-MOESP)的僅輸出可測模態參數識別方法。相比其他算法，該方法可以同時處理過程噪聲與測量噪聲，在周期激勵的作用下，可以準確識別出頻率、阻尼和振型。通過有限元仿真得到結構在周期性激勵作用下的位移響應信號，通過周期性采樣構造Hankel 矩陣，基于Periodic PO-MOESP 算法對響應信號進行識別，得到可觀矩陣，進而推導出系統矩陣并求解系統模態參數。采用一總長為60 m 的圓筒模型作為仿真算例，研究該算法的辨識精度和抗噪能力。最后將此算法應用于懸臂梁點模態試驗，驗證了該算法在周期激勵下進行結構模態參數識別的有效性。

1 理論基礎

1.1 Periodic PO-MOESP算法

運行中的風力渦輪機與停放的風力渦輪機在性質上是不同的，如同直升機和渦輪機等許多重要的系統，必須用線性時間-周期（LTP）運動方程來建模，以正確辨識風機模態參數[26]。

對于線性時間-周期系統，其狀態空間方程描述如下：

式中：x(k)為系統的狀態變量矩陣，x(k)∈Rn；y(k)為系統的測量數據矩陣，由經數學處理后的觀測信號組成，y(k)∈Rm；w(k)為輸入的白噪聲（過程噪聲矩陣），均值為零；ν(k)為輸出白噪聲（測量噪聲矩陣），均值為零；A(k)為系統狀態傳遞矩陣；C(k)為輸出傳遞矩陣。

矩陣A(k)、F(k)和C(k)都是離散周期為k的周期函數，為了把這個系統寫成一種通常用于子空間的形式，定義了輸出數組：

同樣，過程噪聲和測量噪聲也是這樣的數組。參數d為窗口大小，且P≥d＞n。

LTP系統的可觀測性矩陣被定義為

接下來假設系統對于所有k，除以長度d，都是可觀測的，因此有滿列秩，即rank(=n。定義：

=C(k+i-1)Φ(k+i-1,k+j)F(k+j-1)，i=2,…，d；j=1，…，d-1并且i大于j。由于系統的周期性，可以證明和其中指周期的數量。根據這些定義我們可以得到：

因此輸出矩陣可以被記為

Periodic PO-MOESP 目的是利用數據矩陣來識別可觀測矩陣?？勺C：

同時下式也成立：

進行RQ分解：

通過變換可以證明：

觀測矩陣的值域可以通過計算的奇異值分解得到：

利用可觀測矩陣的結構，取其前l行作為系統矩陣。矩陣Uk的i行到j行表示為Uk(i:j,:)，因此輸出矩陣表示為

矩陣A(k)通過求解下列方程得到：

根據模態分析理論可知，求解矩陣A的特征值問題即可得到模態的固有頻率和阻尼比[27]。當采用連續時間狀態空間模型時，已知A的特征值λi和特征向量hi，固有頻率、阻尼比及振型可由式(17)得到：

若采用離散時間狀態空間模型，相應各階模態的固有頻率、阻尼比及振型為

其中：fs為采樣頻率[28]。

1.2 算法實現步驟

(1)對結構施加激勵并采集響應信號；

(2)利用位移響應信號構造Hankel矩陣；

(3)對Hankel矩陣進行QR分解得到R矩陣；

(4)對R矩陣中的R21按式(13)進行奇異值分解；

(5)通過擴展可觀測性矩陣Odk求解A和C；

(6)求A的特征值得到模態固有頻率fi、阻尼比εi；通過矩陣C求解模態振型。

圖1為基于Periodic PO-MOESP算法的流程圖。

圖1 周期力激勵方法流程圖

2 仿真算例

2.1 風機塔架模型

典型的三葉片風力機主要由旋轉葉片、機箱、塔架、基座等部分組成，齒輪箱、傳動機構、發電機和配電裝置都安裝在機箱中。到目前為止，風力機葉片長度最大可達100 米。葉片旋轉時受到風載作用，使得整體結構所承受載荷呈現明顯的周期性。塔架作為整個結構的支撐部件，振動對其的影響最為明顯。因此，風機倒塌往往是由塔架的斷裂引起的。風機塔筒模型結構及測點示意圖如圖2所示。

圖2 風機整體結構示意圖

本算例中將風力機整體結構進行簡化，主要分析塔架的動態特性。以長度為60 m 的圓筒結構作為研究對象，外徑為1.5 m，內徑為1.4 m，沿長度方向等距設置7 個測點，測量相對根部的側向位移。梁結構的有限元模型材料的彈性模量為2×105MPa，密度為7 850 kg/m3，泊松比為0.3，將測點1 固定作為邊界條件。在實際風機模態辨識實驗中，響應信號可通過加速度、速度和位移傳感器進行測量。

通過有限元仿真軟件對風機塔架有限元模型進行模態分析，計算得到結構模型前3 階的模態參數作為參考值，如表1所示。

表1 塔架前3階振動頻率

為了驗證Periodic PO-MOESP 算法的識別效果，現分別采用正弦信號激勵和多頻率諧波信號激勵進行基于瞬態響應信號的結構動態特性分析。

2.2 基于正弦激勵下的模態參數識別

本節以上述建立的風機塔架有限元模型為基礎，給有限元模型的測點7 施加F=10 sin(40π×t)N的正弦激勵力，激勵力隨時間的變化曲線如圖3所示。圖中所示為前5 s的激勵信號，激勵信號施加時間共為10 s。通過Periodic PO-MOESP 方法對結構7 個測點的位移響應信號進行模態參數識別，將識別的模態參數與有限元模型的理論計算值進行對比。

圖3 正弦激勵力隨時間變化曲線

施加簡諧激勵力后，識別結果和仿真結果對比如表2所示。

表2 正弦激勵下的識別結果與計算結果對比

分析可知，在測點7施加正弦激勵后，可識別出結構的前3 階頻率和阻尼比。由對比結果可以得出，算法辨識結果與模態仿真結果的誤差較小。

前3階仿真模態振型與辨識模態振型的對比圖如圖4所示。

圖4 前3階結構振型對比圖

由振型圖對比結果可知，各階模態振型的理論結果和識別結果比較接近，模態振型連續、準確，與精細有限元模態預分析結果具有相似的模態振型特征。將利用算法得到的結果與根據有限元仿真得到的振型結果的MAC 值對比，從圖5 中可以看出，識別出的振型與有限元仿真得到的模態振型匹配度較高，模態置信度能達到0.98以上，表明該算法可以有效識別周期激勵下的風機模態參數。

圖5 識別振型與仿真振型MAC值對比

2.3 基于諧波激勵下的模態參數識別

在實際運行過程中，當結構處于高轉速狀態時風力機會受到明顯的葉輪倍頻和諧波激勵的作用，導致結構振動中的諧波成分占比較大，風機結構的模態頻率與諧波成分較為接近，嚴重影響算法的識別精度。通過在不同時間段施加包含不同頻率諧波成分的信號，模擬風機在實際運行過程中在不同風速下轉子、葉片等結構旋轉引起的周期激勵。諧波激勵力隨時間變化曲線如圖6所示。激勵信號施加時間共為10 s，諧波成分為10 Hz、15 Hz、20 Hz的周期信號。

將圖6的諧波輸入信號作用在測點7上，保持采樣頻率、采樣間隔與正弦激勵時一致，采樣頻率為100 Hz，采樣時間為10 s。通過瞬態響應分析得到7個測點的位移響應信號，數據總長為1 000。利用Periodic PO-MOESP 算法對輸出的數據進行模態辨識，結構模態頻率和阻尼特性的仿真結果和識別結果如表3所示。

表3 諧波激勵下的識別結果與計算結果對比

圖6 諧波激勵力隨時間變化曲線

分析結果可知，在諧波激勵的條件下，算法識別結果與有限元仿真結果吻合，通過辨識得到的前3階頻率和阻尼比的大小與仿真結果誤差較小。前3階仿真模態振型與辨識模態振型的對比如圖7所示。

圖7 前3階結構振型對比圖

由理論與識別結果對比可以看出，此算法在諧波激勵力輸入的條件下對模態的識別很準確。通過MAC值對比圖8可以看出，識別振型與理論振型匹配度較高，模態置信度達到0.98以上。

圖8 諧波激勵條件下辨識振型與仿真振型MAC值對比

2.4 抗噪性能

由于實際工程中不可避免存在噪聲激勵，本節研究5 %白噪聲激勵對模態辨識精度的影響以及Periodic PO-MOESP算法的抗噪性能。

將2.2節中正弦激勵后得到的測點1至測點7的位移響應信號加入5%的白噪聲，其中測點7響應數據加入白噪聲后的信號前后對比如圖9所示。將經過處理后的數據利用Periodic PO-MOESP算法進行模態辨識，得到如表4所示的結果。

圖9 測點7響應數據加入5%白噪聲后前后對比圖

表4 加入5%白噪聲激勵后的識別結果與計算結果對比

分析結果可知，在諧波激勵的條件下，算法識別結果與有限元仿真結果吻合，通過辨識得到的前3階頻率和阻尼比的大小與仿真結果誤差較小。前3階仿真模態振型與辨識模態振型的對比如圖10所示。

圖10 前3階結構振型對比圖

將加入5%白噪聲后利用算法所識別振型與根據有限元仿真得到的振型結果的MAC值對比，見圖11，可知識別振型與仿真振型匹配度較高，模態置信度能達到0.98以上。說明加入5%白噪聲后算法模態參數辨識精度較高，抗噪性能較好。

圖11 識別振型與仿真振型MAC值對比

3 試驗驗證

以一端固支的柔性梁為試驗件，利用激振器在梁頂部沿厚度方向施加周期激勵。梁總長108 cm，截面寬度為3 cm，梁厚度為0.3 cm，采用和風機塔架相同的材料合金鋼，其彈性模量為206 Gpa，泊松比為0.3，密度為7 850 kg/m3。從梁的底端開始，每隔20 cm布置一個電荷式單向加速度傳感器，用以測量加速度響應信號。試驗方案如圖12 所示。試驗件和裝置如圖13所示。

圖12 試驗方案原理圖

圖13 試驗件和裝置

梁結構底端通過夾具固定，激振器采用磁吸座固定。對試驗結構施加40 Hz 的周期加速度信號，大小為1 g。采集加速度傳感器的響應數據，采樣頻率為100 Hz，數據總長為20 000。通過東華DH5922D動態信號測試分析系統識別加速度信號，得到前3 階頻率試驗值為1.369 Hz、7.860 Hz、21.284 Hz。利用本文算法對響應信號進行識別，得到的頻率和阻尼識別結果與試驗值的比較如表5和表6所示。

表5 頻率識別結果與試驗結果對比

表6 阻尼識別結果與試驗結果對比

由識別結果可知，采用本文算法識別梁結構模態參數精度較高，前3 階頻率識別結果與試驗值的誤差范圍均在10%以內，1 階頻率誤差為2.05%；2階頻率的誤差為3.93%；3階頻率的誤差為-7.16%。前3 階模態振型與試驗模態振型的對比如圖14 所示。由識別與試驗結果對比可以看出，此算法在周期激振的條件下對頻率、阻尼和振型的識別都很準確。

圖14 前3階結構振型對比

將基于本文算法識別得到的結果與采用動態信號測試分析系統得到的試驗振型結果進行MAC 值對比，從圖15 中可以看出，算法所識別振型與試驗振型匹配度較高，模態置信度能達到0.96以上。

圖15 試驗振型與識別振型MAC值對比

4 結語

針對風力機運行模態試驗中難以測得激勵信號的問題，提出了一種僅需要位移輸出信號的結構模態辨識方法。仿真算例中采用一60 m長圓筒結構，基于Periodic PO-MOESP算法對其進行模態參數識別，并將此識別結果和有限元正則模態分析結果進行對比。對比正弦激勵和諧波激勵下的模態參數辨識結果，結果顯示不同周期激勵的形式對該算法結構模態參數辨識的影響較小。在響應信號中加入白噪聲，辨識結果證明該方法具有良好的抗噪性能。將此算法應用于一懸臂梁模態試驗，驗證了本文算法在周期激勵條件下識別結構模態參數的有效性。