基于狀態空間模型的大型風力機運行模態及不確定性分析

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

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

0 引言

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

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

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

本文采用一種基于狀態空間模型的運行模態及不確定性分析方法，以解決帶有葉片旋轉的風機結構的模態參數識別問題。最后，本文以實例分析一風機不同塔架方向的模態參數的差異性，以及不同工況下風機模態參數的變化情況，并結合識別結果的不確定性水平，證明該方法的實用性。

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

1.1 隨機狀態空間模型

在環境激勵下，長度為N的結構響應數據Y={y1，y2，…，yN}，可用離散時間隨機狀態空間模型進行模擬

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

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

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

隨機狀態空間模型式（1）的參數化形式具體可表示成

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

1.2 模型參數極大似然估計

系統參數θ可通過觀測響應數據Y進行估計。由于極大似然估計方法具有漸近無偏估計、漸近一致性和漸近有效性等優異特性，本文采用極大似然法來實現隨機狀態空間模型參數θ的估計，進而獲取結構的運行模態參數。

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

式中，

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

根據極大似然法的原理可知，系統參數θ的極大似然估計可表示為

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

1.3 最大期望(EM)算法

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

由貝葉斯公式可知：

兩邊同時取對數可得：

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

式中，ρ為目標收斂率。

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

其中，

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

2 模態參數的不確定性分析

2.1 克拉美羅下界(CRB)

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

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

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

針對隨機狀態空間模型的似然函數式（6），有［28］

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

2.2 模態形式狀態空間方程

由于結構頻率和阻尼比只與A的特征值有關，將一般形式的隨機狀態空間方程轉換為模態形式后計算FIM，可直接獲得模態參數的協方差信息。對隨機狀態空間方程式（1）做相似變換，令xk=，可轉換為模態形式的狀態空間方程：

式（19）中每個復模態分塊矩陣與結構頻率和阻尼的對應關系為

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

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

其中，[1 0]或1代表了每一模態振型最大幅值所對應的位置。

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

此時，模態形式的狀態空間模型式（18）的參數變量可寫成

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

2.3 模態參數不確定性的量化

根據模態形式狀態空間模型式（18），模態參數中的頻率和阻尼比只與系統矩陣有關，而振型只與輸出矩陣有關。因此，模態參數函數Γ()與模型參數θˉ的對應關系為

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

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

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

3 應用舉例

3.1 風電機組現場實測

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

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

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

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

圖4 測量動力響應的SV頻譜圖Fig.4 SV spectrum of measured responses

3.2 濾波處理及初始條件設置

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

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

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

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

通過SSI方法得到的模型參數用作EM算法的迭代初值，然后對似然函數進行多次迭代。盡管理論上，迭代次數越多，其識別結果越精準，但這也將加大線下處理時間。采用目標收斂率ρ=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 模態分析結果

根據上節的濾波處理以及初始條件設置，即可對測得的加速度響應進行模態參數識別，然后再應用CRB界限近似計算模態參數的不確定性水平。

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

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

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

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

風機運行工況下的識別結果如表3所示，可以發現風機風輪轉動的情況下，塔架前后方向的第一階頻率0.50 Hz略大于側向的0.49 Hz。此時，前后方向的第一階阻尼比為2.27%，也高于側向的0.53%。這也符合之前停機時的結論，因為風機運行時槳距角調整是葉片對風，葉片拍打方向轉變為塔身前后方向，所以受氣動阻尼影響，前后方向的一階阻尼比偏大。而對于塔架的第二階模態，從表中可以看出，運行工況下，風電塔前后方向的第二階頻率3.97 Hz，低于側邊方向4.19 Hz，與停機工況下識別出的第二階頻率情況相似。同樣，從前兩階模態參數對應的變異系數都非常小可知，風機運行模態分析的結果具有較好的準確性。

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

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

4 結論

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

（1）風機葉片為細柔不對稱構件，氣動阻尼影響明顯，且氣動阻尼主要影響葉片拍打方向的塔架第一階模態阻尼；

（2）由于風輪轉動的影響，風機運行工況下的模態頻率略大于對應的停機工況下的模態頻率；

（3）無論風機是在停機還是運行工況下，塔架側向振動方向的第二階模態頻率一般大于塔架前后振動方向的第二階模態頻率。