跨活斷層隧道抗位錯精細化分析方法研究

趙伯明，趙天次，周玉書

(1.北京交通大學 城市地下工程教育部重點實驗室, 北京 100044；2.北京交通大學 土木建筑工程學院, 北京 100044)

活斷層在我國分布廣泛，突發性的斷層錯動是產生地震的主要根源，長期以來給人類的生命財產和結構工程的安全性帶來了很大威脅。大量研究表明，活斷層錯動導致的同震地表破裂、土體變形、強地面震動及其他次生地質災害均會給穿越活斷層隧道帶來危害[1-11]。然而在大型工程的選線和設計中不可避免的會穿越或臨近活斷層，因此跨活斷層隧道的抗位錯研究具有重要的理論意義和工程價值。

目前針對穿越活斷層隧道抗位錯分析國內外開展了諸多研究，研究方法包括現場調研、模型試驗、解析計算和數值模擬等。

現場調研可以直接得到最真實的一手資料，通過震后短時間內的儀器監測可得到斷層附近的地表同震位移[12]，運用回歸分析手段將實地調研數據進行處理，建立地表同震位移與震級、上覆土層厚度等指標的定量關系[13-17]。但該方法受限于觀測點的數量和密度，不同區域的地質構造和土層性質差異也會導致回歸關系式的差異[18]，且大多數研究重點關注了地表的同震位移，對隧道等地下工程而言，土體內部的變形對結構受力影響更大，斷層錯動導致的圍巖狀態變化同樣值得關注。

部分學者針對斷層錯動對隧道結構的影響進行了模型實驗研究[19-24]，但在實驗中由于模型尺寸的限制，對斷層錯動的模擬通常簡化為單側底板的抬升，并不能模擬真實的土體變形情況。

在基于彈性半空間的解析解法方面，學者們[25-27]應用簡化的斷層模型和位錯理論，給出了半無限空間中矩形滑動斷層的地震位移場的解析解。該方法推導清晰，操作便捷，其中應用最為廣泛的是Okada提出的均勻半無限空間解析解法，但其缺點是對于多層介質模型或三維復雜情況并不適用。

隨著地質資料和勘測手段的進步，Okada的解析解法已經略顯粗糙，在均勻半無限空間解析解基礎上進一步發展而來的層狀半空間中的靜態彈性變形計算方法考慮了地球的層狀構造，位移矢量需要在柱坐標系下利用柱諧函數來表達，并通過傳播矩陣的算法求解[28-29]，近年來已經逐漸成為主流趨勢。

對于地球層狀構造對同震變形的影響程度，文獻[30-31]的研究表明：在地殼內發生的地震，在震源距不大于100 km的情況下，層狀構造的影響可達到20%。

因此，隨著計算理論的發展，人們對斷層位錯發生后土體的變形了解逐步加深，有限元模型中的輸入也應盡量接近真實的構造環境，以提高計算結果的可靠性，更加精細化的位錯輸入方法有待提出。

本文首先對現有的跨斷層隧道抗位錯計算方法進行誤差分析，宏觀上分析了土體層狀構造和斷層破裂面上位錯分布導致的同震變形場差異，微觀上分析了不同輸入方法導致的土體變形差異，揭示了現有隧道抗位錯計算中同震變形輸入方法的不足。然后基于正交化格林函數方法編制程序，對層狀半無限空間中的靜態彈性變形進行計算，結合Abaqus的input文件將斷層錯動導致的土體變形場以位移邊界的形式施加到有限元模型，并采用Matlab建立了輸入流程的GUI界面，使建模流程大為簡化。最后將該方法應用于成蘭鐵路某隧道的三維抗位錯分析中，實現了隧道工程抗位錯分析中斷層錯動的精細化輸入。

1 現有跨斷層隧道抗位錯機理研究

我國西部地區地質條件復雜，活斷層分布密集，隨著我國西部地區鐵路隧道逐漸增多，部分隧道不可避免的需要穿越眾多的活斷層。當隧道工程處于斷層附近時往往需要進行抗位錯分析，以得到隧道結構在斷層錯動作用下變形效應和破壞機理，隧道工程與活斷層之間的的空間關系見圖1，黑色單箭頭代表斷層的錯動。

圖1 跨斷層隧道空間位置關系

在隧道工程的抗位錯分析中，通常截取有限區域為計算區域，并設置人工邊界來模擬活斷層的位錯[10,32]，但現有的建模過程中往往存在著過度的簡化，簡化分析模型見圖2。

圖2 現有研究中的簡化模型[32]

由圖2可知，上述簡化建模方法中將斷層破裂面直接置于隧道中部，通過斷層上盤底面的整體抬升模擬逆斷層的錯動。而實際地層環境中大部分斷層的破裂面并未直接截斷隧道，對于埋深不為零的隱伏斷層或破裂面未直達地表的斷層而言，該模擬方法顯然是不適用的。此外，上述研究中將斷層破裂面上的位錯量設置為均勻分布，與實際斷層面上的位錯分布存在較大的差別。通過分析可知建模過程中的誤差來源如下：①土體的層狀結構；②破裂面上的位錯分布；③斷層錯動的輸入方法。本節以典型地震為例，通過對比分析得到各誤差源的影響大小，說明了建立精細化位錯輸入方法的必要性。

1.1 土體層狀結構對地表同震變形的影響

土體的層狀構造在同震變形場的計算中存在影響，為量化這一誤差，本節選取典型地震2013年Mw6.7級蘆山地震作為算例，對成層土體模型和均勻土體模型下的地表同震變形進行對比分析，研究土體層狀結構對同震變形場的影響。

由王衛民等[33]給出的蘆山地震反演結果可知，發震斷層為一條埋深2 km的隱伏逆斷層，斷層破裂面并未延伸到地表，斷層參數見表1。斷層參數及坐標系如圖3所示。其中斷層長度為L，寬度為W，傾角為β，斷層面上位錯滑動角為α，滑動量為u。

圖3 斷層空間坐標系

表1 蘆山地震震源參數

成層土體模型的土層劃分參考參考蘆山地震震源區地殼地震波速度模型[34]及前人的研究成果[35-41]，綜合給出了蘆山震區的成層土體模型，自地表向下共劃分了4個不均勻波速層，土層參數見表2。均勻土體模型與成層土體模型的最上層土采用相同的力學參數。

表2 蘆山地震震區地殼結構模型

首先計算均勻土體模型中蘆山地震導致的地表同震變形，豎向位移等值線和水平位移矢量分別見圖4，其中虛線矩形框表示斷層的地表投影。

圖4 均勻位錯模型地表豎向、水平位移矢量圖

地表位移的分布受斷層產狀影響，豎向位移的數值總體上大于水平位移。由圖4(a)可知，地表隆起主要分布在斷層上盤區域即Y軸正半軸，沉降集中下斷層下盤區域，同震變形在斷層上斷線的地表投影附近變化最為劇烈。豎向位移的隆起峰值為0.802 m，沉降峰值為0.062 m，這是由于斷層傾角較小，斷層面上的位錯量經過土層傳導出現了一定的偏移。

由圖4(b)可知，上盤和下盤的土體均向斷層上斷線的地表投影擠壓，隨著與上斷線地表投影線距離的減小，水平位移值逐漸增大，峰值為0.338 m。

基于分層土體模型的蘆山地震地表同震變形計算結果見圖5，從結果來看，地表豎向位移和水平位移的分布和數值都與GPS監測數據[42-43]給出的結果相近。

圖5 分層位錯模型地表豎向、水平位移矢量圖

對比圖5兩種模型的地表同震變形計算結果，考慮地殼層狀構造之后地表的豎向位移和水平位移在分布規律上絕大部分都保持一致，豎向位移隆起峰值為0.775 m，沉降峰值為0.046 m，水平位移峰值為0.395 m，但是在峰值位置上有所差異。均勻土體模型的隆起峰值和沉降峰值都要大于分層土體模型，但水平位移峰值小于分層土體模型，這一差異是由均勻土體模型較分層土體模型更軟的土體力學參數導致的。為清晰的展示兩者差異，將兩種模型的計算結果取差值，見圖6。

圖6 分層模型和均勻模型的地表豎向、水平位移差值

由圖6(a)可知，豎向位移的差異主要集中在上斷線投影附近，地表隆起存在一條沿斷層走向的0.1 m以上分布帶，地表沉降則對稱分布在上斷線投影，量值上僅0.05 m左右。這說明均勻模型具有更大的隆起和沉降峰值，豎向位移最大差值為0.308 m，約占分層土體模型豎向位移最大值的39.7%。

由圖6(b)可知，水平位移的差異同樣體現在上斷線投影處，差異帶寬度10 km左右，最大差值0.236 m，達到分層土體模型水平位移最大值的59.8%。

由蘆山地震震例的計算結果可見，土體的分層結構對地表同震位移的影響是不可忽略的，位移差異達40%～60%，尤其是在靠近斷層上斷線地表投影區域，需要考慮地震發生地區地殼層狀構造的影響。

1.2 破裂面上位錯量分布對地表同震變形的影響

在真實的地震中，斷層破裂面的幾何形狀與位錯量分布都很復雜。在以往的跨斷層隧道計算中，往往假設斷層面位錯是均勻分布的，但實際斷層面上的位錯分布是非一致的。吳忠良[44]通過對多個地震斷層面上的滑動量調查統計發現斷層面上的位錯分布是高度不均勻的；Banerjee等[45]利用41個連續觀測的遠場GPS資料對Sumatra-Andaman地震斷層的滑動量進行計算，發現斷層上的最大滑動量出現在斷層中部附近。

在經典半無限空間的位錯理論中，Okada均勻土體模型的有限矩形位錯源公式是由點源公式在斷層面上積分所得，分層土體模型中有限斷層源是由離散的點源格林函數線性疊加而得，兩種模型均可對斷層進行離散化，使破裂面上的位錯量接近真實的位錯分布，再求得對應的同震變形場。

為分析斷層破裂面上位錯量分布對同震變形場的影響，選取典型地震1989年美國Mw6.95級Loma地震分別計算位錯量均勻分布和真實分布兩種情況的地表同震變形，Loma地震發震斷層參數見表3，Loma地區的土體參數見表4[47]。

表3 Loma地震震源參數

表4 Loma地震震區地殼結構模型

沿走向和傾向將發震斷層劃分為2 km×2 km的子斷層，每個子斷層上的位錯量反演結果見圖7[46]。

圖7 Loma地震斷層破裂面位錯量分布(單位：cm)

圖7中每個子斷層上位錯量的精確反演結果反映了斷層面上位錯的真實分布情況，再求取各個子斷層面上位錯量的平均值作為均勻分布模型的位錯量，采用分層土體模型分別計算Loma地震兩種位錯分布模式下的地表同震位移。

真實位錯分布模式下地表三個方向的位移等值線分別見圖8。

圖8 真實位錯分布下豎向、X向、Y向地表位移等值線(單位：m)

由圖8可知，從結果來看地表位移分布符合斜滑斷層的運動機制。在斷層傾向抬升作用和走向擠壓作用的疊加影響下，豎向位移有整體向左移動的趨勢，豎向位移峰值在斷層中部左側約7 km的位置，峰值量約為0.7 m。X方向位移主要受走滑分量的控制，表現為上盤區域沿走向向左移動，下盤區域沿走向向右移動，X向位移峰值在斷層中部左側約12 km的位置，峰值量約為0.4 m。Y方向位移受到走滑分量和傾滑分量的共同作用，表現為大部分區域向Y正方向移動，只有斷層右側小部分區域向Y負方向移動，Y向位移峰值在斷層中部左側約8 km的位置，峰值量約為0.3 m。

為了量化位錯分布模式對地表同震位移的影響，將均勻分布模式下的結果與真實分布模式的結果取差值進行分析，三個方向上的位移差值分布分別見圖9。由圖9可知，均勻分布下的三方向位移在整個斷層區域內都與真實分布有較大的差別，誤差最大的位置對應位移峰值所在位置。

圖9 位錯均勻分布與真實分布的豎向、X向、Y向地表位移差(單位：m)

考察斷層走向中線剖面位置處的地表同震位移見圖10。由圖10可知，三個方向的位移值都隨著遠離斷層而逐漸減小，隨著與斷層距離的減小，均勻分布模式會產生較大的誤差，具體數值見表5。

圖10 斷層走向中線剖面的地表同震位移對比

表5 不同位錯分布模式的位移峰值誤差

由表5可見，均勻分布模式在水平和豎向都會給地表同震位移帶來誤差，水平同震位移峰值誤差在X方向達到36.6%，在Y方向達到8.4%；豎向同震位移峰值誤差在Z方向達到9.6%，這些誤差放大了斷層錯動的破壞效應。

綜上所述，將斷層面上的位錯量按均勻分布考慮會給計算結果帶來誤差，在計算土體的同震位移場時，采用反演得到的真實位錯分布具有一定的必要性。

1.3 斷層錯動輸入方法對地表同震變形的影響

在求得斷層錯動導致的土體同震變形場之后，將該同震變形場準確合理的輸入截取的計算區域是保證跨斷層隧道抗位錯分析準確性的重要一環。現有研究中通常采用斷層上盤底面的整體抬升模擬逆斷層的錯動，見圖2，這是一種過度簡化的手段，顯然會帶來一定的誤差。本小節選取矩震級Mw7.0級典型地震對斷層錯動導致的模型底面同震變形進行計算，并將模型底面的真實變形與簡化方法的輸入值進行對比，量化了簡化方法的誤差，說明了選擇合理斷層錯動輸入方法的必要性。

為使對比結果更為清晰，忽略土體分層和位錯分布的影響，按均勻土體模型且位錯均勻分布計算，土層參數設置為表2中的第一層土，選取逆斷層進行分析，由斷層破裂參數和矩震級Mw之間的統計公式[46,48-49]得Mw7.0級典型地震的震源參數見表6。

表6 Mw7.0級典型地震震源參數

選取斷層上斷線中軸面為觀察剖面，斷層模型及觀測剖面位置如圖11，分別計算深度為0、50、100、150、200 m處觀察剖面的土體水平和豎向同震變形值，將各深度真實的土體同震變形值和簡化輸入變形值對比，見圖12。

圖11 斷層模型及觀測剖面

由圖12(a)可知，豎直方向土體的真實同震變形與簡化方法中輸入的土體變形值趨勢一致，隨著深度的增加簡化值與真實值的差值逐漸減小，上盤區域差值約0.2 m.下盤區域差值約0.25 m，誤差為真實變形值的25%左右；由圖12(b)可知，而水平方向土體的真實同震變形與簡化輸入變形值差異顯著，簡化方法中認為斷層上盤固定不動僅下盤抬升的輸入方式顯然是不準確的，上盤和下盤的輸入差值都為0.4 m左右，水平變形輸入誤差達80%左右。

2 跨斷層隧道抗位錯精細化計算

由上節分析可知，土體層狀結構、破裂面上的位錯量分布及錯動輸入方法均會導致地表同震變形的計算出現誤差，其中影響最為顯著的是錯動輸入方法，誤差達25%～80%；土體層狀結構影響次之，誤差為40%～60%；破裂面上的位錯量分布影響最小，僅10%左右，該誤差在位錯量反演不明確的條件下可忽略。

因此，本節提出了一套精細化的位錯輸入方法，首先由斷層位錯出發，求得斷層錯動導致的模型邊界上的同震變形場，再將該同震變形場以位移人工邊界的形式施加至有限元模型的邊界上，以實現斷層錯動的精確輸入。以成蘭鐵路跨越龍門山破裂帶某隧道為例，探究隧道在斷層錯動作用下的破壞機理及力學特性，為隧道工程抗位錯措施的建立提供指導。

2.1 精細化位錯輸入方法的實現過程

早期研究中將斷層上盤或下盤的移動簡化為模型單側底面的抬升，這是較不合理的一種簡化輸入方法。本文提出的精細化斷層錯動輸入方法建模思路見圖13。

圖13 精細化輸入方法建模思路

計算流程如下：首先輸入斷層參數和土層參數，求得斷層錯動導致的有限元模型區域的同震變形場，再對計算結果進行提取，獲得模型邊界處的變形值并寫入input文件，最后通過input文件將同震變形場施加到Abaqus有限元模型，以實現斷層錯動的精確輸入，自編程序計算流程見圖14。

圖14 計算流程圖

2.2 跨活斷層隧道數值計算模型

以我國西部成蘭鐵路穿越北川—映秀活斷層某隧道為例進行抗位錯分析，北川—映秀是龍門山中央斷裂帶的主斷裂，也是汶川地震的發震構造之一，本文將其簡化為逆沖斷層，斷層參數見表6。

由圖12可知，斷層位錯在土體產生的位錯分布主要集中在上斷線投影前后500 m范圍內，在這個范圍外的土體位移值雖然很大，但是位移分布一致，位移差值很小，隧道結構的破壞主要是由土體位移差值造成的，所以有限元模型的規模取縱向長為500 m。考慮圍巖影響范圍為3～5倍洞徑，模型寬和高分別取為100 m、100 m，Abaqus模型見圖15，坐標系原點為圖中模型左下角點，模型整體位于第一象限。

圖15 Abaqus有限元模型

為了驗證數值模型的準確性，首先對不加隧道結構只有土體的數值模型進行計算，對該模型施加斷層位錯作用產生的土體同震變形場，提取地表沿縱向中線的水平和豎向位移值，與Okada解析解進行對比，結果見圖16。由圖16可知，從計算結果可以看出，本文精細化輸入方法的數值解與使用Okada位錯理論計算得到的解析解相差極小，驗證了本文方法的準確性。

圖16 有限元模型計算結果驗證

下一步加入隧道結構形成土體-隧道有限元模型，隧道頂板覆土厚度為20 m，隧道斷面選取圓形截面，尺寸參數見圖17，模型中各材料的參數見表7。由于本文重點為為計算方法展示，為節約計算時間隧道襯砌結構采用彈性本構，若有更高的計算需求可替換為各類彈塑性本構，本文方法同樣適用。土體采用Abaqus中的實體單元C3D20R模擬，隧道襯砌采用殼單元S4R模擬，土體-隧道有限元模型網格見圖18，模型X軸為隧道橫向寬度方向，Y軸為隧道縱向長度方向，Z軸為土體的豎向高度方向。

圖17 隧道結構斷面尺寸參數(單位：cm)

表7 有限元模型材料參數

圖18 有限元模型網格劃分

為計算斷層錯動而引發的隧道襯砌結構的受力與變形效應，將上述基于分層土體位錯模型計算得到的土體同震變形場作為模型的位移邊界施加在模型上來模擬斷層的錯動作用，得到隧道襯砌結構的破壞發展過程，斷層錯動前后的隧道變形絕對值云圖見圖19。

圖19 斷層錯動前后的隧道變形絕對值

由圖19可知，變形絕對值沿隧道縱向先減小后增大，隧道中部變形值的變化速率最大，為進一步分析隧道襯砌的受力分布情況，選取典型觀測位置見圖20，繪制隧道各處主應力隨隧道縱向距離變化曲線見圖21。

圖20 隧道橫截面觀測位置分布

圖21 各觀測位置主應力隨隧道縱向距離變化

由圖21可知，隧道的拱腰和拱腳各處整體呈受壓狀態，最大壓應力出現在拱腰附近，拱腳處也具有較大的壓應力峰值；而拱頂和仰拱在斷層上盤靠近上斷線投影處呈受拉狀態，這與斷層上盤在該處的大幅抬升密切相關，最大拉應力出現在仰拱附近。由汶川地震的震害調研可知[24,50-51]，隧道的仰拱和拱頂主要出現抗拉破壞，拱腰和拱腳主要出現抗壓破壞，拱肩損傷較小，這與計算結果相符。因此，本文基于分層土體模型提出的精細化斷層錯動輸入方法可為工程中隧道的抗位錯分析提供一定的指導作用。

3 結論

通過理論分析和數值計算，探究了土體層狀構造、斷層破裂面上位錯分布和錯動輸入方法在同震變形場計算中帶來的誤差，揭示了現有隧道抗位錯計算中同震變形輸入方法的不足，提出了精細化的位錯輸入方法，并將該方法應用于成蘭鐵路某隧道的三維抗位錯分析，得出以下結論：

(1)在斷層錯動導致土體同震變形場的計算中，土體分層結構的影響是不可忽略的，地表同震位移差異達40%～60%。尤其是在靠近斷層上斷線地表投影區域，位移差異更為顯著，需要考慮地震發生地區地殼層狀構造的影響。

(2)斷層面上的位錯量的分布模式同樣給同震變形場的計算結果帶來誤差，但量值僅為10%左右，該誤差在震源位錯量反演不明確的條件下可忽略。

(3)錯動輸入方法導致的地表同震變形誤差最大，傳統簡化方法的誤差可達25%～80%，將斷層上盤固定不動僅下盤抬升的輸入方式顯然是不準確的。

(4)將本文方法的數值解與Okada解析解對比，驗證了本文精細化輸入方法的準確性，并以成蘭鐵路某隧道為例分析了其在典型斷層錯動下的破壞狀態，隧道的仰拱和拱頂主要出現抗拉破壞，拱腰和拱腳主要出現抗壓破壞，拱肩損傷較小，可為跨斷層隧道的抗位錯分析提供參考。