多年凍土場地樁基橋梁地震響應不確定性量化分析

李發達，蘇 雷，萬華平，凌賢長

(1. 青島理工大學 土木工程學院，山東 青島 266520； 2. 浙江大學 建筑工程學院，杭州 310058)

伴隨著國家西部大開發戰略，凍土區規劃和建設的橋梁基礎工程被逐漸提上日程；尤其在凍土區建設的樁基橋梁工程研究越來越多[1-3]。由于我國地處世界兩大地震帶之中，在大陸板塊相互作用的影響下，凍土區的交通基礎設施難免受地震影響，對其進行抗震方面的研究顯得尤為重要。在凍土區交通基礎設施的地震響應特性明顯受輸入參數的影響，這些輸入參數往往存在一定的變異性，分析這些輸入參數變異性及其對交通基礎設施地震響應特性不確定性具有極其重要的意義。在凍土區交通基礎設施的地震響應分析中，定量分析結構體系輸入參數對其整個結構地震響應的影響是一個重要的研究方向[4-5]；定量分析結構體系輸入參數變異性對結構地震響應的影響，即進行結構體系地震響應的不確定性量化分析將是另一個重要研究方向，備受研究者的關注。

目前，在交通基礎設施的抗震方面許多學者進行了廣泛的研究。Soleimani等[6]針對結構不規則和幾何不對稱橋梁的地震響應，采用人工統計學的方法進行相應敏感性分析，研究相關不確定性參數對概率地震需求模型的影響。萬華平等[7]為量化固有頻率的不確定性，使用高斯過程替代模型，計算結構形式復雜、長跨度橋梁結構參數的均值和方差，該方法不局限于橋梁結構參數的概率分布情形，可揭示橋梁結構動力響應特性。王德斌等[8]利用有限元軟件建立近海橋梁引橋段模型，結合一維波動理論與流體動力學方程模擬海底地震動，然后分析海底地震動對近海橋梁易損性的影響。

考慮到大型的樁承結構的非線性地震響應有限元分析耗時長且費用高等特點。Su等[9]在傳統的地震模擬方法的基礎上利用高斯過程替代模型進行大型樁承碼頭結構的地震響應不確定性量化。莫鋈[10]針對我國港口建設中常用的高樁碼頭在其設計年限內極易遭受地震破壞的問題，通過使用有限元法建立了二維高樁碼頭體系數值模型，在考慮樁-土相互作用下，對模型進行了反應譜和隨機地震響應分析。陳清軍等[11]為了探討樁-土-結構相互作用，創新性的以基巖隨機地震作為輸入，采用隨機振動理論與有限元模擬相結合的方法進行了三維隨機地震反應分析。值得注意的是，相關學者在凍土場地橋梁樁基的三維有限元模擬和不確定性量化分析方面研究較少。

一般而言，為了考慮突發事件(如地震、撞擊)對結構造成破壞，針對突發事件對結構造成的破壞進行分析顯得十分必要。特別是現實中對結構體系進行抗震性能分析和預測尤為重要[12-14]。周愛紅[15]針對地震對樁基礎的破壞十分普遍的現象，考慮地震時樁-土-結構之間的相互作用，進行隨機地震響應和可靠度分析。萬華平等[16]針對采用全局敏感性進行橋梁結構動力特性分析時具有計算花費高、耗時長等問題，提出了可以有效、快速的進行橋梁結構動力特性不確定性的敏感性分析方法。

三維有限元數值計算中，高斯過程替代模型具有明顯的優勢，最突出的地方是對有限元計算中衡量構筑物抗震性能的復雜高維積分問題可以簡化為簡單的一維積分。通過這種方式可在復雜的非線性數值模型計算中節省計算成本，而且計算相對復雜的非線性數值模型最明顯的特征是用時少。本文通過高斯過程替代模型在考慮輸入參數變異的影響下對多年凍土場地樁基橋梁地震響應進行不確定性量化分析。

1 基于高斯過程替代模型不確定性量化分析

1.1 基于Sobol序列的采樣方法

為了采用高斯過程替代模型進行地震響應不確定性量化分析，首先基于試驗設計方法構建高斯過程替代模型的訓練樣本，以此生成訓練樣本充滿參數空間的最大數量采樣點。通常，針對三維有限元數值計算中模型的節點多、機時長和模型復雜度高等特點，本文采用高斯過程替代模型進行多年凍土場地樁基橋梁的地震響應不確定性量化分析。一般而言，空間采樣方法主要包括：LHS抽樣、Hammersley序列、Halton序列和Sobol序列。為了最大限度的滿足訓練樣本能夠充滿參數空間的程度，這里使用均勻性良好的Sobol序列進行空間采樣。

Sobol序列屬于低偏差序列，是以一組連續素數為基，最小素數為2的一種空間采樣方法。以最小素數2為基在進行計算時可以避免序列生成高維序列，以此提高模型的計算效率和節省計算時間[17]。

基于一組直接數而構造低差異序列的Sobol序列的表達式為

n=(nmnm-1nm-2…n3n2n1)2=
n1+n221+n322+…+
nm-22m-3+nm-12m-2+nm2m-1

(1)

式中：m=[lnn/lnR]，[·]表示取整運算；nm,…,n1取0或1。

在生成Sobol序列時需要借助不可約的本原多項式

pi(x)=xSi+a1,ixSi-1+a2,ixSi-2+…+aSi-1,ix+1

(2)

式中：pi(x)為第i維的本原多項式；Si為本原多項式pi(x)的度數；a1,i,a2,i,a3,i,…,aSi-1,i取整數0或1。

下面部分通過遞推關系定義以下正數m1,i,m2,i,…,mSi,i

mk,i=2a1,imk-1,i⊕22a2,imk-2,i⊕…⊕
2Si-1aSi-1,imk-Si+1,i⊕2SiaSi,imk-Si,i⊕mk-Si,i

(3)

式中：⊕為二進制按位異或運算，即1 ⊕ 0=0，0 ⊕ 1=0，0 ⊕ 0=1，1 ⊕ 1=1；mk,i(k=1,2,…,Si)為小于2i的奇數。

直接數通過式(4)得到

(4)

通過式(4)可得到以最小素數為基的Sobol序列

xn,i=n0v1,i⊕n1v2,i⊕…⊕nmvm+1,i

(5)

對于基于最小素數2為基的Sobol序列，許多學者提出了不同的生成方法[18-20]。其中，Antonov等[21]提出了基于格雷碼快速生成Sobol序列的方法；針對格雷碼特有的性質，可由式(6)生成Sobol序列的遞歸公式

xk+1,i=xk,i⊕vck,i

(6)

式中，ck為整數k二進制最右零位。

1.2 高斯過程替代模型

通常情況下，三維非線性有限元數值模型節點多、計算量大，對計算機的硬件要求較高且計算成本大，這里采用高斯過程替代模型對多年凍土場地樁基橋梁的地震響應進行不確定性量化分析。高斯過程替代模型通過建立體系參數與地震響應之間的關系，將高維積分簡化為低維積分，在此基礎上進行不確定性量化分析，以此可實現多年凍土場地樁基橋梁的不確定性量化分析。

基于貝葉斯架構的高斯過程替代模型是一種非參數模型，其不同的均值函數和協方差函數可解決不同的研究問題[22]。其中，均值函數M(x)和協方差函數c(x,x′)表征的高斯過程可通過式(7)表示

s(x)～N[M(x),c(x,x′)]

(7)

式中：x={x1,x2,x3,…,xd}為d維參數向量；s為表征多年凍土場地樁基橋梁的地震響應；M(·)為均值函數；c(·,·)為協方差函數。

在高斯過程替代模型中，均值函數一般設為零；協方差函數可通過式(8)表示

(8)

式中：η2為信號變異性；lk為特征長度。

替代模型s(x)的構建需要進行多年凍土場地樁基橋梁的非線性有限元地震響應分析，以此得到所需的訓練樣本集。訓練樣本集用D=[X,S]表示，包含n對矢量參數樣本和標量輸出；其中，X={x1,x2,x3,…,xn}是參數指標，S={s1,s2,…,sn}是多年凍土場地樁基橋梁的地震響應向量。一般而言，使用高斯過程模型中需要預測模型中未知點x*的地震響應s*，即高斯過程替代模型s(x)的特征；使用高斯過程先驗優于模型輸出的假設，其表達式為

p(S)=N[0,C(X,X′)]

(9)

(10)

通過貝葉斯理論，s*的后驗分布表達式為

(11)

由此可知，式(11)具有高斯形式，通過積分后得到的后驗分布是優于預測的輸出，因此可通過下面的公式來表示

(12)

均值和方差的公式分別通過式(13)和式(14)可得

μs*=αTC*

(13)

(14)

式中：C*=C(x*,X)；C=C(X,X)；α=C-1S。

實際計算中高斯過程替代模型采用零均值函數，所以高斯過程替代模型主要由上面提到的協方差函數來確定；因為協方差函數是由其對應的超參數Θ決定，所以超參數中的最優解可通過邊緣釋然函數的最大化來解，其表達式為

(15)

式中，L(Θ)為負對數邊緣釋然函數。

其中，負對數邊緣釋然函數的表達式為

(16)

1.3 不確定性量化分析

使用高斯過程替代模型進行多年凍土場地樁基橋梁地震響應的不確定性量化可由概率分布函數來表征，引起地震響應的不確定性可由統計的均值和方差來表示。由于s(x)遵循高斯分布的特點，因此式(13)和式(14)中的均值和方差同樣具有此特征[23]；表達式為

(17)

(18)

式中：E(·)為期望值；V(·)為變異系數。

針對協方差函數所具有的離散性，式(13)和式(14)可以用如下公式表示

(19)

(20)

將式(19)和式(20)代入式(17)和式(18)中，化簡可得

(21)

V[s(x)]=η2+

(22)

因此，利用以上公式可對高斯過程模型中復雜的高維積分問題轉化為簡單的一維積分問題。

1.4 不確定性量化分析流程

通過高斯過程替代模型計算多年凍土場地樁基橋梁地震響應的統計變量，然后執行不確定性量化。主要思路如下：基于影響多年凍土場地樁基橋梁地震響應的參數，然后確定這些參數所服從分布的特點以此獲得計算的輸入樣本點；其次，針對非線性時程響應的特點，分析不同樣本參數計算所得模型的響應，因此得到樣本參數點的輸入與輸出關系；最后，利用高斯替代過程模型的優勢和特點，對多年凍土場地樁基橋梁的地震響應進行不確定性量化分析，其思路流程如圖1所示。

圖1 地震響應不確定性量化分析Fig. 1 Uncertainty quantification of seismic response

2 多年凍土場地樁基橋梁地震響應不確定量化分析

2.1 數值模擬中p-y曲線的修正

本模擬中的凍土場地樁基橋梁位于多年凍土區，由于凍土區土的特性受溫度的影響十分顯著。當土凍結后，隨著溫度的降低，土的剛度發生顯著的變化。凍土受溫度的作用，土中孔隙充滿冰晶，冰晶與土顆粒之間的黏聚力增加。本模型中所采用的p-y模型是考慮無水情形下的軟巖和硬黏土的p-y模型，借鑒Ge等[24]所采取的方法對本模型中所考慮多年凍土的p-y模型進行修正。由圖2可知，模型中修正的多年凍土p-y曲線一部分為指數方程部分，見式(23)；另一部分為常數部分，見式(24)。

(23)

式中，y≤8ym。

圖2 修正的多年凍土p-y曲線Fig.2 Modified permafrost p-y curve

p=pu

(24)

式中：y>8ym；pu為多年凍土的極限阻力，pu值可通過式(25)求得；ym相當于多年凍土1/2的極限抗力時所對應樁的水平位移。

pu=1.5qub

(25)

式中：qu為模擬中多年凍土的偏應力；b為模擬中樁的直徑。

通常，ym的大小取決于模擬中樁的直徑和多年凍土的剛度，可通過式(26)求得

ym=kmb

(26)

式中,km為常數，對應于多年凍土最大偏應力時應變的50%。

鑒于此，根據地質勘察報告模擬場地的土層剖面自上而下為粉質黏土(0～4 m)、砂礫石(4～11 m)和風化泥巖(11～51 m)。其中在粉質黏土層埋深2 m、砂礫石層埋深6 m和風化泥巖層埋深28 m位置處繪制多年凍土場地的p-y曲線(如圖3所示)。同時，依據多年凍土場地樁基橋梁每月樁的實測溫度曲線，埋深2 m，6 m和28 m位置處所對應的溫度分別取-5 ℃，-2 ℃和0 ℃(如圖4所示)。

2.2 多年凍土場地樁基橋梁數值模擬

多年凍土場地樁基橋梁數值模擬采用開源有限元數值模擬平臺OpenSees進行，前后處理采用MSBridge界面化程序。多年凍土場地樁—土相互作用通過p-y曲線模擬，采取彈性梁柱單元模擬橋面板和蓋梁，橋臺的底部在數值模擬中采用固定連接的方式模擬；基樁和柱墩在地震響應過程中呈現非線性特征，在多年凍土場地樁基數值模擬中采用非線性纖維截面梁柱單元模擬基樁和柱墩(如圖5所示)。

圖3 模型中修正的多年凍土p-y曲線Fig.3 Modified permafrost p-y curves in the model

圖4 8月—11月份監測地溫曲線Fig.4 Monitored ground temperature curve from August to November

圖5 多年凍土場地樁基橋梁有限元模型(m)Fig.5 Finite element model of pile-supported bridge structure in the permafrost region (m)

2.3 輸入參數與關鍵地震響應

基于高斯過程替代模型進行多年凍土場地樁基橋梁地震響應不確定性量化分析，對模擬結果進行整理，選取多年凍土場地樁基橋梁的關鍵部位進行分析。同時，選取有限元模型中排架4進行地震響應分析，選定排架4橋面板的位移和排架4中1號柱墩在埋深1.13 m處的位移、剪力和彎矩作為關鍵的響應。這里選取10種不確定性參數對多年凍土場地樁基橋梁進行模擬。為了量化輸入參數的不同對多年凍土場地樁基橋梁地震響應的影響，基于Sobol序列采樣生成了相應的50個工況的輸入樣本。依據表1中輸入參數的均值和分布，對應得到不同變異參數(coefficient of variation,COV)條件(COV=5%和30%)下50個工況的輸入參數。

圖6為不同變異參數(COV=5%和30%)下采用Sobol序列采樣生成50個工況情況的多年凍土場地樁基橋梁地震響應。圖6(a)為橋梁結構中排架4橋面板的位移響應時程。從圖6(a)可知，變異參數為5%條件下，排架4橋面板的位移在地震作用下隨著加速度的不斷增大，橋梁結構排架4橋面板的位移變化呈現先增大后減小的趨勢；值得注意的是，變異參數為30%條件下，橋梁結構排架4橋面板的位移變化幅度是變異參數為5%時的大約2倍。

圖6(b)為橋梁結構排架4中1號柱墩在埋深1.13 m位置處變異參數為COV=5%和30%條件下的位移響應時程；從圖6(b)可看出，橋梁結構排架4中1號柱墩在埋深1.13 m處位移響應時程的變化與排架4中橋面板的位移變化非常相似；值得注意的是，圖6(a)橋梁結構中排架4橋面板的位移響應時程在不同變異參數下(COV=5%和30%)位移變化幅度大且與此相對應的位移響應時程值高。圖6(c)為橋梁結構排架4中1號柱墩在埋深1.13 m位置處變異參數為COV=5%和30%條件下的彎矩變化響應時程。從圖6(c)可看出，隨著響應時間的增加排架4中1號柱墩的彎矩在變異參數分別為5%和30%情況下彎矩響應幅值先增大后減小。

圖6(d)為橋梁結構排架4中1號柱墩在埋深1.13 m位置處變異參數為COV=5%和30%條件下的剪力響應時程。從圖6(d)可知，橋梁結構排架4中1號柱墩在變異參數分別為COV=5%和30%條件下，隨著響應時間的增加排架4中1號柱墩在埋深1.13 m處的剪力值呈現先增大后逐漸變小并趨于零值，此種現象可進一步表明地震響應的強度會隨著地震作用的時間而變化，隨著地震作用時間的增加剪力值會逐漸變小。

一般而言，從以上分析圖6中橋梁結構排架4中1號柱墩在埋深1.13 m處變異參數為COV=5%和30%條件下的位移、彎矩和剪力時程的變化關系及排架4橋面板的位移響應時程變化關系，可進一步得出橋梁結構中受地震激勵影響的時程響應幅值受參數不確定性的影響較大；當設置橋梁結構參數隨著變異參數由5%增加到30%時，橋梁結構的時程響應幅值也會隨之而增大；由此可得出橋梁結構參數不確定性對多年凍土場地樁基橋梁的地震響應影響顯著。

表1 模型中考慮的不確定性參數

圖6 不同變異參數(COV=5%和30%)下橋梁結構的地震響應Fig.6 Seismic response of bridge structure under different coefficients of variation (COV=5%和30%)

2.4 參數不確定性對體系地震響應的影響

輸入參數不確定性對體系地震響應的影響可通過使用統計分布和變異參數進行量化。具體地，采用6種不同的變異參數(5%，10%，15%，20%，25%和30%)來考查參數不確定性對橋梁結構體系地震響應的影響。圖7為多年凍土場地樁基橋梁地震響應的統計分布。從圖7(a)中可以看出，排架4橋面板的位移隨著變異參數的增大其誤差條逐漸增大，由此可得隨著變異參數的增大對橋梁結構的安全性影響隨之增大。圖7(b)為不同變異參數下橋梁結構排架4中1號柱墩在埋深1.13 m處位移地震響應的統計分布，從圖7(b)可以看出，當結構體系變異參數為5%時，排架4中1號柱墩的位移誤差條較??；隨著結構體系變異參數變大(結構體系變異參數由5%增加到30%)，排架4中1號柱墩的位移誤差條逐漸增大并在結構體系變異參數為30%時達到最大。

從圖7(c)可以看出，隨著結構體系變異參數的增大(由5%增加至30%)，排架4中1號柱墩彎矩地震響應的統計分布變化較明顯；其中，特別是結構體系變異參數為30%時的彎矩誤差條變化最大。由此可表明：結構體系變異參數為30%時相較于5%時對排架4中1號柱墩承受的彎矩變大，此時影響結構的安全系數增大。由圖7(d)可知，隨著結構體系變異參數的遞增，排架4中1號柱墩的剪力誤差條變化的幅度相應的逐漸增大。由此可知，結構體系變異參數為30%其相對應的剪力誤差條變化幅度大，此時1號柱墩承受的剪力值相應的增大；結合圖7(c)和圖7(d)結構體系變異參數由5%～30%時，排架4中1號柱墩的彎矩和剪力的誤差條相應增大。

圖7 多年凍土場地樁基橋梁地震響應的統計分布Fig. 7 Statistics distribution of seismic response of pile-supported bridge structure in the permafrost region

在6種不同的變異參數條件下，橋梁結構體系的地震響應均值會處在一個相對穩定的狀態，而橋梁結構體系地震響應的標準差會隨著橋梁結構體系變異參數的增加而增加，這從側面反映了輸入參數的變異性使地震響應的標準差增大。

通過對圖7變化趨勢的分析與探討，為了進一步更好地詮釋多年凍土場地樁基橋梁地震響應的變異參數對橋梁結構體系參數的不確定性影響，這里采取橋梁結構地震響應的COV和橋梁結構體系參數的COV之間的關系來闡述(如圖8所示)。由圖8(a)可知，隨著橋梁結構體系的COV增加其所對應的排架4橋面板的位移值隨之而增大，結構體系變異參數為30%時排架4橋面板的位移值最大。圖8(b)、圖8(c)和圖8(d)為排架4中1號柱墩在高程-1.13 m處的位移、彎矩和剪力的地震響應COV與結構體系變異參數COV的變化關系；與此同時，隨著結構體系變異參數的增大(由5%增大到30%)，其相對應的地震響應COV值隨之增大(即更大響應的變異性)。從圖8中可以看出，隨著橋梁結構體系COV的增加，其相應的地震響應的COV也是相應的增加?？傊?，從橋梁結構體系參數的COV與響應量COV之間的變化趨勢可進一步表明橋梁結構變異參數的增加直接影響地震響應變異參數的變化。

圖8 橋梁結構體系的COV與地震響應COV關系Fig. 8 Relationship between COV and seismic response COV of pile-supported bridge structure

3 結 論

本文利用高斯過程替代模型，考慮到輸入參數的變異性，進行了多年凍土場地樁基橋梁地震響應不確定性量化分析。首先，根據影響多年凍土場地樁基橋梁地震響應的主要輸入參數，確定這些參數所服從分布的特點；其次，基于Sobol序列的采樣方法，獲得計算時所需的輸入樣本；第三，針對建立的數值模型進行不同樣本參數輸入下模型的響應，從而獲得樣本參數的輸入與輸出關系；最后，利用高斯過程替代模型，對多年凍土場地樁基橋梁的地震響應進行不確定性量化分析。主要得到認識如下：

(1) 高斯過程替代模型可有效表征結構體系輸入參數與地震響應之間的關系，建立輸入與輸出之間相互關系，可以提高分析計算的效率，實現多年凍土場地樁基橋梁不確定性量化分析。

(2) 橋梁結構體系輸入參數隨著COV由5%增加到30%時，橋梁結構的時程響應幅值隨之而增加，橋梁結構參數的不確定性對多年凍土場地樁基橋梁的地震響應影響顯著。

(3) 隨著橋梁結構體系輸入參數的COV增加，相應的地震響應COV隨之增加，橋梁結構體系COV的增加直接影響地震響應COV的變化。