柳振海，李 剛，劉 博，徐明強(qiáng)，呂 晴，蔣 上

1.(中國(guó)電力工程顧問(wèn)集團(tuán)有限公司，北京 100120；2.中國(guó)海洋大學(xué)，山東 青島 266100 )

0 引言

海上風(fēng)電作為一種清潔能源，越來(lái)越受到世界各國(guó)的重視，大力發(fā)展海上風(fēng)電，可優(yōu)化能源結(jié)構(gòu)，促進(jìn)經(jīng)濟(jì)、社會(huì)、環(huán)境協(xié)調(diào)發(fā)展。然而，海上風(fēng)電設(shè)施長(zhǎng)期服役于復(fù)雜的海洋環(huán)境中，一旦發(fā)生故障或損傷，不僅影響風(fēng)電運(yùn)行，也會(huì)造成巨大的經(jīng)濟(jì)損失，開展結(jié)構(gòu)健康監(jiān)測(cè)勢(shì)在必行。獲取準(zhǔn)確的有限元模型是結(jié)構(gòu)健康監(jiān)測(cè)的重要內(nèi)容，但由于建模誤差、參數(shù)不確定性等問(wèn)題，有限元模型不一定能準(zhǔn)確反映結(jié)構(gòu)的真實(shí)情況，需要對(duì)有限元模型進(jìn)行修正，使得修正后的模型能準(zhǔn)確反映真實(shí)結(jié)構(gòu)的動(dòng)力響應(yīng)特征。

有限元模型修正方法可以分為確定性模型修正和不確定性模型修正。前者將結(jié)構(gòu)參數(shù)及響應(yīng)假設(shè)為恒定不變的，修正結(jié)果為唯一確定解。但在實(shí)際工程設(shè)計(jì)建設(shè)中，結(jié)構(gòu)部件的實(shí)際加工組裝可能與設(shè)計(jì)值存在一定誤差，導(dǎo)致結(jié)構(gòu)參數(shù)的不確定性。因此確定性模型修正得到的最優(yōu)解只是不確定性模型修正解中的一個(gè)特例。隨機(jī)模型修正是對(duì)傳統(tǒng)模型修正方法的創(chuàng)新，并考慮到了不確定性因素的影響，利用統(tǒng)計(jì)概率分析量化不確定性因素，將不確定模型修正問(wèn)題轉(zhuǎn)化為均值和標(biāo)準(zhǔn)差的修正問(wèn)題，從而得到與實(shí)際結(jié)構(gòu)動(dòng)力響應(yīng)特征一致的有限元模型[1-3]。

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

本文針對(duì)小樣本數(shù)據(jù)的模型修正提出了一種基于KL 散度的分步型隨機(jī)模型修正方法，通過(guò)數(shù)值模擬驗(yàn)證了此方法的有效性。首先采用靈敏度分析法選擇待修正參數(shù)，通過(guò)正交試驗(yàn)生成訓(xùn)練樣本，利用有限元模型計(jì)算樣本響應(yīng)，從而構(gòu)造以待修正參數(shù)為輸入，以頻率響應(yīng)為輸出的Kriging 代理模型；利用蒙特卡洛模擬得到初始樣本點(diǎn)并計(jì)算其響應(yīng)值，以樣本點(diǎn)響應(yīng)與實(shí)測(cè)數(shù)據(jù)的KL 散度構(gòu)造目標(biāo)函數(shù)，通過(guò)多目標(biāo)優(yōu)化算法對(duì)海上風(fēng)電結(jié)構(gòu)進(jìn)行隨機(jī)模型修正。

1 隨機(jī)模型修正方法

1.1 Kriging模型

相較于其他響應(yīng)面模型(如多項(xiàng)式回歸、神經(jīng)網(wǎng)絡(luò)等)，Kriging 模型在構(gòu)造過(guò)程中只需要較少的樣本點(diǎn)，計(jì)算成本較低，并且處理非線性數(shù)據(jù)效果優(yōu)異[8]，能實(shí)現(xiàn)對(duì)有限元模型響應(yīng)較準(zhǔn)確預(yù)測(cè)，更適合處理小樣本實(shí)測(cè)數(shù)據(jù)的模型修正。

Kriging 模型是一種通過(guò)已知樣本點(diǎn)信息估計(jì)未知試驗(yàn)點(diǎn)信息的無(wú)偏估計(jì)模型，采用Kriging 代理模型的前提是假設(shè)所有數(shù)據(jù)之間都服從n維正態(tài)分布[9]。

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

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

式中：R為數(shù)據(jù)點(diǎn)之間的協(xié)方差矩陣；σ2為方差。

樣本點(diǎn)的似然概率為：

式中：n為樣本數(shù)。

對(duì)式(3)取對(duì)數(shù)：

使得似然概率最大的超參數(shù)β，σ2的估計(jì)值和為：

對(duì)于建立的代理模型，還需要對(duì)其進(jìn)行進(jìn)一步的驗(yàn)證。本文采用決定系數(shù)R2和相對(duì)均方誤差RMSE評(píng)價(jià)Kriging 代理模型精度，R2的取值范圍為0 ～1，R2值越大，表示預(yù)測(cè)效果越好，RMSE越趨近于0 表示預(yù)測(cè)效果越好，其表達(dá)式為

1.2 KL散度

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

在優(yōu)化問(wèn)題中，需要通過(guò)迭代優(yōu)化樣本頻率Q使得兩個(gè)分布之間的KL 散度能夠盡可能小，從而實(shí)現(xiàn)用樣本頻率Q擬合實(shí)測(cè)頻率P，且實(shí)測(cè)頻率分布往往是多峰分布，前向KL 散度求解得到多個(gè)極小值對(duì)應(yīng)的待修正參數(shù)的平均值，不能滿足優(yōu)化要求，而采用反向KL 散度則可以確保至少達(dá)到一個(gè)局部最優(yōu)解，因此本文采用KL 散度計(jì)算實(shí)測(cè)頻率與樣本頻率分布之間的差異程度。

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

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

1.3 目標(biāo)函數(shù)

在優(yōu)化算法中，目標(biāo)函數(shù)提供了對(duì)不同解進(jìn)行評(píng)估和比較的指標(biāo)，決定了算法的搜索目標(biāo)，因此選擇適當(dāng)?shù)哪繕?biāo)函數(shù)對(duì)解決實(shí)際問(wèn)題至關(guān)重要。本文以待修正參數(shù)的均值和方差為修正目標(biāo)，以實(shí)測(cè)頻率與樣本點(diǎn)頻率的KL 散度為目標(biāo)函數(shù)，同時(shí)對(duì)參數(shù)的均值和方差進(jìn)行修正。由于兩個(gè)分布沒有重疊時(shí)，KL 散度沒有意義，無(wú)法進(jìn)行迭代操作，因此將目標(biāo)函數(shù)分為兩個(gè)部分：首先通過(guò)樣本點(diǎn)頻率及實(shí)測(cè)頻率的均值、方差構(gòu)造多目標(biāo)函數(shù)進(jìn)行優(yōu)化計(jì)算，當(dāng)兩個(gè)分布出現(xiàn)重疊時(shí)，改用KL 散度為目標(biāo)函數(shù)，具體目標(biāo)函數(shù)構(gòu)造如下：

式中：nsim為樣本響應(yīng)數(shù)量；nobs為實(shí)測(cè)值數(shù)量。

1.4 優(yōu)化算法

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

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

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

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

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

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

2 海上風(fēng)電結(jié)構(gòu)數(shù)值模型

選取一個(gè)5 MW 海上風(fēng)機(jī)作為研究對(duì)象，風(fēng)機(jī)采用導(dǎo)管架基礎(chǔ)，基礎(chǔ)深入海底，底部高程為-50.001 m，頂部高程為10.262 m，上部塔筒高度100 m，采用變截面，塔筒底部直徑6 m，壁厚0.035 1 m，頂部直徑3.87 m，壁厚0.024 7 m，頂部機(jī)箱和葉片簡(jiǎn)化為偏心質(zhì)量點(diǎn)，共31.452 t。利用ANSYS 建立有限元模型，模型共166 個(gè)節(jié)點(diǎn)，221 個(gè)單元。對(duì)其進(jìn)行模態(tài)分析，得到前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，前三階模態(tài)振型如圖1 所示。

圖1 海上風(fēng)機(jī)前三階模態(tài)振型

3 海上風(fēng)電結(jié)構(gòu)隨機(jī)有限元模型修正

3.1 修正參數(shù)的選擇

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

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

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

表1 修正參數(shù)各階頻率靈敏度

表2 修正參數(shù)初始均值與方差

3.2 建立代理模型

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

表3 正交試驗(yàn)設(shè)計(jì)樣本

采用相對(duì)均方誤差RMSE與決定系數(shù)R2對(duì)建立的Kriging 代理模型進(jìn)行有效性評(píng)價(jià)，計(jì)算結(jié)果見表4 所列，證明構(gòu)造的Kriging 代理模型精度較高，可以代替有限元模型。

表4 Kriging代理模型有效性評(píng)價(jià)

3.3 數(shù)值模擬

假設(shè)有限元模型待修正參數(shù)試驗(yàn)值與初始值見表5 所列，對(duì)試驗(yàn)值進(jìn)行蒙特卡洛抽樣1 000 次得到樣本點(diǎn)，通過(guò)Kriging 模型計(jì)算得到樣本點(diǎn)前兩階頻率響應(yīng)作為實(shí)測(cè)數(shù)據(jù)。首先采用實(shí)測(cè)數(shù)據(jù)與每次迭代生成的隨機(jī)樣本頻率的均值和方差構(gòu)造多目標(biāo)函數(shù)，隨后利用NSGA-Ⅱ多目標(biāo)優(yōu)化算法進(jìn)行求解，得到結(jié)構(gòu)待修正參數(shù)的修正值，見表5 所列。

表5 參數(shù)修正前后誤差

表5 中標(biāo)準(zhǔn)差是實(shí)測(cè)數(shù)據(jù)與其平均數(shù)離差平方的算術(shù)平均數(shù)的平方根，反映實(shí)測(cè)數(shù)據(jù)的離散程度，誤差是試驗(yàn)值與初始值之間的差異。根據(jù)修正參數(shù)最優(yōu)解，通過(guò)構(gòu)造的Kriging 模型進(jìn)行計(jì)算，得到修正后有限元模型前五階頻率的均值和標(biāo)準(zhǔn)差，與初始值和試驗(yàn)值進(jìn)行比較，見表6、表7 所列，修正前后頻率均值的最大誤差由68.59%降低為0.1%，標(biāo)準(zhǔn)差的最大誤差由212.82%降低為17.29%，且前三階頻率標(biāo)準(zhǔn)差最大誤差為0.48%，修正結(jié)果表明本文方法有效。

表6 修正前后頻率均值

繪制修正前后模型與試驗(yàn)?zāi)Ｐ皖l率置信橢圓如圖2 所示，可以看出，修正前實(shí)測(cè)值與初始值頻率置信橢圓大小與中心點(diǎn)位置有明顯差異，修正后實(shí)測(cè)值與修正值的頻率置信橢圓的大小、中心點(diǎn)、傾斜方向幾乎完全相同。圖3給出了結(jié)構(gòu)前兩階頻率修正前后的概率密度曲線，對(duì)比可知，修正前實(shí)測(cè)值曲線與初始值曲線的峰值有明顯差異，兩條曲線重合程度較低。修正后的實(shí)測(cè)值與修正值概率密度曲線形狀相似且?guī)缀跬耆呛?。進(jìn)一步驗(yàn)證了本文所提修正方法的有效性。

圖2 修正前后模型與試驗(yàn)?zāi)Ｐ皖l率置信橢圓

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

4 結(jié)語(yǔ)

考慮到參數(shù)隨機(jī)不確定對(duì)結(jié)構(gòu)響應(yīng)分布的影響以及KL 散度在不確定性度量的優(yōu)勢(shì)，本文采用Kriging 模型代替有限元模型進(jìn)行計(jì)算，并通過(guò)靈敏度分析法選擇合適的待修正參數(shù)，提高了代理模型的精度和修正效率。同時(shí)，選擇能夠度量?jī)蓚€(gè)分布之間差距的KL 散度構(gòu)造目標(biāo)函數(shù)，以海上風(fēng)機(jī)頻率均值和標(biāo)準(zhǔn)差為修正目標(biāo)進(jìn)行隨機(jī)模型修正，修正后的有限元模型響應(yīng)與實(shí)際結(jié)構(gòu)響應(yīng)高度吻合。結(jié)果表明，基于Kriging 模型和KL 散度的隨機(jī)模型修正方法可以用于小樣本實(shí)測(cè)數(shù)據(jù)的海上風(fēng)電結(jié)構(gòu)隨機(jī)模型修正，且修正效果顯著。開展風(fēng)電實(shí)驗(yàn)?zāi)Ｐ托拚龑⑹俏磥?lái)進(jìn)行的工作。