基于主動學習Kriging模型的地下管線抗震可靠度分析

李 娜，侯本偉，杜修力，鐘紫藍，韓俊艷

(城市與工程安全減災教育部重點實驗室(北京工業大學)，北京 100124)

地震作用下地下管線的破壞形式主要有三種[1]：管線接口破壞，管體破壞，在三通彎頭、閥門以及管線與地下構筑物連接處破壞(圖1)。其中承插式管線的接口破壞是最主要的破壞形式[2]。地震對地下管線產生破壞的原因可分為兩類[3]：一類是永久性地面變形，如由斷層運動，砂土液化，滑坡等引起的場地破壞；另一類是地震波傳播過程中導致的土體瞬態變形。永久性地面變形造成的管線平均破壞率較高，但通常發生在局部特殊場地區域；地震波傳播造成的管線平均破壞率相對較低，但地震波動傳播效應對地下管線影響范圍較廣[4]。因此，地震動對研究地下管線結構抗震可靠度具有重要的理論和工程意義。

地震波動作用下管線響應問題的分析方法可分為解析法和有限元法兩大類。共同變位法[5]是最早提出的假定管體與周圍土體同步變形的地下管線抗震解析法；Shinozuka等[6]考慮管-土之間的變形滑動，引入應變轉換系數β來反映場地土體應變與管線應變之間的關系。解析法的優點是概念清晰、方便應用，缺點是地震動荷載以PGV或PGA的函數表征，不能考慮地震動頻譜特性以及不同位置處管線結構、土體參數變化的影響。鑒于此，國內外學者基于有限元法進行地震作用下管線動力時程響應分析。一種是考慮土體存在的精細化有限元模型；另一種是彈性地基梁模型[7]，將埋地管線理想化為彈性地基梁模型，管線與土體之間的相互作用簡化為管-土連接彈簧，從彈簧支座處輸入地震動，求解管線響應。相比于精細化模型，彈性地基梁模型物理意義明確，計算效率高，在地下管線地震分析中被廣泛應用[8-10]。

地下管線一般埋設分布于較大的區域內，很難準確獲取管線模型中每個位置處的外部荷載、管材和土體特性參數等建模參數。因此，進行地下管線的抗震能力分析時有必要考慮模型參數的不確定性。有限元模型進行地下管線抗震可靠度分析時，由于結構響應與模型輸入之間的非線性關系難以使用顯式表達，一般采用Monte Carlo模擬方法(MCS)[11]。當管線失效概率較小時，MCS需要大量的抽樣計算，耗時較長。很多學者使用改進的MCS方法以期達到高效計算的目的。Tee等[12]以及Ebenuwa等[13]提出的改進抽樣方法仍需要大量的樣本點，即有限元模型計算次數很多，計算量減少不明顯；Babu等[14]使用代理模型響應面函數結合一階可靠度方法對地下管-土系統進行可靠性分析，與上述改進抽樣的方法相比，代理模型方法雖然求解功能函數的次數減少，但模型結果的準確性取決于樣本點的選取，容易導致代理模型的計算精度不高。近年來，基于主動學習Kriging代理模型法(active learning based Kriging, AK)，廣泛應用于多種工程可靠度分析[15-17]，主動學習Kriging代理模型方法僅需少量的樣本模型計算便可達到較準確的可靠度分析結果，對于由大量管線構成的實際城市地下管網系統的抗震可靠度計算具有重要的現實應用價值，但該方法在管線抗震可靠度問題的適用性尚待研究。本文采用AK-MCS方法研究地下管線抗震可靠度問題，分析Kriging模型方法的適用性和計算效率，以期為地下管網抗震可靠度計算提供依據。

圖1 地震造成的地下管線破壞Fig.1 Damage to buried pipelines caused by earthquakes

1 地下管線地震響應分析模型

承插式接口地下管線地震反應分析模型如圖2所示。模型由管道、管道接口彈簧和模擬“管-土”相互作用的土彈簧組成，地震動位移時程從土彈簧底部輸入。地震動行波作用下，由于承插式管道接口剛度較小，管道接口承擔了大多數的變形和轉角，而管道結構本身的變形較小。因此，在圖2的分析模型中，地下管道采用彈性梁單元模擬，邊界條件采用與Shi[8]和Kuwate等[9]類似的處理方式，即管線兩端無約束；每根管道單元通過兩個土彈簧與土體連接；在土彈簧底部輸入沿管道軸向的地震動位移時程；采用管道接口相對位移的允許值作為管線結構安全性的判別準則。圖2中的管道長度L通常為3～6 m，管道接口長度與管道單元長度相比可忽略不計，但為了確保有限元模型數值計算的穩定性，參照Shi[8]和Liu等[18]的方法，本文分析模型中設置管道接口長度為ΔL。沿管道軸向的土彈簧采用圖3(a)所示的雙折線模型[19-20]，圖3(a)顯示了土體對單位長度管道的剪切力f和管道-土體相對位移之間的關系，在管-土相對位移≤δs(土彈簧的彈性極限位移)時線性上升到fy(管-土軸向彈性極限摩擦力)，當管-土相對位移>δs后，土體傳遞給單位長度管道的剪切力為常數fy。fy的表達式為[19]：

(1)

式中：K0為土的靜止側向土壓力系數(回填溝管道一般為0.5≤K0≤1.0);γ為土體容重(N/m3);H為管道中心線的埋深(m);D為管道外徑(m);φ為砂土內摩擦角(通常取值范圍為25°～45°);Su為周圍土體的不排水抗剪強度，對于淺埋且回填壓實黏土中的管線,Su的取值范圍為25～60 kN/m2;α是決定管道與土體黏結程度的黏附因子，可由Su的數值估算出來[19]。

圖2 地下管線有限元模型示意Fig.2 Schematic of finite element model of buried pipeline

圖3(b)為反映管道接口軸向拉拔力與相對位移關系的折線模型，當接口兩側管道端點相對位移≤δJ(接口彈性極限位移)時，接口受到的拉拔力呈線性上升；當接口相對位移>δJ后，接口拉拔力為彈性極限拉力常數(FJ)。

圖3 力與相對位移的彈塑性關系Fig.3 Elastoplastic relationship between force and relative displacement

一般認為埋設在土體中的地下管線結構在地震波動作用下的結構響應主要受管線周圍土體變形的控制，管線結構的慣性力可以忽略不計[5]。簡化計算公式得到的地震波動作用下承插式管線接口的張開量可表示為[3]

(2)

式中：εg為土體應變;L為管道長度(m);Va為地震波傳播方向上的最大水平土體質點速度(m/s);Ca為地震波的傳播速度(m/s)[3,10]。本文中地下管線的非一致地震動主要考慮地震動的行波效應，地震波在土體的運動傳播過程中保持波形和振幅不變，沿地震波傳播路徑上距離為ΔX的兩個土體質點振動的加速度、速度和位移時程均只存在一個時滯Δt[10,21]：

(3)

式中ΔX為土體中兩個位置點間的直線距離(m)。

2 地下管線抗震可靠度分析方法

2.1 隨機變量及其分布類型

在地下管線的抗震分析模型中，地震動、接口彈簧、土體類型、管線結構等參數的隨機性導致了管線抗震安全性分析結果的不確定性。Manolis等[22]在研究地下管線在土體中的動力響應時，將土體密度和摩擦角作為隨機變量；Ebenuwa等[13]進行埋地鋼管的抗震可靠度分析時，將土體黏聚力、容重，管道容重、埋深、楊氏模量、壁厚等參數作為隨機變量；Wijaya等[23]進行不確定性條件下埋地分段式管線地震反應分析時，將黏性土的密度、黏聚力和內摩擦角視為相關的隨機變量，同時認為管道滲漏極限接口容許位移為正態分布隨機變量；Jahangiri等[24]對15組不同埋深與管徑比、場地條件和鋼材型號的埋地鋼管進行了地震波動作用下的IDA分析；Shi[8]在研究地震波在埋地分段管道中的傳播效應時，考慮了峰值地面應變、管道接口特性的變化規律，認為管道滲漏極限接口容許位移是接口插入深度的0.52倍，變異系數為10%。綜上，管線結構性能在很大程度上取決于管土相互作用和接口容許位移。綜合上述文獻的參數選取，本文選擇管道參數(接口允許位移、管道埋深)，土體參數(容重、內摩擦角)作為管線模型可靠度分析的隨機變量。

2.2 地下管線結構安全判別準則

圖2所示的管線模型中，在地震動位移時程荷載作用下，不同管道接口的位移響應是隨時間變化的隨機過程g(x,t)：

g(x,t)=R(x)-S(x,t)

(4)

式中：R(x)表示管道接口容許位移值；S(x,t)表示地震動荷載輸入下，t時刻的管接口的相對位移響應值;x為管線結構地震反應分析模型中不確定性參數對應的隨機變量。

若地震動位移時程的持時為Ts，根據動力可靠度的首次超越失效準則，式(4)中隨機過程g(x,t)在時段[0,Ts]內進入失效域g(x,t)≤0的概率表示為

Pf=prob{g(x,t)≤0, ?t∈[0,Ts]}

(5)

Pf=prob{G(x)≤0}=

(6)

式中q(x)為隨機變量集x的聯合概率密度函數。

2.3 基于Monte Carlo模擬的可靠度計算

(7)

式中：x(k) =[x1(k),x2(k),x3(k), …,xn(k)]為第k次模擬抽樣得到的x樣本值；I取值為{0,1}的函數，當G(x(k))≤0時，I[G(x(k))]=1，否則=0。

3 主動學習Kriging代理模型計算管線可靠度

3.1 構建Kriging代理模型

Kriging代理模型一般包含兩部分：回歸部分和隨機誤差[25]。對于k個樣本點初始樣本集(design of experiment, DoE)，輸入變量樣本點X=[x1,x2, …,xk]T及其對應的響應值(輸出變量)Y=[S(x1),S(x2), …,S(xk)]T，采用Kriging模型表示的輸出變量S(x)與輸入變量x之間的表達式為

(8)

式中：每個樣本點中輸入變量的個數為n，x=[x1,x2, …,xn]T，f(x)=[f1(x),f2(x) , …,fm(x)]T為回歸多項式基函數向量，m為回歸多項式的數量；β=[β1,β2, …,βm]T為多項式參數向量；z(x)為服從正態分布N(0,σ2)的隨機過程，其協方差方程為

cov(z(xi),z(xj))=σ2R(xi,xj)

(i,j=1,2,3,…,k)

(9)

式中R(xi,xj)為任意2個樣本點xi和xj的空間相關方程，通常采用高斯相關方程，其表達式為

(10)

(11)

式中R為相關矩陣，元素Rij=R(xi,xj) (i,j=1, 2, …,k)。

基于給定的樣本點，多項式參數向量β與隨機過程方差σ2的估計值計算式分別為

(12)

式中F為由回歸多項式函數值構成的矩陣，其中元素Fij=fj(xi)，且有i=1, …,k；j=1, …,m。

式中u=FTR-1r0-f。本文進行管線結構抗震可靠度分析時，式(8)～(13)中應用Kriging模型構建和模型預測應用的過程，采用MATLAB軟件支持開發的Kriging模型工具包DACE實現[26-27]。

3.2 主動學習更新Kriging代理模型

在采用Monte Carlo方法計算結構可靠度時，采用式(8)～(13)的Kriging代理模型預測抽樣樣本的結構響應，可以減少抽樣樣本結構有限元模型計算次數，從而加速Monte Carlo方法的求解速度。為了保證可靠度計算的準確性，一般需要較多的初始樣本點{X;Y}建立精度較高的Kriging模型；而初始樣本點{X;Y}的建立，仍然需要較多的計算量。為了提高Kriging模型的精度和計算效率，近年來發展了采用主動學習方法逐步提高Kriging模型精度的方法，僅需少量的樣本就可構造精度較高的Kriging模型，也即：根據學習函數選擇最佳樣本點，增補該樣本點后更新擬合Kriging模型，迭代多次更新以提升Kriging代理模型的精度。

為使Kriging模型不斷優化，Echard等[28]提出了類似可靠度指標的定義學習函數U評價樣本點x對于現有Kriging模型的潛在影響：

(14)

(15)

本文采用Echard提出的min{U(x)}≥2作為Kriging模型停止學習準則，通過式(15)得到x*。

3.3 基于AK-MCS的管線抗震可靠度計算

在管線結構系統的有限元模型中，管線接口響應與荷載、場地和管道屬性參數之間的關系難以顯式表達。基于2.3節和3.2節的描述，采用主動學習Kriging代理模型(active learning based Kriging, AK)結合Monte Carlo隨機模擬的方法(AK-MCS)計算地下管線的抗震動力可靠度的流程如圖4所示。具體計算步驟如下：

步驟1：在隨機變量的x取值空間中使用拉丁超立方抽樣方法(LHS)隨機生成包含m個自變量樣本的輸入樣本集X。根據初始樣本的物理參數取值，利用ABAQUS軟件構建管線有限元模型，并計算管線接口的最大位移響應值，作為初始樣本X對應的功能函數響應值Y，此響應值為所有時間范圍內接口位移最大值。其中，初始樣本點的個數m≥(n+1)(n+2)/2[29]，n為隨機變量個數。初始樣本X和響應值Y作為Kriging模型的訓練樣本集DoE。

步驟2：基于訓練樣本集DoE，采用MATLAB的DACE工具箱實現式(8)～(12)的計算過程，構建管線接口位移響應安全狀態的Kriging模型。其中，Kriging模型的回歸部分為一次多項式，相關函數為高斯函數。

步驟3：使用MC抽樣產生N個自變量樣本點XN，利用構建的Kriging模型求得自變量樣本點XN對應的管線接口位移響應預測值和預測方差(式(13))。計算這些樣本點的學習函數值U(式(14))，并從中挑選學習函數最小值(Umin)對應的樣本點x*(式(15))作為最佳樣本點。

步驟4：如果最佳樣本點自變量x*的學習函數值滿足學習停止條件Umin=U(x*)≥2，則構建Kriging模型的主動學習過程結束，轉步驟5；否則，建立最佳樣本點自變量x*對應的管線有限元模型，計算管線接口位移響應真實值S(x*)，將{x*,S(x*)}加入Kriging模型訓練樣本集DoE中，轉步驟2。

步驟5：主動學習過程結束，當前Kriging模型即為滿足精確需求的模型。根據步驟3計算得到XN對應的結構響應值，利用式(7)求得管線地震失效概率Pf的估計值。

圖4 AK-MCS方法求解管線抗震可靠度流程圖Fig.4 Flow chart for seismic reliability evaluation of pipeline by AK-MCS method

4 算例分析

4.1 算例參數

根據圖2所示的管線模型，在ABAQUS有限元軟件中建立了承插式接口球墨鑄鐵管線的地震反應分析模型，沿管線軸向輸入地震動荷載。管線結構參數和周圍覆土的參數見表1，其中管線周圍土體為無黏性回填土，根據式(1)可算得土彈簧的軸向彈性極限力fy=20.5 kN/m。管道外徑61 cm和周向土體抗拉能力0.15 kN/mm，計算得接口的彈性極限拉力FJ=287 kN。土彈簧和接口軸向彈簧建模采用圖3所示的折線模型，土彈簧的彈性極限位移δs=0.3 cm，接口軸向彈簧的彈性極限位移δJ=0.25 cm[8]。根據2.1節的分析，本節進行管線抗震可靠度算例分析時，選擇的模型不確定性參數及其分布特征見表2。

表1 管線有限元模型參數Tab.1 Finite element model parameters of pipeline

表2 管線模型中隨機變量參數及其概率分布參數Tab.2 Random variables and corresponding distribution properties in pipeline model

為分析地震動頻譜特性對管線地震反應的影響，本節沿管線軸向輸入正弦波、EL-Centro記錄以及汶川安縣地震動位移時程記錄。其中，正弦波位移時程的參數為沿管道軸向最大速度Va=30 cm/s(中國地震烈度VIII度)，傳播速度Ca=120 m/s，周期T=3.5 s，持時28 s，此參數與1985年Michoacan地震[30]中產生大量管線破壞的面波特征相似。EL-Centro地震動時程取自1940-05-19 Imperial Valley 南北向地震動記錄，持時40 s，將地震動速度時程調幅至Va=30 cm/s；汶川安縣地震動時程取自2008-05-12 安縣地震動記錄，將地震動速度時程調幅至Va=30 cm/s，積分得到位移時程。為體現非一致地震動輸入的影響，圖2所示的管線不同位置處的土彈簧底部輸入的地震動之間的時間延遲采用式(3)計算。

本節算例中不同工況的模型和荷載參數見表3。為了驗證本文分析模型的正確性，工況1采用與Shi[8]相同的模型參數，即：管線總長3 660 m，單根管道長度4.55 m，地震動位移荷載為沿管線軸向分布的靜力荷載。參數研究表明接口長度對數值模擬結果的影響很小，為了數值模擬的穩定性，模型中接口長度設置為2.5 cm。管道服役年限增加、接口膠圈老化等原因均可能導致管道接口抗拉拔能力下降，Shi[8]和Zhong等[10]研究了管線模型中間存在一個薄弱接口對管線結構響應的影響。本節假定管線模型中間位置處的接口為薄弱接口，工況2～8計算了中間位置處接口的彈性極限拉力為FW={0.2, 0.4, 0.6, 0.8, 1.0}×FJ的不同結果。根據《連續鑄鐵管》[31]，國內承插式接口球墨鑄鐵管道的有效長度一般為6 m，因此工況3～8中管道長度取為6 m。

表3 地下管線地震反應分析工況Tab.3 Analysis for seismic responses of buried pipeline in different cases

4.2 地下管線的地震響應

確定性模型參數工況1的計算結果見圖5。圖5(a)繪制了沿管線長度上各點處的輸入位移和管線位移響應(土彈簧兩個端點的軸向位移)。由圖5(b)可看出管線隨著土體位移一起變形，管道中間位置土節點的位移等于管道節點位移。圖5(c)是沿管線長度方向的管道應變，管道應變從0開始線性增加，直到接口軸向力達到接口的抗拉能力時，管道應變趨于穩定。圖5(d)是沿管線長度方向的接口張開量，接口張開量從0到0.25 cm(接口彈性位移極限)呈線性增加，之后的變化呈正弦曲線趨勢，最大接口位移約1.10 cm。工況1與Shi[8]的模型采用相同的參數，計算結果與Shi[8]相同，驗證了本文管線分析模型的正確性。

本文采用的管線分析模型(圖2)的兩側端點處(0 m和3 660 m)并未設置端部約束，對于考慮位移時程荷載非一致輸入的工況1和工況2，管道接口最大張開量分布如圖6所示。由圖6中工況1的結果可知，接口最大張開量從管道兩端開始逐漸增大，經過一定距離(約20 m)后達到穩定值，管線中部接口最大張開量值基本相同(約為1.10 cm)(圖6)，這與Liu等[18]的結論一致。根據圖6中的結果，可知兩個工況的邊界條件影響區域的接口位移變化規律相同、中間區域的管線接口最大位移響應相等，可認為本文承插式接口管線模型兩側約20 m范圍為邊界條件影響區域。若將管線建模長度選擇在300 m左右，仍可得到不受邊界條件影響的管線接口的位移響應，從而縮小管線模型的計算量。

圖5 工況1計算結果Fig.5 Pipeline response in case 1

圖6 工況1與工況2(FW=FJ)沿管線長度的接口張開量比較Fig.6 Comparison of axial joint opening between case 1 and case 2(FW=FJ)

對于管線中部存在薄弱接口時，圖7給出了工況2～5中薄弱接口最大張開量隨接口剛度變化的關系。由圖可知，薄弱接口張開量隨著FW減小而近似線性增加；與工況2計算結果相比，工況3管道長度6 m對薄弱接口張開量的影響更大；工況3(正弦波)、工況4(EL-Centro波)、工況5(汶川安縣波)的弱接口張開量不相等，而根據式(2)可算得三種工況的解析結果相等，說明接口張開量不僅受地震波峰值影響，頻譜特征不同也會引起接口張開量不同，體現出進行管線地震動力時程響應分析的重要性。

4.3 地下管線的抗震可靠度

分別采用2.3節的標準Monte Carlo模擬(MCS)方法、3.3節的基于主動學習Kriging模型的Monte Carlo模擬方法(AK-MCS)計算工況6～8中管線模型的抗震可靠度。在工況6～8的可靠度算例模型中，隨機變量參數個數為n=4。根據3.3節步驟1的描述，構建Kriging代理模型所需要的最小為初始樣本點為m≥15。為了增加初始樣本點的代表性，本節采用拉丁超立方抽樣產生m=20個初始樣本點建立初始Kriging模型。步驟3中MC抽樣個數為N=1 000，生成自變量樣本點，利用Kriging模型預測響應值(接口最大張開量)。

圖7 工況2～5中間薄弱接口張開量Fig.7 Opening of middle weak joint in cases 2-5

表4給出了在工況7(FW=0.4FJ)中使用AK-MCS方法計算不同樣本點的管線失效概率及直接采用MCS方法的失效概率結果，將MCS方法計算結果視為“精確解”。有限元計算次數是指需要調用ABAQUS有限元軟件執行管線計算的次數。

由表4可知，當采用初始樣本個數m=20，基于此樣本點所建立的Kriging模型進行可靠度評估計算的失效概率為0.112，計算的相對誤差為9.80%，誤差較大，說明所建立的模型精度有待進一步提高。該模型通過主動學習，增加新的樣本點至初始樣本集，重新構建Kriging模型。當樣本點達到23時，其失效概率相對誤差降低至3.92%，此時Kriging模型符合精度要求。為驗證所建立的模型精確性，圖8比較了Kriging模型預測值與真實響應值。圖8(a)是通過初始樣本點構建的Kriging模型，圖8(b)是經過主動學習后的Kriging模型。可以看出，經過主動學習的Kriging模型預測得到的響應值基本接近真實值，兩者響應值的相對差值也逐漸變小，說明構建的Kriging模型是準確的。

表4 MCS與AK-MCS的可靠度結果(工況7(FW=0.4FJ))Tab.4 Reliability of MCS method and AK-MCS method(case 7(FW=0.4FJ))

表5為MCS方法和AK-MCS方法的地下管線抗震可靠度計算結果對比。在FW=FJ時，從圖7可知此時管線接口最大張開量遠小于接口允許位移(5.2 cm)，三種地震波的失效概率均為0。在不同地震波作用下，隨著弱接口抗拉能力的減小，管線的失效概率均逐漸增大。需要說明的是，盡管表5中工況6～8采用相同的地震動參數(Va=30 cm/s和Ca=120 m/s)，由于地震動荷載時程的周期成分不同，導致管線的地震反應和抗震可靠度計算結果差異較大。而按照式(2)的簡化方法計算接口張開量，計算結果保持不變，不能體現管線薄弱接口剛度變化、地震動周期成份對管線地震反應的影響，因此建立地下管線結構的地震動力時程響應分析模型是進行管線抗震可靠度分析的必要支撐。

表5 AK-MCS和MCS計算的管線失效概率Tab.5 Pipeline failure probability evaluated by AK-MCS method and MCS method

由表5可以看出，AK-MCS方法得到的管線失效概率Pf與MCS方法的Pf相對誤差(|MCS-AK-MCS|/MCS)為0.00%～3.92%，而AK-MCS方法的計算時間僅為MCS方法的2%；同時，其計算有限元次數僅為MCS方法的2%左右。因此，對于地下管線抗震可靠度計算問題，AK-MCS方法是一種準確且高效的計算方法。

5 結 論

地下管線結構的抗震可靠度分析是管線抗震安全性評估的重要內容。本文主要考慮土體、管道結構等參數的不確定性，采用主動學習Kriging模型可靠度方法(AK-MCS)對地下承插式球墨鑄鐵管線進行抗震可靠度評價，得出的主要結論如下：

1)本文算例結果表明，AK-MCS法所得的管線失效概率與標準Monte Carlo模擬方法的相對差值在5%以內；AK-MCS法中有限元模型計算次數和計算時間約為標準Monto Carlo模擬方法的2%左右。表明采用AK-MCS法計算地下管線抗震可靠度是可行且高效的。

2)當地震烈度為Ⅷ度(VPG=30 cm/s)時，新建成的承插式球墨鑄鐵管線接口失效概率為0，管線處于安全狀態；隨著管道服役年限增加、接口膠圈老化等，管線接口的抗拉拔能力下降到原接口的60%時，管線失效概率顯著上升。因此，進行現役管線的抗震可靠度分析時，需要注意管線接口模型剛度的合理性。

3)由于地震動時程的頻譜成份不同，不同地震動導致的地下管線地震反應結果差異顯著。僅利用公式簡化方法計算的接口張開量，不能體現管線薄弱接口剛度變化、地震動周期成分對管線地震反應的影響，因此建立地下管線的地震動力時程響應分析模型是進行管線抗震可靠度分析的必要支撐。