基于狀態(tài)空間模型的大型風(fēng)力機運行模態(tài)及不確定性分析

張驚朝戴靠山，2，3施袁鋒，2，*

（1.四川大學(xué)土木工程系，成都 610065；2.深地科學(xué)與工程教育部重點實驗室，成都 610065；3.破壞力學(xué)與防災(zāi)減災(zāi)四川省重點實驗室，成都 610065）

0 引言

風(fēng)能作為一種遍布全球的可再生能源，備受全世界的重視，風(fēng)電已逐步成為最具發(fā)展前景的清潔能源之一［1］。據(jù)全球風(fēng)能理事會（Global Wind Energy Council，GWEC）統(tǒng)計，截至2019年底，全球累計風(fēng)電裝機容量約為650.6 GW，而2014年底約為369.8 GW，5年時間內(nèi)將近增加了一倍［2］。為了更有效捕捉風(fēng)能并提高風(fēng)電經(jīng)濟效益，風(fēng)力機組向著葉片更長和塔架更高的趨勢發(fā)展。然而，大型風(fēng)力機組長期在環(huán)境惡劣和荷載復(fù)雜條件下運行，各種風(fēng)電場的運維問題也隨之而來，尤其是風(fēng)力機組因振動引起的損傷甚至倒塔帶來的影響最為顯著［3-6］。為延長風(fēng)力機組的使用壽命以及減少運維成本，定期對風(fēng)力機組進行結(jié)構(gòu)健康監(jiān)測是重要維護手段之一［7］，其中通過對風(fēng)機結(jié)構(gòu)進行模態(tài)分析獲得其動力特性，是風(fēng)機結(jié)構(gòu)損傷識別和性能評估的主要手段之一［8］。

模態(tài)分析可分為試驗?zāi)B(tài)分析（Experimental Modal Analysis，EMA）和 運 行 模 態(tài) 分 析（Operational Modal Analysis，OMA）［9］。EMA通常用激振器施加并測量人工激勵以及獲取相應(yīng)的結(jié)構(gòu)振動響應(yīng)來開展模態(tài)分析。與傳統(tǒng)的EMA相比，OMA無須對結(jié)構(gòu)施加人工激勵，僅需直接拾取結(jié)構(gòu)在運行或環(huán)境激勵下的振動響應(yīng)來估計結(jié)構(gòu)的模態(tài)參數(shù)，因此可避免人工激勵對結(jié)構(gòu)可能造成的損傷問題。Wang等［10］分別運用接觸式和非接觸式測量，對風(fēng)力機展開OMA；孫溢膺等［11］采用隨機子空間（SSI）法，應(yīng)用于某風(fēng)機模型實測中。盡管學(xué)者們在提出或者采用某種模態(tài)分析方法時，都進行了一些比較來說明識別結(jié)果的準確性，但大部分并沒有真正考慮識別結(jié)果的不確定性問題。在模態(tài)分析中，由于數(shù)據(jù)有限、測量噪聲、建模誤差和激勵未知等各種不確定因素的存在，一般不可能精確地識別出結(jié)構(gòu)的模態(tài)參數(shù)［12］。實現(xiàn)模態(tài)參數(shù)不確定性量化，一方面，可以直接看出每個參數(shù)識別的精度，有利于對模態(tài)分析方法的評估；另一方面，有利于模型的修正，從而減少識別誤差，得到更精準的識別結(jié)果［13］。Mares等［14］采用多元回歸模型計算不確定參數(shù)的靈敏度矩陣，并利用梯度方法計算得到了模型參數(shù)的均值和協(xié)方差；姜東等［15］以三自由度體系和復(fù)合材料板為例，考慮識別結(jié)果的不確定性，提出了一種基于區(qū)間分析的不確定性結(jié)構(gòu)動力模型修正方法。

目前，模態(tài)分析方法可以有效地識別風(fēng)力機組停機工況下的模態(tài)參數(shù)［16-17］，但對于帶有葉片旋轉(zhuǎn)的風(fēng)機運行狀態(tài)下的識別結(jié)果往往存在不精準甚至失效的問題［18-19］。這主要是由于風(fēng)機在運行狀態(tài)下，葉片旋轉(zhuǎn)會產(chǎn)生較大的氣動阻尼、振動信號中帶有明顯的周期性成分以及識別算法中時不變系統(tǒng)和穩(wěn)態(tài)隨機激勵的假設(shè)矛盾，影響了識別結(jié)果的準確性［16，20］。為此，董霄峰［21］利用一種考慮諧波修正的SSI方法，考慮了葉片轉(zhuǎn)動對模態(tài)識別的影響；Dai等［22］提出了一種改進的SSI法，討論了該方法在識別精度和穩(wěn)定上的優(yōu)越性以及對諧波頻率、噪聲程度與諧波幅度等干擾的魯棒性，進而對某在役風(fēng)機停機和運行工況下的模態(tài)參數(shù)進行識別并通過共振校核方法評估風(fēng)機結(jié)構(gòu)具有良好的運行安全性；孟歡等［23］利用快速傅里葉變換以及譜細化法，應(yīng)用于兩座5 MW海上風(fēng)力機的現(xiàn)場實測數(shù)據(jù)中。

本文采用一種基于狀態(tài)空間模型的運行模態(tài)及不確定性分析方法，以解決帶有葉片旋轉(zhuǎn)的風(fēng)機結(jié)構(gòu)的模態(tài)參數(shù)識別問題。最后，本文以實例分析一風(fēng)機不同塔架方向的模態(tài)參數(shù)的差異性，以及不同工況下風(fēng)機模態(tài)參數(shù)的變化情況，并結(jié)合識別結(jié)果的不確定性水平，證明該方法的實用性。

1 基于狀態(tài)空間模型和極大似然估計的運行模態(tài)分析

1.1 隨機狀態(tài)空間模型

在環(huán)境激勵下，長度為N的結(jié)構(gòu)響應(yīng)數(shù)據(jù)Y={y1，y2，…，yN}，可用離散時間隨機狀態(tài)空間模型進行模擬

式中：A∈Rn×n和C∈Rp×n分別為系統(tǒng)的狀態(tài)矩陣和輸出矩陣；xk∈Rn×1為系統(tǒng)在k時刻的n維狀態(tài)向量，并假設(shè)系統(tǒng)初始狀態(tài)向量x1服從正態(tài)分布x1～N(μ1，P1)；yk∈Rp×1為系統(tǒng)在k時刻的p維輸出向量；wk、vk分別為系統(tǒng)的過程噪聲和測量噪聲，可模擬為均值為零的高斯白噪聲，且兩者的協(xié)方差表示為

式中：E(·)為數(shù)學(xué)期望因子；δpq為Kronecker delta函數(shù)。

式（2）考慮了系統(tǒng)過程噪聲和測量噪聲的相關(guān)性。這里，同時假設(shè)x1與(wk，vk)不相關(guān)。

隨機狀態(tài)空間模型式（1）的參數(shù)化形式具體可表示成

式中，θ為模型式（1）的參數(shù)變量，代表了系統(tǒng)所包含的所有自由變化參數(shù)組成的向量形式，系統(tǒng)的特性也完全由θ決定。

1.2 模型參數(shù)極大似然估計

系統(tǒng)參數(shù)θ可通過觀測響應(yīng)數(shù)據(jù)Y進行估計。由于極大似然估計方法具有漸近無偏估計、漸近一致性和漸近有效性等優(yōu)異特性，本文采用極大似然法來實現(xiàn)隨機狀態(tài)空間模型參數(shù)θ的估計，進而獲取結(jié)構(gòu)的運行模態(tài)參數(shù)。

極大似然估計的關(guān)鍵在于似然函數(shù)的建立以及求解，但由于xk、wk以及vk未知，直接通過隨機狀態(tài)空間模型式（1）來構(gòu)建似然函數(shù)難以實現(xiàn)，而結(jié)合卡爾曼濾波可以實現(xiàn)。式（1）的卡爾曼濾波可以表示［24］

式中，

由于卡爾曼濾波式（4）得到的殘差序列{ek(θ)}是不相關(guān)的高斯隨機序列，由此可構(gòu)造出數(shù)據(jù)Y的對數(shù)似然函數(shù)L(θ|Y)，去掉常數(shù)項后得

根據(jù)極大似然法的原理可知，系統(tǒng)參數(shù)θ的極大似然估計可表示為

由于直接通過極大似然估計求解對數(shù)似然函數(shù)L(θ|Y)難以完成，本文應(yīng)用EM迭代算法求解參數(shù)θ的極大似然估計值。

1.3 最大期望(EM)算法

EM算法是一種通過迭代進行極大似然估計的優(yōu)化算法，主要用于包含隱含量的概率模型進行參數(shù)估計，即從觀測數(shù)據(jù)Y中找到使對數(shù)似然函數(shù)L(θ|Y)最大化的估計參數(shù)θ［25］。假設(shè)系統(tǒng)的完整數(shù)據(jù)Z=(X，Y)由X，Y共同構(gòu)成，其中，X為遺漏數(shù)據(jù)（未知），Y為觀測數(shù)據(jù)。

由貝葉斯公式可知：

兩邊同時取對數(shù)可得：

由此，EM迭代算法主要分為E步和M步，以此循環(huán)迭代，直至滿足設(shè)定的收斂條件。①E步：計算條件期望Q(θ)；②M步：求函數(shù)Q(θ|)的最大值，得到新的估計值。迭代收斂條件為

式中，ρ為目標收斂率。

具體應(yīng)用EM算法求解隨機狀態(tài)模型式（1）的參數(shù)θ的極大似然估計值中，我們可將模型的狀態(tài)向量{x1，x2，…，xN+1}視為遺漏數(shù)據(jù)X，離散響應(yīng)數(shù)據(jù)Y={y1，y2，…，yN}為已知數(shù)據(jù)，EM每一步的估計值為［26］

其中，

式（14）中所有條件期望可通過卡爾曼濾波和卡爾曼平滑求得［26］。

2 模態(tài)參數(shù)的不確定性分析

2.1 克拉美羅下界(CRB)

在統(tǒng)計估計理論中，CRB對任何無偏估計量的函數(shù)的協(xié)方差確定了一個下界，其公式可以表示為［27］

式中，Γ(θ)為所考慮模型參數(shù)的函數(shù)，如模型參數(shù)、模態(tài)參數(shù)等；為Γ(θ)對θ的雅可比矩陣，?表示矩陣的偽逆；I(θ)表示系統(tǒng)參數(shù)化模型的Fisher信息矩陣（FIM），其具體元素可表示為

式中，i，j表示FIM的元素位置。

針對隨機狀態(tài)空間模型的似然函數(shù)式（6），有［28］

在計算參數(shù)的CRB時，其關(guān)鍵在于FIM的計算。但由于FIM的維度一般較大，導(dǎo)致運行時間過長，為了在保證計算精度的情況下縮減運行時間，本文采用了由Shi等［29］提出的快速計算FIM的方法。

2.2 模態(tài)形式狀態(tài)空間方程

由于結(jié)構(gòu)頻率和阻尼比只與A的特征值有關(guān)，將一般形式的隨機狀態(tài)空間方程轉(zhuǎn)換為模態(tài)形式后計算FIM，可直接獲得模態(tài)參數(shù)的協(xié)方差信息。對隨機狀態(tài)空間方程式（1）做相似變換，令xk=，可轉(zhuǎn)換為模態(tài)形式的狀態(tài)空間方程：

式（19）中每個復(fù)模態(tài)分塊矩陣與結(jié)構(gòu)頻率和阻尼的對應(yīng)關(guān)系為

式中：fi，ζi為第i個復(fù)模態(tài)的頻率和阻尼比；Ts為采樣步長；j2=-1。

另外，式（18）中模態(tài)形式的Cˉ矩陣為一特殊歸一化的形式，具體表達式為

其中，[1 0]或1代表了每一模態(tài)振型最大幅值所對應(yīng)的位置。

相似變換矩陣T可由式（19）和式（21）的模態(tài)形式Aˉ和Cˉ反推求得，具體參考文獻［29］。

此時，模態(tài)形式的狀態(tài)空間模型式（18）的參數(shù)變量可寫成

式中，λ=[f1，f2，…，fm，ζ1，ζ2，…，ζm]T為系統(tǒng)的模態(tài)頻率和阻尼比組成的向量。

2.3 模態(tài)參數(shù)不確定性的量化

根據(jù)模態(tài)形式狀態(tài)空間模型式（18），模態(tài)參數(shù)中的頻率和阻尼比只與系統(tǒng)矩陣有關(guān)，而振型只與輸出矩陣有關(guān)。因此，模態(tài)參數(shù)函數(shù)Γ()與模型參數(shù)θˉ的對應(yīng)關(guān)系為

利用式（23）并結(jié)合模態(tài)形式狀態(tài)空間模型，然后通過2.1部分的CRB計算，就可以直接求得模態(tài)參數(shù)的協(xié)方差下界。由于極大似然法的漸近無偏估計、漸近一致性和漸近有效性特點，CRB可以用來近似成參數(shù)不確定性水平。

本文采用的基于狀態(tài)空間模型的運行模態(tài)及不確定性分析方法的流程圖見圖1。總的來說，該方法采用隨機狀態(tài)空間（SSI）方法得到狀態(tài)空間模型參數(shù)θ的初值，通過卡爾曼濾波和極大似然估計構(gòu)建對數(shù)似然函數(shù)，利用EM迭代算法更新系統(tǒng)參數(shù)求得滿足收斂條件的極大似然值θ^，從而完成模態(tài)參數(shù)的識別，并利用CRB界限近似計算參數(shù)的不確定性水平，實現(xiàn)了參數(shù)不確定性量化。然而，風(fēng)力機在環(huán)境激勵下采集的振動信號中，除了含有結(jié)構(gòu)的模態(tài)信息外，通常還會存在其他一些虛假模態(tài)，比如激勵模態(tài)、風(fēng)輪轉(zhuǎn)動的諧波模態(tài)以及傳感器電子干擾模態(tài)等。因此，對振動信號進行本文的模態(tài)分析方法識別出的模態(tài)信息，需進一步結(jié)合結(jié)構(gòu)模態(tài)信息的特點，對虛假模態(tài)加以去除。

圖1 基于狀態(tài)空間模型的運行模態(tài)及不確定性分析方法流程圖Fig.1 Flowchart of operating modal analysis and uncertainty quantification based on state-space model

3 應(yīng)用舉例

3.1 風(fēng)電機組現(xiàn)場實測

實測風(fēng)力發(fā)電機組如圖2（a）所示，為Nordex S70 1.5 MW三葉片水平軸風(fēng)力機，風(fēng)機輪轂高度65 m，塔架塔底直徑4.04 m，塔頂直徑2.96 m。采集振動響應(yīng)的四個單軸壓電加速度計安裝在塔內(nèi)的三個維護平臺上，詳細的安裝位置與方向如圖2（b）所示。一般認為，垂直于風(fēng)機掃風(fēng)平面的方向為塔架的前后方向，平行于風(fēng)機掃風(fēng)平面的方向為塔架側(cè)向方向。本次測試的風(fēng)機停機工況，傳感器方向大致與風(fēng)機塔架主軸振動方向?qū)R。考慮風(fēng)電機組兩種典型受力工況下，即停機時和正常運行發(fā)電時的結(jié)構(gòu)振動響應(yīng)，實測加速度響應(yīng)時程如圖3所示，采樣頻率為100 Hz。兩種測試工況下振動信號的奇異值（SV）頻譜圖如圖4所示，只截取前5Hz范圍展示。

圖2 實測風(fēng)力發(fā)電機組和加速度傳感器布置Fig.2 Tested wind turbine and layout of acceleration sensors

圖3 風(fēng)機在停機與運行工況下的實測動力響應(yīng)Fig.3 Measured responses of the tested wind turbine under parked and operating conditions

由文獻［22］可知，該風(fēng)機在停機與運行工況時，其塔架前后方向前兩階模態(tài)的頻率均在0.5 Hz和4.0 Hz附近［22］。從圖3和圖4中可以看出，風(fēng)機運轉(zhuǎn)發(fā)電時的振動幅值明顯大于停機工況，且兩種工況的SV頻譜圖除了塔架前兩階頻率對應(yīng)的0.5 Hz和4.0 Hz附近有明顯的峰值外，其他地方也出現(xiàn)了新的峰值。對于停機工況，1.0～2.0 Hz出現(xiàn)的峰值，這是葉片振動為主對應(yīng)的模態(tài)［30］；而運行工況下出現(xiàn)的新峰值主要對應(yīng)于風(fēng)機葉片旋轉(zhuǎn)造成的諧波頻率，即以1P（葉片旋轉(zhuǎn)頻率），3P，6P，…，3nP等葉片通過頻率出現(xiàn)。其中，3P葉片通過頻率與塔架的第一階振動頻率接近，這給準確估計風(fēng)機的模態(tài)信息帶來了困難。

圖4 測量動力響應(yīng)的SV頻譜圖Fig.4 SV spectrum of measured responses

3.2 濾波處理及初始條件設(shè)置

由于所測得的振動響應(yīng)數(shù)據(jù)點較多，為提高計算效率，對響應(yīng)數(shù)據(jù)先進行濾波處理和重采樣。另外，此處還介紹了在使用本文方法時一些初始條件的設(shè)置。

由于塔架的前兩階模態(tài)在風(fēng)機結(jié)構(gòu)性能評估中顯得更重要，因此本文主要識別塔架兩個主軸振動方向的前兩階模態(tài)。塔架前兩階頻率根據(jù)上一節(jié)已知在5 Hz以內(nèi)，而振動信號的原始采樣頻率為100 Hz，信號范圍帶寬較大。因此，本文采用低通濾波器對測得的加速度響應(yīng)先進行濾波，設(shè)置的截斷頻率為5 Hz，然后對濾波數(shù)據(jù)進行重采樣，重采樣的頻率為10 Hz。

在用隨機狀態(tài)空間模型識別的過程中，模型階次n的確定也影響該方法的準確性和效率。這里，對狀態(tài)空間模型階次通過EM迭代過程來確定。首先，通過查看圖4中停機工況與運行工況的峰值數(shù)，初步確定一個模型階次值；其次，通過對殘差信息ek進行白噪聲驗證，如不滿足，則通過逐步增加模型階次的思路；最后，找出合適的階次值。針對實測風(fēng)機情況，從圖4中可知，風(fēng)電塔停機與運行工況下SV頻譜圖中前5Hz范圍內(nèi)分別有9個、12個峰值，即停機與運行工況的初始模型階次分別取18、24。接著，對識別得到的殘差信息ek進行相應(yīng)的頻譜圖分析，如若其SV頻譜圖中有著明顯的峰值，則依次增加模型階次大小，直至其無明顯峰值，近似于白噪聲在頻域中的表現(xiàn)。結(jié)果發(fā)現(xiàn)，對于停機工況，當n=28時，其殘余信息SV頻譜圖，已無明顯峰值存在，可認為接近白噪聲，如圖5（a）所示；對于運行工況，當n=70時，其殘余信息SV頻譜圖，已經(jīng)近似于白噪聲在頻域中的表現(xiàn)，如圖5（b）所示。

圖5 停機和運行工況下殘差信息SV頻譜圖Fig.5 SV spectrum of residuals under parked and operating conditions

通過SSI方法得到的模型參數(shù)用作EM算法的迭代初值，然后對似然函數(shù)進行多次迭代。盡管理論上，迭代次數(shù)越多，其識別結(jié)果越精準，但這也將加大線下處理時間。采用目標收斂率ρ=1.0×10-5時，兩種工況下的EM迭代過程如圖6所示。

圖6 EM算法迭代曲線圖（目標收斂率ρ=1.0×10-5）Fig.6 EM iteration process（target convergence rate ρ=1.0×10-5）

3.3 模態(tài)分析結(jié)果

根據(jù)上節(jié)的濾波處理以及初始條件設(shè)置，即可對測得的加速度響應(yīng)進行模態(tài)參數(shù)識別，然后再應(yīng)用CRB界限近似計算模態(tài)參數(shù)的不確定性水平。

考慮到在停機工況下傳感器布置方向與塔架主軸方向基本一致，原先為了快速比較兩個方向模態(tài)參數(shù)的差異性，將測得的x方向兩條響應(yīng)數(shù)據(jù)與y方向的數(shù)據(jù)分開分別進行模態(tài)參數(shù)識別（分開識別）。但考慮到風(fēng)機塔架為一整體結(jié)構(gòu)，其x方向與y方向必然存在著耦合作用，且參數(shù)識別時使用的加速度響應(yīng)越多，其識別結(jié)果越精確，因此又將停機工況測得的四條加速度響應(yīng)放在一起進行模態(tài)識別（一起識別），并將兩種識別方式結(jié)果進行了比較。表1和表2列出了上述兩種識別情況的模態(tài)參數(shù)和不確定性結(jié)果。從表1和表2可知，兩種方式識別出的前兩階頻率和阻尼比雖存在一些差異，但數(shù)值比較接近。這說明停機時的傳感器布置方向是與塔架振動方向比較一致的，這種情況下分開識別和一起識別對模態(tài)參數(shù)的估計值影響不大。然而，通過比較這兩種情況的不確定性水平，可以看出一起識別的模態(tài)參數(shù)變異系數(shù)大部分小于分開識別的，尤其表現(xiàn)在阻尼比方面。這也說明了將四條停機數(shù)據(jù)放在一起，包含了更多的振動信息，必然有助于識別的準確性。同時，特別隨著風(fēng)向改變，風(fēng)機對風(fēng)輪轉(zhuǎn)動后風(fēng)機工作，造成傳感器方向與塔架振動方向不一致，需要所有振動信號放在一起識別，并結(jié)合振型才能準確確定風(fēng)機實際的前后振動與側(cè)向振動模態(tài)。因此，風(fēng)機模態(tài)識別結(jié)果應(yīng)采用一起識別的方式識別的為準。

表1 停機工況下分開識別的模態(tài)參數(shù)及不確定性Table 1 Modal parameters and uncertainties identified separately under parked condition

根據(jù)表2的停機工況識別結(jié)果可知，風(fēng)機停機工況下，塔架側(cè)向的前兩階頻率0.49 Hz和4.07 Hz均大于前后方向0.48 Hz和3.83 Hz，且第二階頻率的差異更明顯。側(cè)向的前兩階阻尼比2.37%和0.98%同樣均大于前后方向0.65%和0.72%，且第一階阻尼比的差異更明顯。風(fēng)機塔架振動為主的不同方向的模態(tài)阻尼比的不同，主要與葉片槳距角有關(guān)。受空氣動力阻尼的影響，葉片拍打振動方向的阻尼比往往比揮舞振動方向的阻尼比大一些。另外，從前兩階模態(tài)參數(shù)對應(yīng)的變異系數(shù)可知，阻尼比的不確定性遠大于頻率，這也符合一般模態(tài)分析的現(xiàn)象。

表2 停機工況下一起識別的模態(tài)參數(shù)及不確定性Table 2 Modal parameters and uncertainties identified together under parked condition

風(fēng)機運行工況下的識別結(jié)果如表3所示，可以發(fā)現(xiàn)風(fēng)機風(fēng)輪轉(zhuǎn)動的情況下，塔架前后方向的第一階頻率0.50 Hz略大于側(cè)向的0.49 Hz。此時，前后方向的第一階阻尼比為2.27%，也高于側(cè)向的0.53%。這也符合之前停機時的結(jié)論，因為風(fēng)機運行時槳距角調(diào)整是葉片對風(fēng)，葉片拍打方向轉(zhuǎn)變?yōu)樗砬昂蠓较?，所以受氣動阻尼影響，前后方向的一階阻尼比偏大。而對于塔架的第二階模態(tài)，從表中可以看出，運行工況下，風(fēng)電塔前后方向的第二階頻率3.97 Hz，低于側(cè)邊方向4.19 Hz，與停機工況下識別出的第二階頻率情況相似。同樣，從前兩階模態(tài)參數(shù)對應(yīng)的變異系數(shù)都非常小可知，風(fēng)機運行模態(tài)分析的結(jié)果具有較好的準確性。

表3 運行工況下一起識別的模態(tài)參數(shù)及不確定性Table 3 Modal parameters and uncertainties identified together under operating condition

總的來說，比較表2和表3風(fēng)機停機與運行工況的模態(tài)參數(shù)結(jié)果，可以看出同一工況下風(fēng)機塔架的前后與側(cè)向模態(tài)參數(shù)的確存在著差異性。風(fēng)機運行工況下模態(tài)頻率均比對應(yīng)的停機工況的結(jié)果要大，這是由于風(fēng)葉轉(zhuǎn)動時產(chǎn)生了巨大的慣性力，導(dǎo)致風(fēng)機的有效質(zhì)量減少、葉片的剛度增大，使其整體剛度增大，根據(jù)剛度與頻率的關(guān)系可知，其頻率也會增大。對于運行工況下的塔身前后方向第一階阻尼比2.27%遠大于停機工況下的0.65%，以及對于停機工況下塔身側(cè)向的第一階阻尼比2.37%遠大于運行工況的0.53%，這是因為氣動阻尼主要影響葉片拍打方向的模態(tài)導(dǎo)致的［31］。

4 結(jié)論

本文介紹了一種基于狀態(tài)空間模型的運行模態(tài)及不確定性分析的方法。通過卡爾曼濾波構(gòu)建模型參數(shù)的似然函數(shù)，利用EM迭代算法更新模型參數(shù)求其極大似然值，進而估計模態(tài)參數(shù)，并利用CRB界限近似計算參數(shù)的不確定性水平，實現(xiàn)參數(shù)不確定性量化。在大型風(fēng)力機帶有旋轉(zhuǎn)葉片諧波影響的運行模態(tài)分析中，此方法可有效解決模態(tài)參數(shù)識別困難的問題，并得到模態(tài)參數(shù)估計的不確定性水平，為后續(xù)決策提供科學(xué)依據(jù)。本文提出的針對葉片旋轉(zhuǎn)的風(fēng)機運行模態(tài)分析方法，適用于葉片基本是正常勻速轉(zhuǎn)動的條件下，不適用于啟停階段或葉片轉(zhuǎn)速變化大的運行時段。另外，不同風(fēng)況下風(fēng)機自身的控制策略會調(diào)節(jié)風(fēng)機正常運行的葉片勻速轉(zhuǎn)動的速度，不影響本文方法的應(yīng)用，但是由于氣動阻尼與結(jié)構(gòu)振動的幅度和速度有關(guān)，因此不同風(fēng)況下結(jié)構(gòu)的氣動阻尼會受到影響。通過實際大型風(fēng)力機在停機和運行兩種工況的模態(tài)分析，著重分析了風(fēng)機塔架不同方向模態(tài)參數(shù)的差異性，以及不同工況下模態(tài)參數(shù)及不確定性的變化。分析結(jié)果表明：

（1）風(fēng)機葉片為細柔不對稱構(gòu)件，氣動阻尼影響明顯，且氣動阻尼主要影響葉片拍打方向的塔架第一階模態(tài)阻尼；

（2）由于風(fēng)輪轉(zhuǎn)動的影響，風(fēng)機運行工況下的模態(tài)頻率略大于對應(yīng)的停機工況下的模態(tài)頻率；

（3）無論風(fēng)機是在停機還是運行工況下，塔架側(cè)向振動方向的第二階模態(tài)頻率一般大于塔架前后振動方向的第二階模態(tài)頻率。