丁曦，王佐才,2

（1.合肥工業(yè)大學(xué)土木與水利工程學(xué)院，安徽 合肥 230000 2.土木工程防災(zāi)減災(zāi)安徽省工程技術(shù)研究中心，安徽 合肥 230000）

0 前言

近些年，有限元分析得到廣泛應(yīng)用，但是由于受到有限元模型的網(wǎng)格劃分、邊界條件和材料物理參數(shù)等不確定性因素的影響，建立的有限元模型與真實(shí)結(jié)構(gòu)始終存在一定的差異，因此需通過試驗(yàn)數(shù)據(jù)加以修正，使其盡可能接近實(shí)際結(jié)構(gòu)。

關(guān)于有限元模型修正的方法可以分為三類，即直接類修正算法[1]、迭代類修正算法[2]和基于智能優(yōu)化技術(shù)的修正算法[3]。直接類修正算法是從結(jié)構(gòu)振動(dòng)微分方程出發(fā)，通過直接修改結(jié)構(gòu)剛度與質(zhì)量矩陣中的元素達(dá)到模型修正的目的，但由于修正后的參數(shù)物理意義不明確，此類方法已不常用。迭代類修正算法目前主要運(yùn)用在非線性結(jié)構(gòu)問題中，在實(shí)際選擇參數(shù)方面會(huì)更加的靈活，但是在實(shí)際運(yùn)用過程中存在計(jì)算量巨大的問題無法解決。

在實(shí)際運(yùn)用有限元模型計(jì)算的過程中，由于直接調(diào)用有限元軟件進(jìn)行模型修正存在著計(jì)算效率低、不方便二次開發(fā)等缺點(diǎn)，學(xué)者提出了采用代理模型并結(jié)合智能算法進(jìn)行模型修正的方法，此類方法可以克服修正過程中大量調(diào)用有限元模型計(jì)算繁瑣的不足。目前常用的代理模型包括響應(yīng)面模型[4]、徑向基函數(shù)模型與Kriging代理模型等。

其中，Kriging代理模型方法[5]也稱空間局部估計(jì)法，由南非地質(zhì)學(xué)者Krige于1951年提出，最初是用來確定礦產(chǎn)儲(chǔ)量分布。Sacks[6]等第一次將Kriging方法應(yīng)用于結(jié)構(gòu)設(shè)計(jì)。Rom ero[7]等將其應(yīng)用于求解結(jié)構(gòu)可靠度問題中。Kriging代理模型是一種半?yún)?shù)化的插值模擬建立模型的方法，考慮設(shè)計(jì)變量在空間上的相關(guān)特性，利用某一點(diǎn)周圍信息模擬某一未知點(diǎn)的響應(yīng)等信息[8]。

本文以濱湖路鋼桁架橋?yàn)檠芯繉?duì)象，選擇橋梁結(jié)構(gòu)的參數(shù)作為修正對(duì)象，基于Kriging代理模型和橋梁荷載試驗(yàn)數(shù)據(jù)，以修正橋梁的有限元模型，使有限元模型更接近橋梁的實(shí)際情況。

1 Krig ing模型修正理論

1.1 Kriging模型

Kriging模型方法作為線性回歸分析的改進(jìn)技術(shù)，由線形回歸模型與隨機(jī)過程相加組合而成，其基本形式為：

式中：x是維度為d的設(shè)計(jì)參數(shù)[x1，x2，…，xd]；y（x）為多項(xiàng)式模型；z（x）為正態(tài)分布函數(shù)，其不獨(dú)立但同分布；f（x）為類似響應(yīng)面模型，可以根據(jù)需要選擇零階、一階及二階多項(xiàng)式，m是多項(xiàng)式的數(shù)目。

f（x）負(fù)責(zé)模擬全局近似計(jì)算，z（x）負(fù)責(zé)模擬局部近似計(jì)算。z（x）是隨機(jī)過程，其均值為零且協(xié)方差不為零。在Kriging代理模型中，由于z（x）這一隨機(jī)過程的靈活運(yùn)用，其建立的代理模型具有很強(qiáng)的靈活性，這樣一來，模型可以更好的運(yùn)用于求解復(fù)雜情況的非線性結(jié)構(gòu)，解決了傳統(tǒng)基于多項(xiàng)式函數(shù)所建立的代理模型在精度上的不足。

z（x）的協(xié)方差矩陣具有如下形式：

式中：

R[R（ xi，xj）]為ns個(gè)樣本點(diǎn)中，任意兩個(gè)點(diǎn)xi、xj間的空間相關(guān)函數(shù)，其對(duì)模型的精度有主導(dǎo)作用，其核函數(shù)具有多種形式可以選擇，在實(shí)際運(yùn)用過程中，高斯函數(shù)通?？梢缘玫捷^好的計(jì)算結(jié)果，其形式如下：

1.2 Kriging模型參數(shù)的確定

對(duì)于參數(shù)θk的估計(jì)，通常用已知結(jié)構(gòu)響應(yīng)的線性組合來估計(jì)，引入以下矩陣：

由上式可得出結(jié)構(gòu)響應(yīng)估計(jì)值的殘值：

若上式所示殘差的平均值為0，則得到的結(jié)構(gòu)響應(yīng)值則為無偏估計(jì)值，可得到下式所示關(guān)系：

可以求得結(jié)構(gòu)響應(yīng)估計(jì)值的方差為：

式中：r為相關(guān)矩陣R中的某一列。通過求解結(jié)構(gòu)響應(yīng)估計(jì)值方差的最小值，并使其滿足無偏估計(jì)的約束，可建立優(yōu)化目標(biāo)函數(shù)：

式中：λ為拉格朗日乘子向量，通過優(yōu)化求解可得以下結(jié)果：

將上述結(jié)果引入結(jié)構(gòu)響應(yīng)估計(jì)值可得：

式中：空間相關(guān)矩陣R為：

上述兩式β*與σ2均是參數(shù)θk的函數(shù)，通過極大似然估計(jì)使擬函數(shù)取得極大值為：

求解關(guān)于以上優(yōu)化問題，可得到參數(shù)θk的估計(jì)值，即可建立Kriging模型。

基于Kriging模型的有限元模型修正方法也屬于響應(yīng)面方法的一種，通過不同函數(shù)建立的響應(yīng)面模型其適用性也不盡相同，所以，就需要對(duì)所建立的響應(yīng)面模型進(jìn)行精度驗(yàn)算。驗(yàn)證響應(yīng)面模型精度的方法有很多，如：殘差的正態(tài)分布檢驗(yàn)、殘差的均值、EISE檢驗(yàn)、R2檢驗(yàn)[10]及相對(duì)均方根誤差(Root Mean Squared Error)[9]檢驗(yàn)方法。在較為復(fù)雜的模型中，R2檢驗(yàn)和相對(duì)均方根誤差(RMSE)檢驗(yàn)可更好的檢驗(yàn)?zāi)Ｐ途?。其?shù)學(xué)表達(dá)式如下式所示：

相對(duì)均方根誤差(RMSE)：

判定系數(shù)：

式中：

y和yreg分別是設(shè)計(jì)空間上各點(diǎn)的真值和響應(yīng)面模型的值為設(shè)計(jì)空間上各點(diǎn)真值的均值；N為設(shè)計(jì)空間上檢驗(yàn)點(diǎn)的數(shù)量。相對(duì)均方根誤差(RMSE)的大小R2代表了響應(yīng)面的精度，可根據(jù)實(shí)際情況確定，其值越接近0則表示響應(yīng)面模型越精確。代表了響應(yīng)面與真值之間的差值，其值在0～1之間，當(dāng)其值為1時(shí)，表示二者完全一致。

2 初始有限元模型建立

2.1 橋梁概況

濱湖路橋位于宛溪河上，其主體結(jié)構(gòu)為兩跨鋼桁架橋其跨徑為50+50m，主桁架采用雙片桁架，邊支點(diǎn)桁高為4.5m，中支點(diǎn)桁高為 17.5m。上弦桿線形為二次拋物線形式，主桁架下弦桿采用焊接整體節(jié)點(diǎn)結(jié)構(gòu)形式，主橋橋梁橫斷面共33m，全橋材料為Q345qD鋼材。

橋型布置如圖1所示。

圖1 濱湖路橋橋型布置圖

2.2初始有限元模型的建立及計(jì)算結(jié)果

按照設(shè)計(jì)圖紙，使用有限元計(jì)算軟件M idas建立全橋有限元模型，上部結(jié)構(gòu)采用Q345qD鋼材，設(shè)計(jì)強(qiáng)度200MPa，容重取，考慮到加勁肋、橫隔板、橫肋的重量，橋梁結(jié)構(gòu)的自重系數(shù)取1.2。每側(cè)護(hù)欄取，每側(cè)燈柱取，橋面鋪裝取，人行道自重簡(jiǎn)化成荷載：，人群荷載集度?。?，汽車荷載取公路-Ⅰ級(jí)。

橋面系模型采用4片主梁，主桁架采用雙桁架模型進(jìn)行分析，全橋模型均采用梁?jiǎn)卧M(jìn)行模擬。全橋共離散成1112個(gè)單元。其有限元計(jì)算模型見圖2所示。

圖2 全橋有限元計(jì)算模型

濱湖路橋理論計(jì)算自振頻率與振型 表1

表1給出利用有限元模型計(jì)算得到的結(jié)構(gòu)前三階豎彎自振頻率與振型。

3 現(xiàn)場(chǎng)脈動(dòng)試驗(yàn)及結(jié)果

3.1 脈動(dòng)試驗(yàn)內(nèi)容

脈動(dòng)試驗(yàn)主要測(cè)量主橋的自振頻率、振型和阻尼比。試驗(yàn)選用壓電式加速度傳感器，記錄橋梁結(jié)構(gòu)在環(huán)境激勵(lì)下，如風(fēng)、水流、地脈動(dòng)等引起的橋梁振動(dòng)，然后對(duì)記錄下來的橋梁振動(dòng)時(shí)程信號(hào)進(jìn)行處理,并進(jìn)行時(shí)域和頻域分析，求出橋梁結(jié)構(gòu)自振特性。全橋共設(shè)置6個(gè)加速度傳感器測(cè)點(diǎn)，試驗(yàn)中加速度傳感器測(cè)點(diǎn)布置如圖3所示。

對(duì)環(huán)境激勵(lì)下橋梁的響應(yīng)信號(hào)進(jìn)行多次功率譜的平均分析，可得到橋梁的各階自振頻率，再利用各個(gè)測(cè)點(diǎn)的振幅和相位關(guān)系，可求得橋梁各階模態(tài)相應(yīng)的振型，利用幅頻圖上各峰值處的半功率帶寬或時(shí)域上的自相關(guān)確定各階模態(tài)阻尼比。利用脈動(dòng)試驗(yàn)測(cè)出全橋結(jié)構(gòu)前三階豎彎振型以及頻率。

圖3 全橋脈動(dòng)試驗(yàn)測(cè)點(diǎn)布置圖

3.2 脈動(dòng)試驗(yàn)結(jié)果分析

圖4 1號(hào)測(cè)點(diǎn)加速度時(shí)程響應(yīng)

圖5 基于SSI的豎向振動(dòng)頻率穩(wěn)定圖

濱湖路橋自振頻率與振型實(shí)測(cè)結(jié)果 表2

相關(guān)參數(shù)的初始值 表3

本文中，使用隨機(jī)子空間方法(SSI)方法，處理所采集的響應(yīng)信號(hào)，識(shí)別結(jié)構(gòu)的模態(tài)參數(shù)。圖4為1號(hào)測(cè)點(diǎn)的加速度時(shí)程響應(yīng)。圖5為基于SSI的豎向振動(dòng)頻率穩(wěn)定圖。表2列出了該橋?qū)崪y(cè)自振頻率以及振型圖。

4 有限元模型修正

4.1 Kriging代理模型建立的試驗(yàn)設(shè)計(jì)

利用Kriging代理模型進(jìn)行有限元模型修正，即通過隱式函數(shù)關(guān)系來表示隨機(jī)輸入變量與輸出變量之間的關(guān)系，并采用相關(guān)試驗(yàn)來擬合假定的函數(shù)關(guān)系，利用其替代原復(fù)雜結(jié)構(gòu)的有限元求解過程?；贙riging代理模型來進(jìn)行有限元模型修正通過四個(gè)步驟來實(shí)現(xiàn)：①試驗(yàn)設(shè)計(jì)，設(shè)計(jì)一定的試驗(yàn)方案，確定結(jié)構(gòu)的輸入與輸出信息；②確定Kriging模型的回歸模型與相關(guān)模型函數(shù)的具體形式；③建立代理模型，通過相應(yīng)的數(shù)據(jù)擬合出響應(yīng)面模型，并檢驗(yàn)?zāi)Ｐ偷臄M合精度；④優(yōu)化與預(yù)測(cè)，利用實(shí)測(cè)響應(yīng)值與參數(shù)通過Kriging模型進(jìn)行分析與預(yù)測(cè)，得到優(yōu)化后的結(jié)構(gòu)參數(shù)以及其對(duì)應(yīng)的有限元模型。

基于上述理論，首先對(duì)橋梁脈動(dòng)試驗(yàn)進(jìn)行分析，選定濱湖路橋結(jié)構(gòu)相關(guān)參數(shù)作為修正影響因素分析。本文，選取結(jié)構(gòu)的彈性模量、泊松比、線膨脹系數(shù)與容重作為需要修正的影響因素。選取結(jié)構(gòu)的各參數(shù)的初始值如下表所示，在此基礎(chǔ)上，選擇相關(guān)參數(shù)變化空間為初始參數(shù)的0.75～1.25倍。表3為相關(guān)參數(shù)的初始值。

本文基于均值設(shè)計(jì)，樣本設(shè)定為U12（124），即4參數(shù)，12組樣本。表4和表5為12水平均勻設(shè)計(jì)表與12水平均勻設(shè)計(jì)表的使用表，均勻設(shè)計(jì)法可以很好的將參數(shù)分布在整個(gè)設(shè)計(jì)參數(shù)空間內(nèi)，并且可以大大地減少試驗(yàn)次數(shù)。

根據(jù)表6與表7均勻設(shè)計(jì)表及其使用表，設(shè)計(jì)如下12組Kriging代理模型試驗(yàn)。

根據(jù)12水平均勻設(shè)計(jì)表及其使用表格，設(shè)計(jì)各組試驗(yàn)，一共12組，通過改變?cè)O(shè)計(jì)的參數(shù)，利用建立的有限元計(jì)算模型，得到所需頻率值，其結(jié)果如表7所示。

12水平均勻設(shè)計(jì)表 表4

12水平均勻設(shè)計(jì)表使用表 表5

Kriging代理模型試驗(yàn)設(shè)計(jì) 表6

4.2 代理模型的建立

通過上述均勻設(shè)計(jì)的樣本點(diǎn)作為已知信息輸入矩陣；將有限元計(jì)算軟件得到的相應(yīng)頻率值作為輸出矩陣。由于實(shí)際工程問題中非線性問題程度較高，本文中回歸模型選用線性函數(shù)；對(duì)于相關(guān)函數(shù)方程的選擇，由于高斯函數(shù)相對(duì)于其他幾種相關(guān)函數(shù)優(yōu)化模型的優(yōu)化效果最好，本文中相關(guān)函數(shù)選用高斯函數(shù)。

運(yùn)用M atlab軟件中的DACE工具箱構(gòu)建Kriging代理模型，本文中，采用粒子群優(yōu)化算法對(duì)參數(shù)進(jìn)行優(yōu)化，用確定的參數(shù)確定最優(yōu)相關(guān)模型。

基于Kriging理論所建立的代理模型為一個(gè)隱式函數(shù)，在本文中不給出具體的函數(shù)表達(dá)式，僅以響應(yīng)面圖形表示所建立的Kriging代理模型。圖6～圖8分別表示各主要結(jié)構(gòu)參數(shù)與前三階固有頻率值之間的Kriging模型。圖(a)為彈性模量和泊松比分別與前三階固有頻率之間的Kriging模型，圖(b)為線膨脹系數(shù)和容重與前三階固有頻率之間的Kriging模型。

將前三階實(shí)測(cè)頻率輸入構(gòu)建的Kriging代理模型中用于回歸預(yù)測(cè)，便可得到預(yù)測(cè)值，將預(yù)測(cè)所得的設(shè)計(jì)參數(shù)代入原有限元模型中以達(dá)到修正此有限元模型的目的。通過求解修正后的有限元模型便可得到修正后的理論振動(dòng)頻率。預(yù)測(cè)參數(shù)值如下表8所示，表9給出了通過Kriging代理模型修正后的理論振動(dòng)頻率與實(shí)測(cè)振動(dòng)頻率以及修正前的原始有限元模型計(jì)算所得振動(dòng)頻率。

由表9可知，經(jīng)Kriging代理模型修正后計(jì)算所得結(jié)構(gòu)頻率與真實(shí)值相差均在以內(nèi)。

4.3 Kriging代理模型精度檢驗(yàn)

驗(yàn)證Kriging代理模型所得預(yù)測(cè)值的精度，本文引入前文所述的檢驗(yàn)及相對(duì)均方根誤差(RMSE)對(duì)修正后的模型進(jìn)行評(píng)估。表10表示了所建立的Kriging代理模型的精度計(jì)算值。

圖6 各主要參數(shù)與第一階固有頻率Kriging模型

圖7 各主要參數(shù)與第二階固有頻率Kriging模型

圖8 各主要參數(shù)與第三階豎彎頻率Kriging模型

設(shè)計(jì)試驗(yàn)所得頻率值 表7

經(jīng)Kriging代理模型修正后的參數(shù)值 表8

實(shí)測(cè)頻率與修正前后頻率結(jié)果對(duì)比表 表9

各階次Kriging代理模型精度分析 表10

從上表可知，相對(duì)均方根誤差(RMSE)值均很接近0，這表明所建立的Kriging代理模型與有限元模型計(jì)算值之間的差異很??；在值計(jì)算中，各階次的Kriging代理模型精度的計(jì)算結(jié)果均在0.95以上，這表明建立的Kriging代理模型可以準(zhǔn)確地描述系統(tǒng)輸入與輸出之間的關(guān)系。此時(shí)建立的Kriging代理模型已經(jīng)可以很好的進(jìn)行有限元參數(shù)優(yōu)化。

5 結(jié)論

本文介紹了Kriging代理模型在有限元模型修正中的應(yīng)用，并通過濱湖路橋?qū)嵗f明Kriging代理模型結(jié)合有限元軟件在橋梁工程中的實(shí)際應(yīng)用。

①舉例說明Kriging代理模型在橋梁工程中的具體應(yīng)用，利用實(shí)際測(cè)量的橋梁結(jié)構(gòu)動(dòng)力特性去修正橋梁實(shí)際的結(jié)構(gòu)設(shè)計(jì)參數(shù)，并驗(yàn)證了修正后參數(shù)的具體物理意義；

②驗(yàn)證了利用橋梁荷載試驗(yàn)，建立Kriging代理模型優(yōu)化原有限元模型方法的可行性，這種方法可以減小橋梁試驗(yàn)對(duì)原橋梁結(jié)構(gòu)的影響；

③通過對(duì)既有橋梁進(jìn)行預(yù)測(cè)并與實(shí)測(cè)結(jié)果對(duì)比，結(jié)果表明，修正后的有限元模型可以更好的反映橋梁結(jié)構(gòu)的真實(shí)情況。