基于正交試驗的巖質邊坡動力穩定性分析

付宏淵，劉建華，張立，呂東濱

(1.長沙理工大學 交通運輸工程學院，湖南，長沙，410004；2.湖南華罡規劃設計研究院有限公司，湖南，長沙，410076)

我國是一個多山的國家，特別是西南山區存在大量的自然巖質邊坡。目前，該區域地震十分活躍和頻繁，正處于地震活動高峰時期，地震誘發的巖體邊坡破壞(滑坡與崩塌)是常見的地震災害[1?3]，在山區和丘陵地帶，地震誘發的邊坡滑動和坍塌往往分布廣、數量多、危害大，成為影響該地區公路建設，特別是確保震后生命線工程暢通所需要迫切解決的核心問題[4?5]。地震的波動震蕩造成的邊坡滑坍往往具有啟程劇動、行程高速等特點，從而使得地震誘發的邊坡失穩無論是在規模，還是在造成的經濟財產損失上都較其他因素引起的邊坡失穩影響更為巨大。如2008年汶川MS8.0級大地震所觸發的滑坡、崩塌、碎屑流等總數達3～5萬處，其中對震后人員安全和臨時安置構成直接威脅的災害隱患點就達12 600余處(四川省39個重災縣統計)，規模大于107m3的巨型滑坡達數十處[6]。汶川大地震雖然給我們造成了巨大的災難，但同時也給我們留下了豐富的震害信息，地震觸發的大量滑坡、崩塌等地質災害為我們研究邊坡動力荷載作用下的穩定性提供了寶貴的資源，進而為可能發生逆斷層型地震地區的地震危險性分析、地震動衰減分及沿線公路建設的邊坡防護等工作提供了依據。在此，本文作者在總結分析震區巖質邊坡失穩影響因素的基礎上，采用FLAC3D方法對湖南省對口援建項目即四川理縣至小金公路工程中豹子嘴邊坡建立了三維模型，并采用2008年“5·12”汶川M8.0級地震中木卡站實測的地震波，在對地震波進行校正的基礎上，對該邊坡進行動力數值模擬和穩定性分析，最后采用正交試驗設計，對影響該邊坡動力穩定性主要因素進行敏感性分析。

1 邊坡動力穩定性與地形地貌的關系

1.1 水系的影響

地震觸發的滑坡、崩塌等地質災害在區域分布上的一個顯著特征是沿河流水系呈線狀分布。這些水系河流與發震斷裂一起控制了地震災害的區域性分布。其中，發震斷層的控制更具宏觀性，而河流水系的控制則更具有部位的具體性。一般而言，這些河流水系所處的地貌構造單元多為中低山?高山地區，高差一般為500～3 500 m，由于受到構造運動和斷裂活動的控制，導致這些地區河谷深切，山體陡峻，形成眾多的河曲和斷層崖，客觀上形成了地表地質災害易發的條件。此外，在河流凹岸坡腳處由于受到河流垂向和側向的侵蝕沖刷，巖體應力釋放，在地震荷載作用下，更加有利于產生大規模的滑坡和崩塌。再者，交通線常常沿著河谷進行展布，造成了大量的人工邊坡，破壞了邊坡的自身穩定性，這也是造成大量地質災害沿著河谷發育的重要原因。因此，汶川地震中觸發的滑坡、崩塌等地質災害具有顯著的沿河流水系呈線狀分布的特征。

1.2 坡形的影響

邊坡的坡形是影響邊坡動力穩定性的最直觀的因素。自然界中邊坡坡面的形狀千變萬化，大體可分為直線坡、凸坡和凹坡3種，其他的形狀都是由這3種組合變化而來。研究表明：(1)在 3種坡形中，除凹型坡外，地震作用下最大剪應變和最大位移均出現在離坡腳 1/8～1/5高程處，失穩均在臨近坡腳的地方先出現較大剪應變破壞之后引起的淺層滑動；(2)在較小地震作用下，凹型最利于邊坡穩定，而凸型最不利。但在高烈度地震作用下(Ⅷ以上)，在凹坡變坡點稍上位置會出現較大的變形[7?8]。

通過對都江堰到汶川一帶的實地調研并結合其他學者對汶川地震震后地質災害的調查資料也得出了相似的結論：在地震動的作用下，不同坡形的邊坡具有不同的動力響應，各自的動力穩定性狀況也大為不同，直線型的斜坡發生滑坡和崩塌的概率相對較小，凹坡和凸坡更容易發生滑坡和崩塌，尤以凹坡上發生滑坡和崩塌的概率最高，并且大多數是發生在邊坡坡度由緩變陡的過渡轉折部位。這與在靜力作用下邊坡的穩定性有很大的區別。在靜力作用下凹坡的穩定性大于平坡的穩定性，平坡的穩定性又大于凸坡的穩定性。

1.3 坡度的影響

一般來說，地震動作用下邊坡地形坡度越陡則越容易引發滑坡、崩塌。就地震滑坡而言，在某一區域條件下通常存在著一個觸發滑坡的優勢角度范圍，超過這個范圍一般會發生崩塌。但對于不同的區域條件下，觸發滑坡的優勢角度范圍也不盡相同，表現出明顯的區域差異性。2008年的汶川地震，絕大部分的地質災害點(汶川地震后16個極重災縣或重災縣范圍內)都分布在20°～50°的范圍內，數量占總數的86.8%，其中以 40°～50°范圍分布密度最大；同時，從坡面上單位塊體物質的切應力(τ)、摩擦因數(f)與坡度間的關系圖(見圖1)也可以得出：易發生滑坡的斜坡坡度多數集中在 25°～45°[9?10]。

1.4 發育位置的影響

地震動反應幅值和頻譜隨著地形高度而變化，這種現象在很多已有的強震觀測和物理實驗的結果中得到了驗證。高野秀夫斜坡地震效應的觀測結果[11]表明：(1)斜坡上的地震烈度相對于谷底大約增加1°；(2)在角度超過 15°的圓錐狀山體上部點的位移幅值與下部點的位移幅值相比，其局部譜段值增加高達 7倍；王存玉[12]通過振動模型實驗發現：邊坡頂部對振動的反應幅值與邊坡底部相比存在明顯的放大現象(即存在垂直向放大效應)，邊坡的邊緣部位對振動的反應幅值與內部(處于同一高度上的 2點)相比也存在放大現象(即存在水平向放大效應)；雖然現有的研究還不能定量的給出地震動響應隨著邊坡高度變化的范圍，但是至少可以得出這樣一個趨勢即對于同一邊坡而言，其上部的地震動反應幅值相對于下部的大，存在放大效應。這就表明了在地震動作用下邊坡的上部相對于下部應該更容易發生破壞，滑移災害破壞集中在邊坡的中上部發育。這種現象在汶川地震造成的滑坡、崩塌災害中表現也很顯著。

圖1 邊坡坡度與切應力及摩擦因數的關系Fig.1 Relationship among slope angle with shear stress and friction coefficient

通過對汶川地震后汶川沿線公路的317個滑坡、崩塌的開始滑動的位置進行了統計和分析，結果見表1?？梢娺吰略诘卣饎幼饔孟掳l生的滑坡、崩塌災害大多數發育集中在邊坡的中上部。

表1 滑坡與發生位置的關系Table1 Relationship between landslide and occurrence position

2 巖質邊坡動力數值分析

2.1 地震波的校正

本文依托工程項目為湖南省對口援建項目——理縣至小金公路，其地形起伏大，沿線路地面相對高差較大，地貌類型為高山沖溝地貌，成因為斷裂構造引起的斷塊山，山頂呈三角形狀，地形坡度一般為30°～70°。地質構造屬于四川西部地槽區金湯弧形褶皺東翼的薛城—臥龍“S”型構造，地質構造復雜，該構造體系東臨茂汶斷層和九頂山華夏系構造體系斜接，褶皺在這一帶表現特別緊密。根據《中國地震動峰值加速度區劃圖》(GB 18306—2001)和國家質量技術監督局《中國地震動參數區劃圖》(2001)，以及“5·12”汶川大地震分析，邊坡所處區域屬強震區，地震動峰值加速度為0.15g，反應譜特征周期為0.40 s。

本文采用的地震波為汶川MS8.0級地震主震臺站理縣木卡站實測的離散的未校正的加速度記錄，在施加之前有必要進行基線校正和濾波處理，其基線校正的目的是在動力結束時模型底部不會出現繼續的速度和殘余的位移。通過對理縣木卡站實測的加速度時程進行了低通濾波處理，在加速度時程上增加一個低頻率的波形(本文采用三次多項式)，過濾了地震波中頻率大于15 Hz的成分，使得最終位移和速度均為零，從而實現了加速度時程的基線校正[9]，調整后的傅里葉振幅譜見圖2。

圖2 校正后的加速度的傅里葉震幅譜Fig.2 Fourier amplitude spectrum of correcting acceleration

FLAC3D可以在模型邊界或內部節點施加動荷載來模擬材料受到外部或內部動力作用下的響應，允許的動力載荷輸入可以是加速度時程、速度時程、應力(壓力)時程和集中力時程。本文采用的是黏滯邊界，通過在模型的法向和切向分別設置自由的阻尼器，從而實現吸收入射波，阻尼器提供的法向黏滯力和切向黏滯力的計算公式分別為：

由于黏滯邊界上的作用力是根據邊界上的速度分量計算得到的，如果再施加速度載荷就會使黏滯邊界失效。因此，在黏滯邊界上只能輸入應力時程[13]。首先通過數值積分把加速度時程(見圖3)轉化成速度時程(見圖4)，然后，利用式(3)和式(4)將速度時程轉化成可以直接輸入的應力時程。

式中：σn和σs分別為施加在黏滯邊界上的法向和切向應力；vn為模型邊界上法向速度分量；vs為模型邊界上切向速度分量；ρ為介質密度；CP為P波波速；CS為S波的波速。

2.2 動力數值模型

分析的邊坡為理縣至小金公路K2+640—K2+824段豹子嘴邊坡，位于梭羅河左岸一級階地，階地寬10～20 m，邊坡位置河岸高出河床約15.00 m，地基巖層主要為第四系更新統塊石、滾石和三疊系中風化板巖，下覆巖層為反傾向巖層，對邊坡穩定性有利，各主要巖土層物理力學參數取值如表2所示。

表2 計算模型中各主要巖土層物理力學參數Table2 Physical mechanical parameters of rock layers in calculation model

FLAC3D計算模型范圍，順河方向(X軸方向)長184 m，即順河方向包括整個邊坡的長度；垂直河谷方向(Y軸方向，指向邊坡內為負)寬度為100 m，模型最大高度(Z軸方向，垂直向上為正)110 m。邊坡的三維模型如圖5所示。

圖3 校正后的加速度時程Fig.3 Correcting acceleration time-history

圖4 速度時程Fig.4 Velocity time-history

圖5 邊坡三維計算模型Fig.5 Three-dimensional calculation pattern of slope

根據邊坡的工程地質條件、地貌特征等因素，其主計算坡面的尺寸如圖6所示。

圖6 計算主剖面Fig.6 Calculation of main profile

2.3 巖質邊坡數值模擬結果分析

圖7和圖8所示分別為邊坡在地震作用結束后的最小大主應力和最小主應力云圖。從圖7和圖8可以看出：邊坡在地震荷載作用下并未出現明顯的拉應力區，基本上以壓應力為主，即假如邊坡發生破壞，將以“壓?剪”破壞模式為主。

在坡面附近主應力等值線與坡面近于平行且較為平滑，而在邊坡坡腳處則出現了應力集中帶，形成了一個最大剪應力增高帶。此外，由于巖體分界面的存在，使得其附近區域的巖土體不均勻和不連續，造成分界面附近出現應力阻滯現象，等值線出現突變，發生應力集中。由于邊坡內部分布的應力并不大，因此，不足以降低邊坡強度。

圖7 震后邊坡最小主應力云圖Fig.7 Minimum principal stress contour of slope after earthquake

圖8 震后邊坡最大主應力云圖Fig.8 Maximum principal stress contour of slope after earthquake

3 動力穩定性影響因素敏感性分析

3.1 正交試驗設計原理

在實際問題中考慮1個因素或2個因素對試驗結果的顯著性分析可以選用一元或二元方差分析，而巖質邊坡動力穩定性需考慮多個因素對其受力與變形的影響，可采用正交試驗的方法進行分析。

根據給定需要考察的因素及各因素的水平，選擇與之相適應的正交表Ln(r1×r2×…×rm)。式中：L為正交表；n為正交表行數(即可安排n次試驗)，而m為該正交列數(即試驗最多可安排的因子數)，且第j個 因素有rj個水平。常用的是等水平正交表，即：r1=r2=…=rm=r，簡記為：Ln(rm)。

設A,B, …表示不同的因素；r為各因素的水平數；Ai表示因素A的第i個水平(i=1, 2, …,r)；Xij表示因素j的第i水平的值(i=1, 2, …,r；j=A,B, …)。

在Xij下進行試驗得到因素j第i水平的試驗結果指標Yij，Yij是服從正態分布的隨機變量。在Xij下進行n次試驗可得到n個試驗結果Yijk(k=1, 2, …,n)。有關計算參數如下：

評價因素顯著性的參數為極差Rj，公式為：

極差越大，說明該因素的水平改變對試驗結果影響也越大，極差最大的因素也就是最主要的因素。極差越小的因素雖然不能說是不重要的因素，但至少可以肯定當該因素在所選用的范圍內變化時，對該指標影響不大[14?15]。

3.2 動力穩定性正交試驗設計

本次試驗選的主要因素有：黏聚力、摩擦角、坡比、水平地震加速度和豎直地震加速度。每個因素的取值范圍按一般邊坡工程和工程類比確定，并將其概化為4個水平。在這里，為突出主要矛盾，對其他因素進行適當的簡化。地震區域巖質邊坡的永久位移是邊坡穩定性分析與抗震設計中較合理的參量。尤其是對于數值計算而言，永久位移是對地震作用下數值模型的網格中各單元永久應變進行積分而得到的，計算結果會更加科學，故選擇邊坡表面水平永久位移的最大值作為邊坡在地震作用下穩定性的評價指標，即本次正交試驗的指標。

確定的因素水平表如表3所示。對于5因素4水平正交試驗，最少試驗次數為16次，記為L16(45)。

表3 巖質邊坡動力穩定性正交試驗設計Table3 Orthogonal test design for rock slope dynamic stability

3.3 敏感度分析

選用5因素4水平正交試驗表L16(45)，并按表3所確定的4水平因素，采用FLAC3D程序對該段邊坡進行非線性動力分析，各水平正交因素水平計算結果見表4。

表4 正交試驗計算結果Table4 Orthogonal test calculation result

根據 Wilson等[16]在洛杉磯地區的滑坡災害評估中，巖質邊坡的臨界位移為10 cm，而則該處巖質邊坡在強震動力作用下最大永久位移為4.2 cm，總體上仍為較穩定邊坡。

對表4中的試驗結果進行極差分析，如表5所示。根據極差由大至小順序可排出因素的主次為：坡比、水平地震加速度、豎直地震加速度、摩擦角、黏聚力。

表5 多因素敏感性分析Table5 Sensitivity analysis for multi-variations

4 結論

(1)地震引起的巖質邊坡失穩與地形地貌關系密切，其滑坡、崩塌等地質災害在區域分布上沿河流水系呈線狀分布；直線型的斜坡發生滑坡和崩塌的概率相對較小，凹坡和凸坡更容易發生滑坡和崩塌，尤以凹坡上發生滑坡和崩塌的概率最高，并且大多數是發生在邊坡坡度由緩變陡的過渡轉折部位；滑坡多發生在坡度在 20°～50°范圍內的邊坡上；滑坡、崩塌開始滑動的位置大多發育在邊坡的中上部。

(2)所計算的巖質邊坡在地震荷載作用下并未出現明顯的拉應力區，主要以“壓?剪”破壞模式為主。在坡面附近主應力等值線與坡面近于平行且較為平滑，而在邊坡坡腳處則出現了應力集中帶，形成了一個最大剪應力增高帶。

(3)采用正交試驗設計，對影響該段邊坡動力穩定性的主要因素進行正交試驗設計，得到黏聚力、摩擦角、坡比、水平地震加速度和豎直地震加速度等影響因素對巖質邊坡動力穩定性的敏感性由大到小的順序為：坡比、水平地震加速度、豎直地震加速度、摩擦角、黏聚力。

[1]趙明華, 劉建華, 陳炳初, 等.邊坡變形及失穩的變權重組合預測模型[J].巖土力學, 2007, 28(增刊): 553?557.ZHAO Ming-hua, LIU Jian-hua, CHEN Bing-chu, et al.Variable-weight combination forecasting model of slope deformation and instability[J].Rock and Soil Mechanics, 2007,28(s1): 553?557.

[2]黃理興, 陳奕柏.我國巖石動力學研究狀況與發展[J].巖石力學與工程學報, 2003, 22(11): 1881?1886.HUANG Li-xing, CHEN Yi-bo.Rock dynamics in China: Past,present and future[J].Chinese Journal of Rock Mechanics and Engineering, 2003, 22(11): 1881?1886.

[3]李海波, 蔣會軍, 趙堅, 等.動荷載作用下巖體工程安全的幾個問題[J].巖石力學與工程學報, 2003, 22(11): 1887?1891.LI Hai-bo, JIANG Hui-jun, ZHAO Jian, et al.Some problems about safety analysis of rock engineering under dynamic load[J].Chinese Journal of Rock Mechanics and Engineering, 2003,22(11): 1887?1891.

[4]Shou K J, WANG C F.Analysis of the Chiufengershan landslide triggered by the 1999 Chi-Chi earthquake in Taiwan[J].Engineering Geology, 2003, 68: 237?250

[5]Keefer D V.Landslides caused by earthquakes[J].Geological Society of America Bulletin, l984, 95(4): 406?421.

[6]黃潤秋.汶川 8.0級地震觸發崩滑災害機制及其地質力學模式[J].巖石力學與工程學報, 2009, 28(6): 1239?1249.HUANG Run-qiu.Mechanism and geomechanical modes of landslide hazards triggered by Wenchuan 8.0 earthquake[J].Chinese Journal of Rock Mechanics and Engineering, 2009,28(6): 1239?1249.

[7]Richter C F.Elementary seismology[M].San Francisco: W.H.Freeman and Co, 1958: 23?67.

[8]Keefer D V.Landslides caused by earthquakes[J].Geological Society of America Bulletin, 1984, 95(4): 406?421.

[9]Rodriguez C E, Bommer J J, Chandler R J.Earthquake-induced landslides: 1980—1997[J].Soil Dynamics and Earthquake Engineering, 1999, 18: 325?346.

[10]李鐵漢, 潘別桐.巖體力學[M].北京: 地質出版社, 1980:47?64.LI Tie-han, PAN Bie-tong.Rock mechanics[M].Beijing:Geology Press, 1980: 47?64.

[11]Gazetas G.Seismic response of earth dams: Some recent developments[J].Soil Dyn Earthq Engrg, 1987, 6(1): 3?47.

[12]王存玉.地震條件下二灘水庫岸坡穩定性研究[C]//巖體工程地質力學問題(七).北京: 科學出版社, 1987: 67?81.WANG Cun-yu.Study of slope stability for Ertan reservoir under earthquake condition[C]//Engineering and Geological Mechanical Question for Rock (Seven).Beijing: Science Press,1987: 67?81.

[13]Itasca Consulting Group Inc.FLAC-3D (Fast Lagrangian analysis of continua in 3 dimensions), version 2.00, users manual(Vol Ⅴ)[Z].Minneapolis, USA: Itasca Consulting Group Inc,1997: 101?128.

[14]方開泰, 馬長興.正交設計與均勻設計[M].北京: 科學出版社, 2001: 1?17.FANG Kai-tai, MA Chang-xing.Orthogonal and homogeneous design[M].Beijing: Science Press, 2001: 1?17.

[15]陳希孺.數理統計引論[M].北京: 科學出版社, 1997: 5?22.CHEN Xi-ru.Statistics introduction[M].Beijing: Science Press,1997: 5?22.

[16]Wilson R C, Keefer D K.A real limits of earthquake-induced landsliding[C]//Evaluating Earthquake Hazards in the Los Angeles Region.Washington: US GPO Edition Press, 1985:316?345.