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

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

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

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

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

1 理論基礎(chǔ)

1.1 Periodic PO-MOESP算法

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

對于線性時間-周期系統(tǒng)，其狀態(tài)空間方程描述如下：

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

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

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

LTP系統(tǒng)的可觀測性矩陣被定義為

接下來假設(shè)系統(tǒng)對于所有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。由于系統(tǒng)的周期性，可以證明和其中指周期的數(shù)量。根據(jù)這些定義我們可以得到：

因此輸出矩陣可以被記為

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

同時下式也成立：

進行RQ分解：

通過變換可以證明：

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

利用可觀測矩陣的結(jié)構(gòu)，取其前l(fā)行作為系統(tǒng)矩陣。矩陣Uk的i行到j(luò)行表示為Uk(i:j,:)，因此輸出矩陣表示為

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

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

若采用離散時間狀態(tài)空間模型，相應(yīng)各階模態(tài)的固有頻率、阻尼比及振型為

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

1.2 算法實現(xiàn)步驟

(1)對結(jié)構(gòu)施加激勵并采集響應(yīng)信號；

(2)利用位移響應(yīng)信號構(gòu)造Hankel矩陣；

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

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

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

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

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

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

2 仿真算例

2.1 風(fēng)機塔架模型

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

圖2 風(fēng)機整體結(jié)構(gòu)示意圖

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

通過有限元仿真軟件對風(fēng)機塔架有限元模型進行模態(tài)分析，計算得到結(jié)構(gòu)模型前3 階的模態(tài)參數(shù)作為參考值，如表1所示。

表1 塔架前3階振動頻率

為了驗證Periodic PO-MOESP 算法的識別效果，現(xiàn)分別采用正弦信號激勵和多頻率諧波信號激勵進行基于瞬態(tài)響應(yīng)信號的結(jié)構(gòu)動態(tài)特性分析。

2.2 基于正弦激勵下的模態(tài)參數(shù)識別

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

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

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

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

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

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

圖4 前3階結(jié)構(gòu)振型對比圖

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

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

2.3 基于諧波激勵下的模態(tài)參數(shù)識別

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

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

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

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

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

圖7 前3階結(jié)構(gòu)振型對比圖

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

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

2.4 抗噪性能

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

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

圖9 測點7響應(yīng)數(shù)據(jù)加入5%白噪聲后前后對比圖

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

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

圖10 前3階結(jié)構(gòu)振型對比圖

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

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

3 試驗驗證

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

圖12 試驗方案原理圖

圖13 試驗件和裝置

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

表5 頻率識別結(jié)果與試驗結(jié)果對比

表6 阻尼識別結(jié)果與試驗結(jié)果對比

由識別結(jié)果可知，采用本文算法識別梁結(jié)構(gòu)模態(tài)參數(shù)精度較高，前3 階頻率識別結(jié)果與試驗值的誤差范圍均在10%以內(nèi)，1 階頻率誤差為2.05%；2階頻率的誤差為3.93%；3階頻率的誤差為-7.16%。前3 階模態(tài)振型與試驗?zāi)B(tài)振型的對比如圖14 所示。由識別與試驗結(jié)果對比可以看出，此算法在周期激振的條件下對頻率、阻尼和振型的識別都很準(zhǔn)確。

圖14 前3階結(jié)構(gòu)振型對比

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

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

4 結(jié)語

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