基于高斯過程回歸的橋梁多變量地震易損性分析

閆業祥， 孫利民,2,3

(1.同濟大學 橋梁工程系，上海 200092；2.同濟大學 土木工程防災國家重點實驗室，上海 200092；3.上海期智研究院，上海 200232)

地震易損性分析表達了在給定地震動強度的條件下工程結構達到或超過某個性能極限狀態的概率，已廣泛應用于土木基礎設施的震前風險預測、震后快速評估以及為抗震加固提供依據。在由太平洋地震工程研究中心發展的新一代基于性能的地震工程框架中，地震易損性分析是極其重要的一環，涵蓋了需求分析與損傷分析兩大模塊[1]。地震易損性分析的一個重要結果是地震易損性曲線，定義為結構在不同的地震強度指標(intensity measure, IM)下對各性能極限狀態的超越概率函數。但是，當前廣泛使用的地震易損性曲線建立在單一的參數(即一個標量IM，譬如峰值地面加速度PGA，或結構第一階周期處的譜加速度)之上，其不足主要體現在兩個方面：① 單一的IM不足以表征地震動對結構的破壞能力，而矢量IM可能更有效[2]；② 結構相關參數中存在的不確定性對易損性有不可小覷的影響[3]。因此，傳統的基于單參數的地震易損性曲線無法彌補上述不足，亟需一個新的工具來考慮多重變量對易損性的影響，這正是本文發展多變量地震易損性分析方法的背景與目的。

當前已有一些學者對多變量地震易損性分析作了有益的探索[4-6]，其中基于代理模型的Logistic回歸是一種較為有效的方法。該方法主要利用代理模型(也稱元模型或近似模型，譬如響應面模型、神經網絡回歸等)來代替有限元“黑箱”分析，通過代理模型預測的大樣本地震需求來構造生存-失效二值矢量(0代表安全，1代表超越了極限狀態)，并在此基礎上擬合Logistic回歸獲得構件的多變量地震易損性函數。上述方法的主要缺點在于流程上的繁瑣以及單一的IM仍被使用，而且易損性函數的精度極大地依賴于代理模型的預測精度以及Logistic回歸的擬合精度。不同于上述方法，本文提議了一個新的基于地震需求均值和需求方差模型的多變量地震易損性分析方法，并可以直接包含矢量IM與多重結構參數。其中均值與方差采用高斯過程回歸來擬合，高斯過程回歸是一個非參數機器學習方法，最近一些研究已將其應用于結構的抗震性能評估中[7-9]。作為多變量易損性函數的直接擴展，本文發展了均值與包絡易損性曲線的概念以充分考慮各種不確定參數的影響。最后，識別了不同的IM對易損性曲線離散性的影響，并確定了最優的IM形式。

1 基于高斯過程回歸的多變量地震易損性模型

1.1 多變量概率地震需求模型

當前的地震易損性評估方法主要采用非線性動力分析進行，即在一系列結構參數與地震動時程的隨機組合下進行彈塑性動力分析，然后基于結構地震需求(或損傷指標)和輸入參數進行后續的易損性分析。當采用多元回歸的方式構造地震需求與輸入參數之間的數學關系時，可建立如下表達式

lnyD=g(lnX)+ε

(1)

式(1)實質上是使用回歸模型來代理隱式的非線性有限元模型，它將地震需求分為了兩部分：均值部分與誤差部分。均值部分呈現了對平均地震需求的預測；而誤差部分則主要來源于如下兩個方面：① 有限的IM難以全面呈現地震動之間的變異性；② 回歸模型本身的擬合能力難以完全代替真實的有限元模型。因此，本文建議對均值與誤差分別建立回歸模型來預測結構的地震需求

mD=gm(X)

(2)

(3)

r2=(yD-gm(X))2

(4)

當然，一個常數方差估計也可利用下式得到

(5)

比較式(3)和式(5)不難發現，常數方差僅僅平均了殘差平方，而式(3)則是一個更加靈活的回歸模型，可以捕捉到可能存在的異方差模式，即方差隨著輸入參數的變化而出現一定的趨勢。

1.2 高斯過程回歸

式(1)提供了一個多變量需求模型的框架，其中回歸模型可以是任意的。本文使用高斯過程回歸作為特定的代理模型考慮到其強大的非線性擬合能力以及靈活的超參數求解技術。

所謂的高斯過程回歸即是假定地震需求在所有的N組樣本上服從一個均值矢量為mN和協方差矩陣為CN的高斯過程先驗

yD～GP(mN,CN)

(6)

在實際使用時，mN可以設置為零元素矢量(即響應數據減去其均值)，而協方差矩陣CN∈RN×N中的每個元素可以表達為

(7)

(8)

(9)

式中,κ*為核函數矢量，具有如下形式

(10)

(11)

結合式 (9)～(11) 可以看出在新樣本上的預測即為利用既有的N組地震響應樣本進行加權求和，距離新樣本點越近，則對預測的貢獻越大。

1.3 多變量地震易損性函數

(12)

式中：X*是由x*按行組成的矩陣；K*則是κ*按行構造的矩陣；下標m和β分別用來指代均值與方差的高斯過程回歸預測。式(12)之所以是多參數的，是因為新的樣本點中可以包括任意的IM和結構參數矢量。之所以稱作易損性函數，是因為IM可以設置為逐漸增加的形式。

式(12)可直接用于橋梁構件層次的易損性模型，此時每個構件對應一個地震需求矢量yD。為了獲得一個參數化的系統層次易損性函數，本文提議一個基于多元高斯Copula函數的方法。首先根據結構可靠度理論，多重構件的串聯系統的易損性函數[11]可以表示為

(13)

式中：m為構件的數目；Fi,Fj,…Fm分別為各構件的易損性函數；P(FiFj…Fm)為所有構件之間的聯合易損性函數，其余項類似。式(13)的關鍵點在于聯合易損性函數的求解，用高斯Copula函數可以近似的代替為

P(FiFj…Fm)=CGauss(F)=

ΦR(Φ-1(F1),…,Φ-1(Fm))

(14)

式中：CGauss(·)為高斯Copula函數；ΦR(·)為協方差矩陣等于其相關系數矩陣R的多元高斯累計分布函數；Φ-1(·)為單變量標準正態分布的逆累積分布函數。式(14)的前提假定是將各構件的易損性函數看作某個變量的累積分布函數，這與經典易損性模型中的假定是一致的[12]。將式(12)和式(14)代入式(13)，即可獲得系統層次的易損性函數，且保留了多變量和參數化的特性。應當注意本文僅分析了基于高斯Copula函數的串聯系統多元易損性模型，至于其他的系統形式(如并聯和混合串并聯等)，可直接擴展式(13)；對于其他Copula函數類型(如Clayton、Frank、Gumbel等形式)，可直接擴展式(14)。更多的發現將有賴于詳細的比較研究，限于篇幅本文暫未考慮。

1.4 均值與包絡易損性曲線

將任意給定的參數組合代入式(12)～(14)即可獲得一個特定極限狀態的超越概率，因此在連續變化的參數組合下其表現為高維空間中的易損性曲面。為了對應經典的單變量易損性曲線，可通過下式的高維積分得到

Fm(im)=FIM1[FIM(X)]=

(15)

式中：EIM1[·]為向標量IM1求取期望的算子；IM*為除去IM1以外的地震強度指標；fΘ和fIM*分別為結構參數Θ與IM*的聯合概率密度函數?？梢娛?15)為一個以積分求期望的過程(先積分結構參數再積分強度指標)，但是當變量維度很高或者變量間有相關性時變得十分復雜，甚至不可解。因此本文提議使用如下的非線性曲線擬合方式以獲得易損性曲線

(16)

圖1顯示了式(16)的擬合過程，首先將易損性投影到F-IM1平面內，然后通過非線性最小二乘擬合出一個參數為λim和βim的對數正態分布的累積分布函數。式(16)所求的易損性曲線實質上是考慮了所有參數隨機性的一個平均，因此將其稱作均值易損性曲線。

圖1 均值易損性曲線與包絡易損性曲線

既然投影到F-IM1平面內的易損性是離散的，且式(16)代表了其均值，那么其包絡易損性曲線也可以容易獲得。首先提取出離散易損性點的包絡點，對上界和下界包絡點分別應用式(16)即可獲得包絡易損性曲線，如圖1所示。

2 多變量地震易損性分析步驟

本文提議的多變量地震易損性分析的具體步驟,如圖2所示。

步驟1：選擇數量為N0的地面運動記錄，對于具有兩個正交水平分量的記錄可取其幾何均值作為新的記錄。確定用來縮放地面運動的地震強度指標IMs(譬如PGA或Sa(T1))以及縮放目標區間[IMsmin,IMsmax]。

圖2 多變量地震易損性分析流程圖

步驟2：選定隨機結構參數Θ。

步驟3：利用拉丁超立方抽樣(LHS)技術抽取行數為N的結構參數矩陣。注意將上述IMs設置成[IMsmin,IMsmax]內的均勻分布，和結構參數一起進行抽樣。并以生成的IMs樣本為目標對N0條地震波進行縮放。如N0

步驟4：利用生成的N個結構參數組合與N條縮放后的地震波代入到有限元模型進行非線性動力分析，得到每個構件的地震需求矢量yD。

步驟5：將yD和X進行對數變換后進行高斯過程回歸，分別建立均值需求模型和方差需求模型。

步驟7：假定系統為串聯，利用高斯Copula函數建立系統的多變量地震易損性模型Fsys(X*)。

步驟8：分別擬合構件與系統層次的均值易損性曲線與包絡易損性曲線。

需要強調的是步驟3中將縮放目標IMs與結構參數一起進行拉丁超立方抽樣的目的是使生成的結構參數樣本和目標IM保持空間填充特征，這類似于增量動力分析(IDA)，區別在于IDA的縮放IM是離散的，而本文是在區間[IMsmin,IMsmax]內均勻連續分布的。因此，這種均勻縮放的方法能夠以盡量少的動力分析次數獲得盡可能高的易損性精度，而這也得到了后文的蒙特卡羅方法的驗證。

3 橋梁結構與有限元模型

3.1 有限元模型與本構關系

本文的算例為一座60 m+100 m+60 m三跨鋼筋混凝土連續梁橋，采用單柱式橋墩，兩墩等高16 m，空心矩形截面，尺寸為縱橋向3 m，橫橋向8 m，壁厚0.8 m。上部結構為變截面連續箱梁，寬17 m，高3 m。支座為盆式橡膠支座。由于本文目的是發展易損性分析方法，為簡化問題忽略了基礎對橋梁動力響應的貢獻。利用地震工程仿真平臺OpenSees建立了如圖3所示的橋梁有限元模型。上部結構采用彈性梁單元模擬；橋墩采用基于力的纖維截面梁柱單元模擬，其中混凝土纖維的本構采用Concrete01材料模擬，核心區混凝土采用Mander模型來計算箍筋的約束效應，鋼筋纖維采用Steel01材料模擬。橋臺與填土之間的被動力-位移關系使用Hyperbolic Gap材料模型，橋臺與主梁之間的沖擊則利用Impact材料模擬。對于可滑動支座，使用理想彈塑性模型來模擬其摩擦作用。圖3(b)列舉了這些本構模型的參數，其具體解釋請參照OpenSees的用戶手冊，關于參數的計算方法可參見文獻[13]。利用該有限元模型，并將第3.2節中的隨機參數設置為其均值，可以計算出橋梁的動力特征。由計算得橋梁的前兩階模態周期分別為0.71 s和0.42 s，其振型均包含橋墩的縱向彎曲。因此為節省篇幅，本文僅分析縱橋向上的易損性，并利用這兩階振型的頻率計算Rayleigh阻尼矩陣的質量和剛度比例系數，用于隨后的動力計算。

(a) 橋梁有限元模型

3.2 隨機參數

既然式(12)所示的易損性函數可以包含多個變量，因此可以將橋梁參數設置為隨機變量以考慮其不確定性對易損性的影響。本文共選取21個與橋梁地震響應有關的隨機參數，如表1所示。這些參數包括4個無約束混凝土本構參數(fcc-峰值受壓強度，εcc-與峰值強度對應的應變，fcu-極限受壓強度，εcu-極限受壓應變)；1個與所有混凝土構件相關的參數ρc-混凝土質量密度；4個受約束混凝土的本構參數，其參數符號意義與無約束下的相同；3個與鋼筋本構相關的參數(fys-屈服應力，Es-彈性模量，bs-屈服后剛度比)；3個與支座摩擦有關的參數(Fy1-橋墩上支座的屈服力，Fy2-橋臺上支座的屈服力，δy-屈服位移)；3個與橋臺-主梁沖擊有關的參數(K1-初始剛度，K2-二次剛度，xy-屈服位移)；2個與橋臺-填土被動擠壓有關的參數(Kmax-初始剛度，Fult-極限被動土壓力)；1個橋梁阻尼比參數。

表1 各隨機結構參數的概率分布特征

3.3 構件的損傷指標與極限狀態

本文共對五個構件進行易損性分析，包括橋臺后回填土，橋臺上支座以及墩上支座，縱向固定支座下的橋墩墩底和墩頂截面。橋臺和回填土之間的峰值壓縮位移，支座的縱向峰值滑動位移以及橋墩截面的曲率Park-Ang指標分別作為其損傷測度。曲率形式的Park-Ang指標定義如下

(17)

式中：φy和φu分別為等效屈服和極限曲率，My為等效屈服彎矩，上述參數可由截面的彎矩-曲率分析得到；φm是地震產生的截面峰值曲率；ET是由滯回耗散的能量；因子βPA用來呈現循環荷載對結構損傷的影響，本文取為0.294。

共定義輕微、中等、嚴重和完全破壞四種極限狀態(limit states, LS)，每個構件的每種極限狀態所對應的能力概率分布如表2所示。所有構件的抗力被假定為對數正態分布，并為輕微和中等破壞采用較低的變異系數0.25，而為嚴重和完全破壞狀態指定更高的變異系數0.5，關于確定分布參數的細節參見文獻[12-14]。

表2 各構件極限狀態與能力的均值

3.4 地震動記錄

本文選取了報告FEMA P695[15]中整理的50條地震動記錄，包括22條遠場記錄，14條近場無脈沖記錄，以及14條近場含脈沖記錄。這50條地震波將根據第2章中的設計樣本總數N進行縮放。譬如N=500，則先將50條地震波復制10次，然后根據拉丁超立方抽樣生成的500個IMs樣本進行縮放，這樣既保證了地震波的隨機縮放特征，又保證了和結構參數一起的空間填充特征，因此共進行N次非線性動力分析。關于地震波的細節參見FEMA P695，本文不再列出。

對于本文，考慮到PGA更加容易縮放，將其作為縮放IMs并在[0,2g]內采樣。為了充分考慮到其余IM對易損性精度的貢獻，本文增加了另外3種地震動強度指標，分別為峰值地面加速度PGV，結構第一階自振周期處的譜加速度Sa(T1)，以及周期區間[0.2T1,2T1]內的幾何平均譜加速度Saavg。為了呈現縮放后的地震波的上述IM的分布特征，圖4顯示了N=500時各IM的分布，其中PGA服從[0,2g]內的均勻分布，這正是本文所采用的均勻縮放所期望的結果，而其余三個IM擬合成了對數正態分布。對于算例橋梁T1=0.71 s，Sa(T1)與Saavg的均值分別為1.57和1.19g，而PGV與結構無關，均值為1.44 m/s。

4 分析結果

4.1 地震需求均值與方差模型

按照第2章中的流程即可獲得各構件的地震需求矢量yD與輸入參數矩陣X，對數變換后應用高斯過程回歸即可獲得結構的地震需求均值與方差模型。圖5顯示了構件1和構件5的地震需求以及殘差平方與PGA的關系，可以看出即使經過了對數變換地震需求與PGA之間的關系仍然呈現出了非線性，意味著傳統的線性回歸是不適用的，而高斯過程回歸則可以較好地捕捉到非線性模式。此外，構件5的殘差明顯表現出了異方差模式，即隨著PGA的增加方差逐漸增大。若使用常數方差估計，顯然會在較高的PGA處低估其方差，而再次使用高斯過程回歸進行方差預測則可以靈活地適應方差的變化。另外需要強調的是構件1的異方差模式不顯著，則高斯過程回歸直接給出了常數方差估計，即如果從殘差中學習不到任何模式，其最終會給出一個均值估計。因此，高斯過程回歸可以靈活地適應任何方差模式，這也是本文發展需求方差模型的原因。其余構件與此類似，不再贅述。

4.2 樣本尺寸與模型精度驗證

實際上，樣本尺寸會影響高斯過程回歸的預測性能，進而影響到由此生成的多變量易損性函數的精度。為了研究其影響以及驗證所提方法的精度，使用蒙特卡洛仿真來直接生成各構件的易損性。具體做法為將PGA在區間內分成以0.1g為增量的20個離散的水平，在每一個PGA水平利用第2節中的抽樣技術抽取5 000個樣本進行非線性動力分析，各極限狀態下的超越概率可利用失效的樣本數直接除以總數(即5 000)得到。因此基于蒙特卡洛的易損性總共實施了10萬次動力分析。為了測量所提方法生成的均值易損性曲線與蒙特卡洛之間的誤差，定義如下的指標

(18)

式中：Fi為所提議方法生成的均值易損性曲線在每個PGA水平處的值，Fmc,i則對應蒙特卡洛計算的失效概率，因此指標Error測量了二者在所有PGA水平上的平均易損性偏差。

圖6為各構件與系統層次的易損性誤差隨樣本尺寸N的變化情況，在N=50時所有構件和系統有較大的誤差，這意味著較少數量的樣本無法為高斯過程回歸提供充足的訓練。但是當N=100時易損性誤差迅速減小，在N=300趨于穩定，即不再隨著N的不斷增加而減小。事實上，這些不能被縮減的誤差是模型誤差，來源于高斯過程回歸的預測，以及在生成易損性函數時的正態性假定。若想繼續減小這些誤差，需要發展更加合理的代理模型以及易損性生成方法。

圖6 易損性精度與樣本尺寸的關系

對于本文算例來說，N=300可以認為是一個最優的樣本尺寸，即在保證精度的同時最大程度地降低有限元分析的代價。圖7顯示了N=300時所有構件與極限狀態的誤差，這些誤差都在0.02以下，這對于地震作用下的高失效概率來說是可接受的。此外，觀察到系統層次的易損性誤差也是較低的，表明本文所提議的基于多元高斯Copula函數的系統易損性函數生成方法是可靠的。至于本文所提議方法的效率，由于蒙特卡洛法的總分析次數為10萬，因此效率比為300/(1×105)=0.003，可見本文所提議的方法具有非常高的計算效率。

圖7 N=300時各構件與系統的易損性誤差分布

4.3 IM的影響

4.2節驗證了所提議方法的精度與效率，即生成的均值易損性曲線收斂于蒙特卡洛的解，不過使用的是PGA，使用PGA的原因是其作為縮放IM。但是其他的IM也可以用來生成均值易損性曲線，圖8顯示了由不同IM生成的均值易損性曲線的結果，每個IM的區間由圖4推斷得到?？梢钥闯龈鲌D最大的區別在于易損性投影點的離散性，其中Saavg的最小，其次Sa(T1)和PGV，離散性最大的是PGA。越高的離散性說明還有更重要的參數影響著易損性，而其本身的影響則相對較小。因此，Saavg是一個更好的測度用來生成均值易損性曲線。盡管這樣，Saavg仍具有一定的離散性，說明其它參數(包括其余IM以及結構參數)的影響也是不可忽略的，這也意味著任何單個參數都不足以完全呈現易損性，佐證了本文發展多變量易損性函數的必要性。此外，需要強調的是盡管PGA的離散程度較大，但是上節證明了其均值易損性曲線是精確的，因為其平均了所有不確定參數的影響，而離散性反映的是其不確定程度，可以由包絡易損性曲線來呈現。

4.4 均值和包絡易損性曲線

根據4.3節的分析，圖9顯示了所有構件以及系統層次以Saavg生成的均值與包絡易損性曲線。由圖可知幾乎所有構件的所有極限狀態的均值易損性最終都達到了1，而且包絡易損性曲線的面積也幾乎縮減至0，意味著失效是必然的，沒有任何不確定性存在。這里Saavg的區間設置為[0,4g]，這由圖4中的均值和方差推斷而來，即μ+3σ=1.19+3×0.88=3.83g，這說明雖然用來縮放的PGA限制在[0,2g]內，但是Saavg卻可以高達4g左右，對于本文算例來說是很高的地震強度。此外，上下包絡易損性之間的距離衡量了該極限狀態易損性的不確定性程度，可以發現這種不確定性呈現了兩端小(端部為0)，中間大的特點，而且不同的極限狀態之間會出現重疊現象，暗示著易損性曲線表征破壞狀態之間的模糊現象，需要未來更深入地探索。

表3列出了所有均值與包絡易損性曲線的中位值，即當超越概率為0.5時所需的Saavg強度，其值越小表明更加易損。表3顯示橋臺填土和支座比橋墩更加容易遭受地震損傷；由于固定支座下橋墩與上部結構表現為剛架形式，故而上端截面也承受較高的損傷概率，但稍低于底部截面；由于本文假定橋梁系統為串聯形式，所以其中位值Saavg低于所有構件。

表3 均值與包絡易損性曲線的中位值

最后需要強調雖然Saavg被建議用來擬合均值與包絡易損性曲線，但是其依賴于結構的自振周期，因此PGA等與結構無關的指標仍然被建議用來縮放地震波，考慮到縮放與生成均值、包絡易損性曲線是相互獨立的環節。

5 結 論

本文提議了一個基于高斯過程回歸的多變量地震易損性分析框架，可以將易損性表達為由地震動強度指標和結構參數組成的矢量輸入的多元函數。生成的均值與包絡易損性曲線能夠全面考慮地震波與結構參數中的不確定性。以一座三跨鋼筋混凝土連續箱梁橋為算例，結果表明本文所提議的方法與直接蒙特卡洛仿真具有相一致的精度，并且具有極高的效率，大大減少了有限元分析數量。通過比較研究發現，平均譜加速度Saavg是一個生成均值與包絡易損性曲線更加合適的指標，考慮到其具有最小的易損性離差，而且在區間端部不確定性縮減為零。此外，研究結果顯示任何單個參數均不能完全呈現結構的易損性，需要在分析中包含多重變量的影響。