淤地壩減輕坡溝系統滑坡侵蝕的數值模擬

鄒兵華，袁 潔，李占斌，李 鵬

（1.中國水電顧問集團成都勘測設計研究院，四川 成都610072；2.西安理工大學 西北水資源與環境生態教育部重點實驗室，陜西 西安710048；3.西安理工大學 水利水電學院，陜西 西安710048）

黃土高原有大小溝壑27萬多條，溝壑縱橫，支離破碎，溝壑密度3～5km／km2，溝壑面積占土地總面積的20%～40%，其支離破碎的地貌主要是由溝蝕造成的，溝蝕的發展使溝壑面積日益擴大，耕地面積日趨縮小［1］，是世界上水土流失最嚴重的地區之一。黃土高原的各類溝壑中以溝頭前進、溝底下切、溝岸擴張3種形式的溝蝕危害最為嚴重［1］。溝蝕加劇了面蝕的發展，造成了更多的陡峭臨空面，加劇了重力侵蝕［2－3］。重力侵蝕作為黃土區主要侵蝕類型，其研究一直處于滯后狀態［4］。理論研究方面，重力侵蝕的微觀及宏觀力學機理尚不明晰，計算模型的構建剛起步；試驗及觀測方面，至今仍未形成重力侵蝕系統觀測和試驗研究體系［4］。根據黃河水利委員會西峰、天水、綏德3個水土保持站在典型小流域的調查，黃土高原溝壑區的西峰南小河溝，重力侵蝕面積占流失面積的9.1%，重力侵蝕量占總流失量的57.5%；黃土丘陵溝壑區第Ⅲ副區的天水呂二溝，重力侵蝕面積占流失面積的30.7%，重力侵蝕量占總流失量的68.0%；黃土丘陵溝壑區第Ⅰ副區的綏德韭園溝，重力侵蝕面積占總流失面積的12.9%，重力侵蝕量占總流失量的20.2%。在黃土塬區和丘陵區，溝頭前進多以土體崩塌形式進行，溝岸擴張是崩塌與滑坡共同作用的結果［5－7］。

淤地壩系作為治理水土流失的主要工程措施，在蓄水攔沙、防 洪 保 收 等 方 面 起 到 了 重 要 作 用［1，5，8－10］，有機統一了當地致富和治河的關系，深受黃土高原地區群眾的喜愛，同時又為治河部門所關注。打壩淤地、蓄水攔沙是流域水土保持綜合治理的一項重要措施，目前國內外有關淤地壩的研究主要集中在壩系相對穩定、減水減沙和生態環境效應等方面［10－11］，而對淤地壩減輕坡溝系統重力侵蝕特別是滑坡侵蝕方面的研究則少有人提及。淤地壩通過抬高溝底侵蝕基準面，提高了溝坡的穩定性，減少了溝坡重力滑坡侵蝕發生的可能性。本研究采用數值模擬方法對隨著壩地逐漸淤高溝道坡溝系統的穩定性、滑塌概率以及坡溝系統滑坡侵蝕破壞的部位等方面進行研究，以期為坡溝系統水土保持工程措施及生物措施的配置提供有益的參考，并為評價坡溝系統的穩定性提供一定可靠度的依據。

1 研究區概況

本研究選取的關地溝是王茂溝流域左岸的一條典型支溝。王茂溝是陜北綏德縣韭園溝中游左岸的一條試驗性治理支溝，是無定河的一條2級支溝，地理坐標為東經110°20′26″—110°22′46″，北緯37°34′13″—37°36′03″，海 拔 高 度 940～1 188m，流 域 面 積5.97km2，主溝長3.75km，溝道平均比降為2.7%，溝谷地面積2.97km2，占流域總面積的46.7%，溝壑密度4.3km／km2。流域內共有溝道46條，其中＜0.1km2的溝道24條，0.1～0.5km2的溝道18條，0.5～1.0km2的溝道2條，＞1.0km2的溝道2條。地面坡度一般在20°以上，坡度＜10°的面積僅占2.2%，坡度＞45°以上的面積占34.7%。地形地貌主要由梁、峁以及分割梁峁的溝谷組成［12］，地形破碎，溝壑縱橫，坡陡溝深。梁頂和峁頂面積不大，均略成穹形，坡度在8°～10°，梁、峁坡的坡度一般在20°～35°。梁、峁坡以下為溝谷，溝谷橫斷面在上游及支溝均呈“V”字型，在下游略呈“U”字型。溝谷坡極陡，一般都在35°以上。按土壤侵蝕區劃，屬黃土高原丘陵溝壑區第Ⅰ副區。

流域內地質構造比較單純，表層多被質地勻細、組織疏松的黃綿土覆蓋，厚度20～30m，在梁、峁頂，梁、峁坡上均有分布。其下為紅色黃土，厚度50～100m，多出露于溝谷。再下為基巖，基巖主要為三迭紀的砂頁巖，巖層傾角甚小，接近水平，其露頭只在干溝中下游溝床及其兩側有所見。

流域屬北溫帶干旱大陸性氣候，年平均氣溫10.2℃，無霜期160d左右，夏季多東南風，春秋多西北風。流域多年平均降水量513.1mm，降雨量的年際變化率大，年最大降雨量是年最小降水量的3.5倍；降雨的年內分配極不均勻，年內降雨量主要集中在汛期7—9這3個月，汛期降水量占年降水量的73.1%，且多以暴雨形式出現，一次暴雨產沙量為全年產沙量的60%以上，土壤侵蝕以水蝕及重力侵蝕為主，治理前流域年平均侵蝕模數為18 000t／（km2·a），屬劇烈侵蝕區，在黃土丘陵溝壑區具有一定的代表性。

2 坡溝系統滑坡侵蝕的數值分析原理與方法

2.1 崩塌滑坡的蒙特卡洛概率計算

對于土質邊坡，根據其土體結構、破壞機理和受力狀況，可以建立坡體地質條件和環境因素的狀態函數：Z＝g（X1，X2，…，Xm）其中 X1，X2，…，Xm為m個具有一定分布、獨立統計的隨機變量，假定它們的統計量已知。如果把狀態函數定義為安全系數，且隨機地從諸隨機變量Xi的全體中抽取同分布的變量X′1，X′2，…，X′m，則由上述狀態函數可求得安全系數的一個隨機樣本。如此重復，直到達到預期精度的充分次數N，則可得到N個相對獨立的安全系數Z1，Z2，…，Zn。安全系數所表征的極限狀態為Z＝1，可構造一個隨機變量R：當Z≤1時，R＝1；當Z＞1時，R＝0。設在N次隨機抽樣的試驗中，出現R＝1即Z≤1的次數為M，則坡體破壞的概率pf為［13－14］：

式中：pf——坡體破壞的概率；M——出現安全系數Z≤1的次數；N ——隨機抽樣試驗的次數。

此式即為直接蒙特卡洛法計算坡體破壞概率的公式。顯然當N足夠大時，由安全系數的統計樣本Z1，Z2，…，Zn可比較精確地近似安全系數的分布函數G（z），并估計其分布參數，其均值和標準差分別為：

式中：zi——安全系數的樣本值；μz——安全系數的均值；σz——安全系數的標準差。

進而可根據G（z）擬合的理論分布，通過積分方法求得破壞概率。本文中N取值為10萬次，滑坡計算選用Spencer法。

2.2 有限元強度系數折減法的基本原理

有限元強度系數折減法的基本原理是將坡體強度參數黏聚力c和內摩擦角φ，同時除以一個折減系數ω，得到一組新的c′，φ′，其中：

式中：c——巖土體原本的黏聚力；φ——巖土體原本的內摩擦角；ω——巖土體強度折減系數；c′——巖土體折減后的黏聚力；φ′——巖土體折減后的內摩擦角。

然后作為新的資料參數輸入，再進行試算，利用相應的穩定判斷準則，程序可以自動根據彈塑性計算結果得到破壞滑動面，確定相應的ω值為坡體的最小穩定安全系數，此時坡體達到極限狀態，發生剪切破壞，同時又可得到坡體的破壞滑動面［15－17］。有關研究［15－19］表明，有限元強度折減法的安全系數在本質上與傳統方法是一致的。鄭穎人等［15］通過多種比較計算說明有限元折系數法用于分析土坡穩定問題是可行的，但必須合理地選用屈服條件以及嚴格地控制有限元法的計精度。

2.3 強度系數折減法失穩判斷準則與屈服準則

采用有限元強度折減法分析邊坡穩定性的一個關鍵問題，是如何根據有限元計算結果來判別邊坡是否處于破壞狀態。目前的失穩判據主要有兩類：第一類是在有限元計算過程中采用力和位移的不收斂作為邊坡失穩的標志［18］；第二類以廣義塑性應變或等效塑性應變從坡腳到坡頂貫通作為壩坡破壞的標志［19］，以上兩種判據得到的安全系數相差不大［15］。本文采用節點不平衡力的不收斂作為判據，便于計算和可視化，同時根據塑性區的范圍及其連通狀態確定相應的安全系數，以此評價坡體的穩定性。當前流行的有限元軟件中采用的均是廣義米賽斯準則［15］，由于強度折減有限元法計算的安全系數與選用的屈服準則密切相關，鄭穎人等［15，18］等研究發現不等角六邊形外接圓（DP1）屈服準則計算的安全系數較傳統的莫爾—庫侖準則會偏大許多，而鄭穎人［15］研究提出的莫爾—庫侖等面積圓準則的計算精度與Spencer的平均誤差在5%左右，可滿足工程要求，并給出了不等角六邊形外接圓（DP1）屈服準則求得的安全系數與莫爾—庫侖等面積圓準則求得的安全系數之間的轉換關系式（5）。本文把在外接圓（DP1）屈服準則下求得的安全系數轉換為莫爾—庫侖等面積圓屈服準則下的安全系數：

式中：η——外接圓（DP1）屈服準則下安全系數轉換為莫爾—庫侖等面積圓屈服準則下安全系數的轉換系數；φ——巖土體的內摩擦角。

3 坡溝系統滑坡侵蝕的有限元數值模擬與分析

3.1 坡溝系統概化模型及有限元計算模型

根據相關地質資料，關地溝地層構造主要是馬蘭黃土（Qeol3），梁、峁頂、峁坡均有分布，厚度20～30m，其下為離石黃土（Qeol2），厚度50～100m，多出露于溝坡上，再下主要是三疊紀砂頁巖層，基本接近水平，多出露于干溝、支溝的下游及其兩側。通過實測數據統計分析得知，坡溝系統峁坡坡度在21°～25°所占比例最大，達到50%，而實測坡度在21°～25°的平均值為23.1°，取峁坡坡度24°；溝坡坡度在36°～45°所占比例最大，達到43.6%，而實測坡度在36°～45°范圍內的平均值為40.8°，取溝坡坡度41°。通過對關地溝基本坡溝地貌特征的綜合分析，建立了坡溝系統概化模型。該重力滑坡侵蝕模型取峁坡坡度24°，溝坡坡度為41°。溝坡概化模型土層從上到下分別為馬蘭黃土（Qeol3）和離石黃土（Qeol2），厚度分別取為土層厚度的平均值25m和75m。

關地溝流域典型坡溝系統的有限元網格見圖1。坡溝系統的計算采用平面應變條件下摩爾—庫侖不等角六邊形外接圓屈服準則（DP1）和非關聯流動的理想彈塑性模型［20－21］。坡溝系統二維模型基巖底邊界為固定約束，左邊界為水平約束，其他自由。彈塑性有限元計算采用大變形靜態分析選項，最大迭代次數為1 000。

圖1 坡溝系統三維有限元網格

3.2 計算參數取值

根據邊坡工程經驗［13－14］、現場資料分析、現場及室內巖土物理力學試驗和有限元計算的需要，溝坡概化模型各土層材料特征取值如表1所示。

表1 坡溝系統的計算參數

3.3 溝坡系統穩定性隨壩地淤高度的變化

由于關地溝4號壩從1959年修建到1987年被洪水沖毀的28a運行過程中，累計於高4.86m，年均淤高0.18m；垮壩后由于群眾對壩體的零星修復填筑及2005年的修復加高，壩地被填高，壩地已不再是原來垮壩前的自然淤高，截至2007年壩地累計高10.8m。由于關地溝4號淤地壩跨壩前壩地的年均淤高僅為0.18m，在跨壩前的淤高模擬計算中以5a為1個時間跨度，計算壩地每淤高0.9m時的坡溝系統穩定特征參數；垮壩后為了與垮壩前保持一致性，壩地以0.9m遞增。在有限元程序中，模擬淤地壩壩地按0～0.9，0.9～1.8，1.8～2.7，2.7～3.6，3.6～4.5，4.5～4.86，4.86～5.4，5.4～6.3，6.3～6.48，6.48～7.2，7.2～8.1，8.1～9.0，9.0～9.9，9.9～10.8m的淤積過程每淤積（增高）0.9m其相應的穩定系數、滑塌概率和滑塌量的散點圖如圖2—4所示。

圖2 隨壩地淤高坡溝系統穩定系數的變化規律

圖3 隨壩地淤高坡溝系統滑坡侵蝕滑塌概率的變化規律

圖4 隨壩地淤高坡溝系統滑坡侵蝕滑塌量的變化規律

由圖2—4分別擬合出相關系數最高的隨壩地淤高坡溝系統的穩定系數、滑塌概率和滑塌量變化規律的相關方程如表2所示。

表2 隨壩地淤高坡溝系統的穩定系數、滑塌概率和滑塌量的變化規律

當壩地分別淤高到11.7，12.6，3.5，14.4和15.3m時，由程序計算得到坡溝系統的穩定系數、滑塌概率和滑塌量與通過上述擬合出的相關方程計算得到的對應值的比較如表3所示。

由表2—3可以得出，隨壩地淤高坡溝系統的穩定系數和滑塌概率滿足二項式分布規律，而滑塌量則滿足線性減少分布規律。穩定系數和滑塌概率的相關系數最高，達到0.980；滑塌量的相關系數為0.977，也相當的高。通過壩地累計淤高11.7，12.6，13.5，14.4和15.3m時坡溝系統的穩定系數、滑塌概率和滑塌量的程序計算值和擬合方程計算值的比較，穩定系數的最大相對誤差為1.61%，可見穩定系數的變化是很有規律的，而且穩定性好；滑塌概率的相對誤差較大，最大達到15.67%。這是由于滑塌概率本身就是另一個概率計算的結果，還有別的影響因素。滑塌量的最小誤差7.92%，最大誤差20.09%，這可能是由于隨著壩地的淤高，坡溝系統形態的改變，滑坡侵蝕發生過程中最危險和最可能滑動面變化，導致滑坡侵蝕中量的波動比較大，故而滑塌量的相對誤差較大的原因。總體來說，三者的相對誤差均在滿足的范圍內，擬合出的方程精度較高，可應用于該溝道重力滑坡侵蝕的定量分析和計算。

表3 坡溝系統的穩定系數、滑塌概率和滑塌量軟件計算值和相關方程計算值的比較

3.4 穩定系數最小、坡溝系統瀕臨破壞時的應力和位移分布

由圖2可得，在壩地相對原始溝道累計淤高0m的時候，其穩定系數最小，滑塌的概率最大。為了優化配置坡溝系統的水土保持工程措施和生物措施，達到水土保持措施的最佳防治效果，從力學機理上對坡溝系統的應力場和位移場進行分析，得出坡溝系統應力場和位移場的最危險區域，有針對性地進行措施強化，保證坡溝系統的土壤侵蝕減到最小。由于坡溝系統的近似對稱性，這里只考慮溝道的1／2來進行有限元分析。關地溝4號壩上游典型坡溝系統的有限元網格見圖1。由有限元強度折減法得摩爾—庫侖不等角六邊形外接圓屈服準則（DP1）下求得此時的安全系數為1.507，將其按式（5）轉換為莫爾—庫侖等面積圓屈服準則下的安全系數為1.224，其與滑坡分析采用Spencer法求得此時的安全系數1.255的相對誤差為2.42%，在鄭穎人等［16，18］研究指出莫爾—庫侖等面積圓屈服準則的安全系數與Spencer法的平均誤差5%范圍內，精度滿足要求。

由計算可得，坡溝系統中X方向的位移均為正，在溝坡中下部位垂直向坡體內約15m的范圍為X方向位移的最大值區域，X方向的其他位移從右向左、從溝坡向峁坡遞減，直至減小為最小值0。Y方向均為下沉位移，最大值在從峁頂向左向右各延伸10m左右的扇形區域，其他下沉位移從峁頂向溝坡坡腳遞減，直至減小為最小值0。

坡溝系統中除峁坡和溝坡坡緣線向坡體內一定范圍內出現了局部拉應力外第1主應力和第3主應力幾乎全為壓應力，且第1主應力的拉應力區域幾乎包括在第3主應力拉應力區域內。由于土體是抗壓不抗拉的，出現拉應力就會使土體發生破壞，進而發生侵蝕，這就從力學機理上說明了坡溝兼治的正確性和峁坡和溝坡進行生物措施和工程措施治理的必要性。坡溝系統中第1主應力的最大拉應力值區域位于從峁頂到溝坡坡腳的整個坡緣線在峁坡向坡體內延伸6m左右、在溝坡向坡體延伸3m左右的帶形區域內；第3主應力的最大拉應力值在從峁頂經峁坡向溝坡坡緣線向坡體內延伸9m左右的帶形區域內。

4 結論

基于邊坡分析軟件和有限元技術，通過建立關地溝典型坡溝系統的概化模型，得到并驗證了隨淤地壩壩地淤高坡溝系統的穩定系數和滑塌概率滿足二項式分布規律，而滑塌量則滿足線性減少分布規律，三者的精度均較高，可用于該溝道重力滑坡侵蝕的預報、滑坡侵蝕的定量計算和指導淤地壩的分期加高。

應用有限元強度折減法對坡溝系統的應力場和位移場進行了分析，指出了穩定系數最小、滑坡概率最大、坡溝系統瀕臨滑坡侵蝕時的最大應力和位移分布區域。（1）X方向的最大位移：溝坡中下部位垂直向坡體內約15m的范圍；（2）Y方向的最大位移：從峁頂向左向右各延伸10m左右的扇形區域；（3）拉應力最大值區域：從峁頂經峁坡向溝坡坡緣線向坡體內延伸9m左右的帶形區域。這些區域是坡溝系統最容易發生侵蝕的危險區域，在配置坡面水土保持工程措施（水平溝和魚鱗坑等）和生物措施的時候要有針對性地進行重點防護，以使坡溝系統的土壤侵蝕降低到最低限度。

［1］唐克麗.中國水土保持［M］.北京：科學出版社，2004.

［2］Williams J R，Berndt H D.Sediment yield prediction based on watershed hydrology［J］.Transaction of the ASAE，1977，20（6）：1100－1104.

［3］Nachtergaele J，Poesen J，Vandekerck H L，et al.Testing the ephemeral gully erosion model（CEGEM）for two mediter ranean environment［J］.Earth Surface Process and Landforms，2001，26（1）：17－30.

［4］薛海，王光謙，李鐵鍵.黃河中游區重力侵蝕研究綜述［J］.水科學進展，2009，7，20（4）：599－606.

［5］崔靈周，李占斌，朱永清，等.流域侵蝕強度空間分異及動態變化模擬研究［J］.農業工程學報，2006，22（12）：17－22.

［6］高鵬，劉作新，鄒桂霞.丘陵半干旱區小流域土地資源定量化評價研究［J］.農業工程學報，2003，19（6）：298－301.

［7］于國強，李占斌，李鵬，等.黃土高原小流域重力侵蝕數值模擬［J］.農業工程學報，2009，25（12）：74－79.

［8］陳浩.黃土丘陵溝壑區流域系統侵蝕與產沙關系［J］.地理學報，2000，55（3）：354－363.

［9］郭文元，賈志軍.晉西淤地壩試驗研究文集［M］.鄭州：黃河水利出版社，2001：102－165.

［10］鄭寶明.黃土高原丘陵溝壑區淤地壩建設效益與存在問題［J］.水土保持通報，2003，23（6）：32－35.

［11］黃自強.黃土高原地區淤地壩建設的地位及發展思路［J］.中國水利，2003（17）：8－11.

［12］鄭寶明，王曉，田永宏，等.淤地壩實驗研究與實踐［M］.鄭州：黃河水利出版社，2003：152－367.

［13］王家臣.邊坡工程隨機分析原理［M］.北京：煤炭工業出版社，1996：68－139.

［14］祝玉學.邊坡可靠性分析［M］.北京：冶金工業出版社，1993：45－96.

［15］鄭穎人，趙尚毅.有限元強度折減法在土坡與巖坡中的應用［J］.巖石力學與工程學報，2004，23（19）：3381－3388.

［16］Zienkiewicz O C，Taylor R L.The finite element method［M］.New York：McGraw－Hill，1989.

［17］Zienkiewicz O C，Humpeson C，Lewis R W.Associated and nonassociated visco－plasticity in soil mechanics［J］.Geotechnique，1975，25（4）：671－689.

［18］趙尚毅，鄭穎人，時衛民，等.用有限元強度折減法求邊坡穩定安全系數［J］.巖土工程學報，2002，24（3）：343－346.

［19］吳海真，顧沖時.有限元強度折減法在土石壩邊坡穩定分析中的應用［J］.水電能源科學，2006，8（4）：54－58.

［20］石長.趙新銘.用有限元強度折減法分析邊坡穩定［J］.隧道建設，2006，26（3）：5－8.

［21］張彩雙，李俊杰，胡軍.有限元強度折減法的邊坡穩定分析［J］.中國農村水利水電，2006（5）：72－74.