基于粒子濾波算法的風力發電塔地震動力響應預測

徐亞洲， 任倩倩， 于明陽， 時文浩

(西安建筑科技大學 土木工程學院，西安 710055)

隨著電子計算機技術的發展，有限元數值分析模型和結構動力試驗成為結構動力學中愈加成熟的領域，因此動力學模型更新和對動態系統結構性能的合理預測仍是非常重要的研究方向[1-2]。利用測量系統響應進行動力模型更新，在結構振動控制、結構健康監測以及可靠性評估等方面都有廣泛的應用，但在建立結構計算模型及其激勵的數學模型的過程中，可能存在由于外部原因導致的建模誤差及參數不確定性影響，因此能夠準確和適當地量化此類不確定性對于結構動力模型的更新至關重要[3-5]。土木工程結構中的有限元模型往往與許多來源的不確定性和建模誤差有關，尤其是對于復雜的土木工程結構，在建模過程中有限元模型相比實際結構進行了大量的假設和簡化，當此類有限元模型用于預測一些感興趣的量時，這些不確定性因素會影響土木工程結構的設計和性能評估[6-8]。這些不確定性因素促使研究人員通過概率有限元模型更新方法來考慮潛在的結構不確定性。

基于經典概率的理論分析模型是目前工程中常用的不確定性分析模型[9]，但在統計信息不充分的條件下，經典概率分析往往具有一定的局限性。近年來，貝葉斯估計理論由于克服了經典概率理論中難以處理的小樣本問題成為概率不確定性分析研究的新熱點[9]，貝葉斯的基本思路是根據狀態的先驗信息和測量信息來不斷遞推更新狀態的后驗概率密度函數[10]。毛文貴等[11]運用基于遺傳智能采樣技術的改進貝葉斯理論，可以很好地處理軸承轉子系統不平衡識別中的模型不確定性參數的影響。Pepi等[12]建立了位于特爾尼的彎曲斜拉橋有限元模型，該有限元計算模型考慮了模型不確定性以及對連接和構件的有限認識，在貝葉斯框架下進行環境振動試驗以確定用于模型更新的主要動態參數，并用于評價基于有限元模型計算結果的準確性。Sedehi等[13]在貝葉斯層次模型的基礎上，提出了漸近近似和最大后驗估計方法，該方法可以有效解決由外界因素引起的參數不確定性的影響。

在貝葉斯理論估計計算框架中，由于粒子濾波(particle fliter,PF)方法具有較好的處理復雜系統的能力而被廣泛研究及應用[14]。萬春風等[15-17]基于粒子濾波算法對單自由度和多自由度結構進行參數損傷識別研究，指出了粒子濾波算法在參數識別方面相較于其他濾波方法的優越性。一些學者[18-21]采用粒子濾波方法在線預測疲勞裂紋擴展情況并更新剩余壽命，并通過試驗結果驗證了粒子濾波對疲勞裂紋擴展和剩余壽命在線預測的可行性和合理性。樊學平等[22-23]采用粒子濾波算法建立了隨時間更新的動態模型，并結合極值監測數據實現橋梁結構可靠度指標的動態預測。

風力發電塔對動力荷載的高敏感性引起了研究者的關注[24-26]，且目前關于風力發電塔的抗震性能研究多采用有限元模擬的手段，但在構造結構數值分析模型的過程中，總是存在簡化模型所帶來的建模誤差和參數不確定性。振動臺試驗作為目前獲得結構動力響應分析的主要手段，被廣泛應用于結構性能評價與分析中，試驗結果對于校正能夠準確描述結構響應的有限元數值分析模型至關重要。因此本文提出了一個用于結構動力響應預測的概率貝葉斯估計計算框架，即首先在貝葉斯估計原理基礎上，建立地震激勵下結構有限元分析模型。然后采用粒子濾波方法，充分利用結構先驗信息并結合振動臺試驗觀測的風力發電塔動力響應數據信息選擇重要性函數，實現風力發電塔動力響應預測，并以有限元數值解與試驗值之間的差值來量化有限元模型的不確定性誤差。最后通過風力發電塔振動臺試驗結果，對計算方法的合理性及有效性進行驗證。

1 粒子濾波算法

在存在結構建模不確定性的情況下，需要確定并建立合理的結構模型，并基于實測結果進行加權修正。粒子濾波算法具有顯著處理參數不確定性估計的優越性，目前被廣泛用于結構損傷識別和結構不確定性分析中。該算法主要是應用蒙特卡洛方法產生大量的隨機樣本(粒子)，然后結合觀測值利用樣本權重預測下一時刻的系統狀態或參數概率分布，最后根據算法不斷遞推更新得到結構參數的最優貝葉斯估計值。

與其他濾波方法一樣，粒子濾波算法需要基于估計問題構建狀態空間模型，主要由狀態方程、觀測方程和初始狀態信息組成。狀態方程描述了系統狀態量隨時間的演變過程，觀測方程描述了系統適時輸出與系統當時狀態的模型，函數表達式如下

狀態方程為

Xk=f(Xk-1)+ωk

(1)

觀測方程為

Zk=h(Xk)+νk

(2)

初始狀態信息為

p(X0Z0)=p(X0)

(3)

式中:Xk和Xk-1分別為k時刻和k-1時刻的狀態值；Zk為k時刻的觀測值；f(·)為表達系統當前時刻與前一時刻的變化關系的已知非線性函數；h(·)為狀態值和觀測值之間的關系非線性函數；wk為過程噪聲，描述了系統的各種不確定性影響因素；vk為觀測噪聲，代表了觀測過程中的不確定性觀測誤差。

1.1 貝葉斯估計理論

貝葉斯的基本步驟包括預測和更新[27]，預測為利用模型前一時刻的先驗已知信息對當前時刻進行預測，更新是結合當前時刻的實際觀測進行遞推，從而獲得待估值在當前時刻的后驗概率密度，在貝葉斯更新過程中，狀態變量的更新是為了更好地表示真實的結構，其中更新過程是由先驗已知信息和所研究的結構測量中包含的測量信息完成的。狀態量中包含了所有的不確定性，因此，假設前一時刻的后驗概率p(Xk-1Z1:k-1)已知，則基于貝葉斯估計算法，可得到后驗概率密度p(XkZ1:k)的估計式

狀態預測方程

(4)

狀態更新方程

(5)

式(4)和式(5)組成了貝葉斯估計的計算框架。首先若已知k-1時刻的狀態值Xk-1，則可通過式(1)預測得到k時刻Xk的先驗概率密度p(Xk-1Z1:k-1)，基于觀測方程(式(2))可得到似然概率p(ZkX1:k)。在貝葉斯計算中，直接計算后驗概率密度函數牽涉到復雜的積分計算，對于非線性動態系統的狀態估計問題求解困難。因此，引入了基于隨機抽樣計算的蒙特卡洛方法來計算后驗概率。

1.2 蒙特卡洛原理

基于大數定理，蒙特卡洛方法[28]將所求解的隨機問題發生的概率用大量試驗中發生的頻率來估計，當統計數據比較大時，可以認為該隨機事件發生的頻率能無限接近期望值，從而可使貝葉斯后驗概率中的積分運算用大量樣本點的平均回報來逼近。

如對本研究來說，可以在狀態的后驗概率密度中大量的隨機采樣，則狀態的后驗概率密度可表達如式(6)。作為一種復雜問題有效的處理求解方式，該方法回避了工程應用中積分計算的困難，用狀態發生的頻率來表示概率，這在積分運算中有著非常重要的意義。

(6)

式中，δ(·)為狄拉克函數。然后待估狀態的期望值計算如式(7)所示

(7)

1.3 粒子濾波算法

(8)

作為本文研究問題所要求解的對象，故不能在待求解的后驗概率中采樣。一般采取從已知重要性密度函數q(XkZ1:k)中選取的N個樣本近似表示

(9)

其中,ω(Xk)為樣本權重，即

(10)

依據式(10)，式(9)可進一步化簡為

(11)

其中，歸一化后的權重值為

(12)

基于馬爾可夫假設與簡化，歸一化后粒子的權值(重要性權值)在每個時間步長上有如下遞推關系

(13)

通過權重值的不斷計算，這些權重被附加到每個相應的粒子上。但是這些粒子并沒有被修改，只是它們的相對權重被改變了，也就是說，當所有的粒子偏離觀測值時它們不會被拉回觀測值，重要性采樣只改變了它們的相對權重。這個缺點可能會導致單個粒子的權重占比非常大，而其他粒子的權值很小。因此，需要進行重采樣步驟，定義有效粒子數量為

(14)

當有效粒子的數量小于采樣數量N時，則表示粒子樣本退化嚴重，需要進行重采樣，即高重要性權值的粒子會被多次復制，低權值的粒子少量復制或被去除以保持粒子數量不變。重采樣后每個后驗粒子的權值相同，均為1/N。

重要性密度函數一般選取為狀態轉移密度函數，即

(15)

將式(15)代入式(13)，可重新得到權值更新方程的表達式

(16)

2 動力響應預測的粒子濾波實現

建立有限元模型進行風力發電塔動力響應分析的過程中，由于理想假設或簡化容易導致參數不確定性或模型不確定性，繼而導致有限元模擬值與實際測量值之間存在較大區別。因此本章基于第1章關于粒子濾波算法的計算原理的介紹，利用粒子濾波算法中的狀態空間模型將這些不確定信息考慮進來，作為先驗已知信息實現結構動力響應的預測過程。

2.1 力學模型

對于建筑結構，常用的力學模型分為層模型、桿模型及桿系-層間剪切模型。根據風力發電塔有限元分析方法的需要，目前有4種比較常見的建模方式。第一種為IEC規范[29]中所建議的單自由度模型簡化方法，即將塔體頂部葉輪、機艙的質量及塔筒50%的質量集中到塔筒頂部；第二種為我國GB 50011—2010《建筑抗震設計規范》[30]中常用的建議多自由度模型簡化方法，即將塔體頂部葉輪、機艙的質量集中到塔體頂部質量點處，各個筒段的質量分別集中到對應質量點處；第三種為簡化的有限元模型[31]，即葉輪、機艙采用剛體單元，將質點位置按照剛體運動耦合到塔筒頂部中心；第四種為精細化有限元模型[32]，即將葉片按照實際結構原尺寸建模從而考慮葉片與塔體相互作用對地震作用下結構性能的影響，每個葉片單元的質量集中在一個參考點，輪轂和機艙用集中質量點表示，通過剛性連接連接到塔頂部分。本節采用的第二種建模方法，即將風力發電塔簡化為多自由度層剪切結構進行分析。

以剪切型動力學模型為例，地震激勵下其動力運動方程為

(17)

結構的阻尼采用瑞麗阻尼，即

C=αM+βK

(18)

式中，α，β為阻尼系數。

2.2 動力方程求解

在結構動力方程的求解過程中，Newmark-β法應用較為廣泛，其主要思想是假定ti～ti+1時間段內的加速度變化呈線性規律，其假定為

(19)

式中，a，b為可調整參數。當a=1/2，b=1/6時，即為線性加速度法；當a=1/2，b=1/4時，即為平均常加速度法；當a=1/2，b=0時，即為中心差分法。本節采用線性加速度法計算。由式(19)可求得t+Δt時刻的速度和加速度

(20)

其增量表示形式為

(21)

將式(21)代入式(17)，可得到

(22)

通過求解式(22)可得到t+Δt時刻的位移，將其代入式(20)即可求得結構加速度及速度響應。

2.3 結構動力響應預測實現

基于試驗觀測數據，粒子濾波算法的結構動力響應預測過程的主要步驟如下。

步驟1首先需要構造狀態空間模型，根據Newmark-β方法求得有限元簡化模型結構的動力響應數值解作為趨勢項，用來構造本文的狀態方程；使用監測儀器測得結構的實際動力響應，然后利用狀態量及振動臺試驗動力響應觀測量的信息構造觀測方程；初始狀態先驗信息基于有限元數值計算結果進行概率統計給出。

步驟2基于狀態空間模型，根據初始先驗信息的概率分布函數，利用蒙特卡洛抽樣方法可得到狀態變量的初始參數樣本，即初始先驗粒子{x0,i,i=1,2,3,…,N}。

步驟3在每一個時間步長上進行粒子更新，利用狀態方程可獲得下一個時刻的后驗粒子

(i=1,2,3,…,N)

步驟4基于k時刻的實際測量值Zk，可通過似然概率密度函數更新每個后驗粒子的權值，此權值衡量了該粒子產生的模擬數據接近實際觀測數據的近似程度。似然概率密度函數可近似為正態分布[33]，即

(23)

式中:Zmn,i為第i個粒子的模擬觀測數據；R為量測噪聲方差，可通過狀態量及觀測值的數據信息近似估計。

步驟5根據式(12)進行權重歸一化。

步驟6根據式(14)計算有效粒子的數量，當有效粒子的數量小于初始采樣數量N時，則進行重采樣，重采樣后每個粒子的權值相同，均為1/N。

以此不斷遞推更新，求得結構響應的最優預測值。

具體流程圖如圖1所示。

圖1 粒子濾波算法預測過程流程圖

3 試驗驗證

本文開展了2 MW陸上風力發電塔振動臺試驗，基于狀態空間模型和試驗觀測值對有限元計算結果的不確定性進行修正，說明本文所提用于結構動力響應預測的粒子濾波框架的合理性。

3.1 試驗模型信息

根據量綱協調原理和結構動力方程[34-35]確定相似關系，按照相似原則及抗彎剛度等效設計模型。試件與原型結構的幾何相似比為1/20，模型材料為Q345鋼。塔呈錐形，直徑從底部的200 mm線性減小到頂部的130 mm。塔筒為4段，段與段之間用10 mm厚法蘭板連接，每段高度從下到上分別為1 m,0.85 m,1 m,1 m。由于薄板在軋制和焊接過程中制作困難，根據抗彎剛度等效原則，4個管段對應的管壁厚度分別為4 mm,3 mm,3 mm和3 mm。機艙和轉子的總質量為157.5 kg，4段塔筒的質量分別為85 kg，58 kg，51 kg，35 kg。試驗模型的示意圖如圖2所示。

按照抗震設計規范中輸入地震波的選取原則，從太平洋地震工程研究中心[36]數據庫選取了Chi-Chi波、Westmorland波、El-Centro波和Taft波作為模型結構振動臺輸入地震激勵，將垂直于葉片的方向定義為地震激勵輸入的方向，圖3給出了4條地震動記錄水平方向加速度時程。將地震動加速度峰值(peak ground acceleration，PGA)統一調幅至0.07g，其加速度反應譜如圖4所示。

(a) 試驗模型

(a)

(c)

地震波的輸入工況分別設置為0.07g，0.2g和0.4g。在不同幅值地震動輸入前后采用PGA為0.03g的白噪聲激勵來識別試驗模型的動力特性。本試驗加速度相似比為2，頻率相似比為6.324 6，故實際臺面輸入地震激勵需按照相似關系進行調整。

在風力發電塔結構抗震設計中，塔頂的動力響應是抗震需求的重要指標，因此通過合理布置測點測量了塔頂加速度響應及位移響應。不同地震動強度幅值輸入前后利用白噪聲掃頻，并對結果進行處理得到模型結構的動力特性。利用傳遞函數對試驗模型的固有頻率和振型進行識別，得到結構前兩階自振頻率分別為3.22 Hz,22.46 Hz，結構一階阻尼比為1.5%。

圖4 地震動記錄的加速度反應譜

3.2 模型動力響應預測

基于本文方法，對模型結構動力響應進行預測。

有限元數值分析：基于2.1節的介紹，本研究將風力發電塔動力系統簡化為5個自由度的剪切動力學模型來描述，如圖2(b)所示，根據質量點處的截面尺寸計算出截面抗彎剛度。基于MATLAB有限元軟件，計算5個自由度簡化模型的質量和抗彎剛度矩陣。根據結構實測的阻尼比及自振頻率，采用瑞麗阻尼計算阻尼矩陣，最后利用2.2節中介紹的Newmark-β法求得有限元模型的塔頂加速度響應和位移響應。

狀態空間方程建立：根據12個地震工況作用下求得的5個自由度結構的動力響應數值計算結果，選用塔頂加速度響應和位移響應為狀態待估量，可估計出每個地震工況作用下的狀態噪聲方差；利用數值計算結果及對應工況下的振動臺試驗動力響應觀測量信息估計出測量噪聲方差；基于數值計算結果進行概率統計給出狀態量的初始信息，并基于初始狀態信息選取200個粒子。利用第2章中介紹的粒子濾波算法框架進行風力發電塔動力響應預測。

以振動臺試驗測得的觀測值為參考，基于粒子濾波算法對有限元模型結果進行修正預測。定義考慮和不考慮觀測值修正后的計算值與試驗實測值之間差值的絕對值作為動力響應的計算偏差。則考慮觀測值修正后(PF)，和不考慮觀測值修正(finite element method,FEM)的加速度和位移響應的偏差最大值計算結果，如圖5所示。結果表明有限元計算的動力響應偏差隨地震動PGA幅值的增大而增大。以Chi-Chi地震動激勵下塔頂加速度響應為例，當PGA值為140 gal時，加速度偏差為2.21 m/s2，而PGA值為800 gal時，加速度偏差為14.89 m/s2。考慮觀測值修正之后有限元計算結果偏差顯著減小，以塔頂位移響應為例，考慮觀測值修正之后位移偏差減小近一倍。

(a) 塔頂加速度響應偏差

(b) 塔頂位移響應偏差

為進一步說明本文方法的有效性及可行性，選取了El-Centro波和Chi-Chi波激勵作用下剪切動力學模型結構的加速度響應時程進行分析，圖6給出了粒子濾波算法修正前后加速度時程曲線對比。同一地震激勵下，隨輸入地震動幅值的增大，有限元計算結果相對于試驗值的偏差逐漸增大。考慮觀測值修正后粒子濾波預測的結果精度較好，預測值與試驗實測值近似相等，具有很好的一致性，進一步驗證了本文所提校正不確定有限元分析模型的概率貝葉斯框架和粒子濾波預測算法的合理性，表明粒子濾波算法可以利用有限元分析結果不確定性的所有先驗信息，從而實現對包含變異性結構動力響應的穩步預測。

4 結 論

本研究基于粒子濾波算法，提出了一個用于校正不確定有限元數值模型的概率貝葉斯計算框架對有限元計算的動力響應結果進行預測。利用某2 MW風力發電塔縮尺模型的振動臺試驗驗證了此框架的有效性和可行性。從分析結果可以得到以下結論：

(1) 有限元計算結果具有顯著的不確定性，且其計算誤差隨地震輸入幅值的增大而增大，考慮觀測值修正之后預測結果不確定性顯著減小。

(a)

(d)

(2) 粒子濾波模擬預測的動力響應值精度較高，模擬預測值與試驗實測值的時程曲線近似相等，變化規律具有很好的一致性。

(3) 粒子濾波算法與試驗觀測手段相結合，可以完成高效的結構動力響應預測，繼而為結構不確定性及可靠度分析提供參考。

結構不確定性分析應成為工程結構分析的一個重要部分，粒子濾波算法可以利用有限元分析結果不確定性的所有先驗信息，從而實現對結構動力響應的穩步預測。