考慮參數不確定性的海上風電結構模型修正方法

柳振海，李 剛，劉 博，徐明強，呂 晴，蔣 上

1.(中國電力工程顧問集團有限公司，北京 100120；2.中國海洋大學，山東 青島 266100 )

0 引言

海上風電作為一種清潔能源，越來越受到世界各國的重視，大力發展海上風電，可優化能源結構，促進經濟、社會、環境協調發展。然而，海上風電設施長期服役于復雜的海洋環境中，一旦發生故障或損傷，不僅影響風電運行，也會造成巨大的經濟損失，開展結構健康監測勢在必行。獲取準確的有限元模型是結構健康監測的重要內容，但由于建模誤差、參數不確定性等問題，有限元模型不一定能準確反映結構的真實情況，需要對有限元模型進行修正，使得修正后的模型能準確反映真實結構的動力響應特征。

有限元模型修正方法可以分為確定性模型修正和不確定性模型修正。前者將結構參數及響應假設為恒定不變的，修正結果為唯一確定解。但在實際工程設計建設中，結構部件的實際加工組裝可能與設計值存在一定誤差，導致結構參數的不確定性。因此確定性模型修正得到的最優解只是不確定性模型修正解中的一個特例。隨機模型修正是對傳統模型修正方法的創新，并考慮到了不確定性因素的影響，利用統計概率分析量化不確定性因素，將不確定模型修正問題轉化為均值和標準差的修正問題，從而得到與實際結構動力響應特征一致的有限元模型[1-3]。

由于能夠代替復雜的有限元模型分析，代理模型被廣泛用于模型修正。許澤偉[4]等提出一種基于多項式混沌展開和KL 散度的隨機有限元模型修正方法，以三維桁架為例，對彈性模量和密度的均值和標準差進行修正。冷建成[5]等以某海洋平臺結構為研究對象，提出了Kriging模型與多目標遺傳算法優化相結合的動力學模型修正方法，實驗顯示能夠顯著改善Kriging 模型精度。孫永朋[6]等針對隨機模型修正精度和效率低的問題，提出一種基于Kriging 模型和小波包能量譜的隨機有限元模型修正方法。丁雅杰[7]等提出一種基于貝葉斯推理的非線性結構模型修正方法，同時考慮激勵的隨機性，建立了復合隨機振動系統的動力可靠度分析方法。

本文針對小樣本數據的模型修正提出了一種基于KL 散度的分步型隨機模型修正方法，通過數值模擬驗證了此方法的有效性。首先采用靈敏度分析法選擇待修正參數，通過正交試驗生成訓練樣本，利用有限元模型計算樣本響應，從而構造以待修正參數為輸入，以頻率響應為輸出的Kriging 代理模型；利用蒙特卡洛模擬得到初始樣本點并計算其響應值，以樣本點響應與實測數據的KL 散度構造目標函數，通過多目標優化算法對海上風電結構進行隨機模型修正。

1 隨機模型修正方法

1.1 Kriging模型

相較于其他響應面模型(如多項式回歸、神經網絡等)，Kriging 模型在構造過程中只需要較少的樣本點，計算成本較低，并且處理非線性數據效果優異[8]，能實現對有限元模型響應較準確預測，更適合處理小樣本實測數據的模型修正。

Kriging 模型是一種通過已知樣本點信息估計未知試驗點信息的無偏估計模型，采用Kriging 代理模型的前提是假設所有數據之間都服從n維正態分布[9]。

已知樣本點x= [x1,x2, ...,xm]，xi是n維向量，對應響應為y(x) = [y(x1),y(x2), ...,y(xm)]，由于x服從正態分布，因此y服從均值為μ、協方差C的多維高斯分布。y(x)可以用式(1)表示：

式中：P代表基函數的個數；f(x)=[f1(x)，f2(x)，…，fl(x)，…，fp(x)]∈Rp×p為基函數矩陣；β=[β1,β2…,βl,…,βP]為基函數權重；fT(x)β為線性回歸模型，z(x)為服從高斯分布的隨機過程；不同變量之間的協方差矩陣C為：

式中：R為數據點之間的協方差矩陣；σ2為方差。

樣本點的似然概率為：

式中：n為樣本數。

對式(3)取對數：

使得似然概率最大的超參數β，σ2的估計值和為：

對于建立的代理模型，還需要對其進行進一步的驗證。本文采用決定系數R2和相對均方誤差RMSE評價Kriging 代理模型精度，R2的取值范圍為0 ～1，R2值越大，表示預測效果越好，RMSE越趨近于0 表示預測效果越好，其表達式為

1.2 KL散度

KL 散度可以用來度量兩個概率分布之間的差異[10]，需要注意的是，KL 散度并不滿足對稱性，即DKL(P||Q) ≠DKL(Q||P)。假設實測頻率與樣本頻率分別為二維正態分布P和Q，此時DKL(P||Q)為前向KL 散度，主要表現為在任何P(x) ＞ 0 的位置使得Q(x) ＞ 0，即使得Q能夠最大程度上覆蓋P；DKL(Q||P)為反向KL散度，主要表現為在P(x)趨近于0 時，Q(x)也盡可能趨近于0。

在優化問題中，需要通過迭代優化樣本頻率Q使得兩個分布之間的KL 散度能夠盡可能小，從而實現用樣本頻率Q擬合實測頻率P，且實測頻率分布往往是多峰分布，前向KL 散度求解得到多個極小值對應的待修正參數的平均值，不能滿足優化要求，而采用反向KL 散度則可以確保至少達到一個局部最優解，因此本文采用KL 散度計算實測頻率與樣本頻率分布之間的差異程度。

假設P、Q服從均值為μ0，μ1，協方差矩陣為Σ0，Σ1的高斯分布，則兩分布之間的KL散度為：

式中：tr表示矩陣跡；det 表示某矩陣的行列式；k為Σ0的維度。

1.3 目標函數

在優化算法中，目標函數提供了對不同解進行評估和比較的指標，決定了算法的搜索目標，因此選擇適當的目標函數對解決實際問題至關重要。本文以待修正參數的均值和方差為修正目標，以實測頻率與樣本點頻率的KL 散度為目標函數，同時對參數的均值和方差進行修正。由于兩個分布沒有重疊時，KL 散度沒有意義，無法進行迭代操作，因此將目標函數分為兩個部分：首先通過樣本點頻率及實測頻率的均值、方差構造多目標函數進行優化計算，當兩個分布出現重疊時，改用KL 散度為目標函數，具體目標函數構造如下：

式中：nsim為樣本響應數量；nobs為實測值數量。

1.4 優化算法

在確定好目標函數后，需要選擇合適的優化算法求待修正參數的最優解。NSGA-Ⅱ[11]是目前最常見的多目標優化算法之一，在原始算法的基礎上，引入了快速非支配排序、擁擠度等概念，減少了計算的復雜度，并提高了算法效率，其基本思想為：

1)隨機生成規模為N 的初始種群Pt，通過非支配排序、選擇、交叉、變異，生成子代種群Qt；

2)將初始種群Pt與子代種群Qt合并得到新的種群Rt；

3)進行快速非支配排序，并對每個非支配層的個體進行擁擠度計算，根據非支配關系及個體的擁擠度選擇合適的個體組成新的父代種群Pt+1；

4)采用遺傳算法生成新的子代種群Qt+1，并將Pt+1與Qt+1合并生成新的種群Rt；

5)依此類推，直到滿足程序結束的條件。

2 海上風電結構數值模型

選取一個5 MW 海上風機作為研究對象，風機采用導管架基礎，基礎深入海底，底部高程為-50.001 m，頂部高程為10.262 m，上部塔筒高度100 m，采用變截面，塔筒底部直徑6 m，壁厚0.035 1 m，頂部直徑3.87 m，壁厚0.024 7 m，頂部機箱和葉片簡化為偏心質量點，共31.452 t。利用ANSYS 建立有限元模型，模型共166 個節點，221 個單元。對其進行模態分析，得到前6 階頻率為：0.490 16 Hz、0.490 28 Hz、1.661 4 Hz、1.661 8 Hz、3.774 2 Hz、3.774 3 Hz，前三階模態振型如圖1 所示。

圖1 海上風機前三階模態振型

3 海上風電結構隨機有限元模型修正

3.1 修正參數的選擇

考慮到海上風機長期受腐蝕、沖刷、溫度、濕度等因素影響，其中材料參數選擇彈性模量、密度為待選參數，幾何參數選擇基礎腿柱外徑、斜撐外徑、腿柱壁厚、斜撐壁厚為待選參數。對待選參數進行靈敏度分析，選擇高靈敏度修正參數：將各初選參數值增加2%，得到修改后的頻率值，計算各階頻率變化率，除以參數變化率即為各參數前兩階頻率的靈敏度S，如式(11)所示。

式中：f i為第i階初始頻率，Hz；為參數值增加后的第i階初始頻率，Hz；Pm為初始參數值；Pmc為增加后的參數值。

表1 為修正參數各階頻率靈敏度，由表1可知，在材料參數中，彈性模量與密度靈敏度相差不大，在幾何參數中，腿柱外徑靈敏度較高。由于采用材料參數可以有效地提高方差修正效果，因此最終彈性模量(E)、密度(Des)及腿柱外徑(OL)為待修正參數。假設待修正參數服從高斯分布，根據施工圖紙及工程經驗確定其初始均值、方差見表2 所列。

表1 修正參數各階頻率靈敏度

表2 修正參數初始均值與方差

3.2 建立代理模型

根據各參數的初始均值和方差選擇95%置信區間為修正參數取值范圍，即：彈性模量取值范圍為[205.941，206.059]，密度取值范圍為[7 849.4，7 850.6]，腿柱外徑取值區間為[1.191，1.309]。通過正交試驗設計選擇，取三因素五水平共25 組設計樣本參數，見表3 所列，將樣本參數代入有限元模型中計算得到前兩階頻率，由此構建Kriging 模型。

表3 正交試驗設計樣本

采用相對均方誤差RMSE與決定系數R2對建立的Kriging 代理模型進行有效性評價，計算結果見表4 所列，證明構造的Kriging 代理模型精度較高，可以代替有限元模型。

表4 Kriging代理模型有效性評價

3.3 數值模擬

假設有限元模型待修正參數試驗值與初始值見表5 所列，對試驗值進行蒙特卡洛抽樣1 000 次得到樣本點，通過Kriging 模型計算得到樣本點前兩階頻率響應作為實測數據。首先采用實測數據與每次迭代生成的隨機樣本頻率的均值和方差構造多目標函數，隨后利用NSGA-Ⅱ多目標優化算法進行求解，得到結構待修正參數的修正值，見表5 所列。

表5 參數修正前后誤差

表5 中標準差是實測數據與其平均數離差平方的算術平均數的平方根，反映實測數據的離散程度，誤差是試驗值與初始值之間的差異。根據修正參數最優解，通過構造的Kriging 模型進行計算，得到修正后有限元模型前五階頻率的均值和標準差，與初始值和試驗值進行比較，見表6、表7 所列，修正前后頻率均值的最大誤差由68.59%降低為0.1%，標準差的最大誤差由212.82%降低為17.29%，且前三階頻率標準差最大誤差為0.48%，修正結果表明本文方法有效。

表6 修正前后頻率均值

繪制修正前后模型與試驗模型頻率置信橢圓如圖2 所示，可以看出，修正前實測值與初始值頻率置信橢圓大小與中心點位置有明顯差異，修正后實測值與修正值的頻率置信橢圓的大小、中心點、傾斜方向幾乎完全相同。圖3給出了結構前兩階頻率修正前后的概率密度曲線，對比可知，修正前實測值曲線與初始值曲線的峰值有明顯差異，兩條曲線重合程度較低。修正后的實測值與修正值概率密度曲線形狀相似且幾乎完全吻合。進一步驗證了本文所提修正方法的有效性。

圖2 修正前后模型與試驗模型頻率置信橢圓

圖3 修正前后頻率概率密度曲線

4 結語

考慮到參數隨機不確定對結構響應分布的影響以及KL 散度在不確定性度量的優勢，本文采用Kriging 模型代替有限元模型進行計算，并通過靈敏度分析法選擇合適的待修正參數，提高了代理模型的精度和修正效率。同時，選擇能夠度量兩個分布之間差距的KL 散度構造目標函數，以海上風機頻率均值和標準差為修正目標進行隨機模型修正，修正后的有限元模型響應與實際結構響應高度吻合。結果表明，基于Kriging 模型和KL 散度的隨機模型修正方法可以用于小樣本實測數據的海上風電結構隨機模型修正，且修正效果顯著。開展風電實驗模型修正將是未來進行的工作。