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

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

黃土高原溝壑縱橫，支離破碎，有大小溝壑27萬多條，溝壑面積占土地總面積的20%～40%，是世界上水土流失最嚴(yán)重的地區(qū)之一，其支離破碎的地貌主要由溝蝕造成［1－2］。黃土高原的各類溝壑中以溝頭前進、溝底下切、溝岸擴張三種形式的溝蝕危害最為嚴(yán)重［1－2］。溝蝕造成了更多的陡峭臨空面，加劇了重力侵蝕［2－4］。重力侵蝕是指地表物質(zhì)在重力作用下，分散或成塊，急速或緩慢向下移動的過程，主要表現(xiàn)形式有土的蠕動、松散物的滑動、松散物的流動和崩塌等［5］。重力侵蝕直接或間接地向河道輸送了大量的泥沙，是流域土壤侵蝕和河流泥沙來源的主要物理過程之一，并在區(qū)域地貌演化中扮演了十分活躍的角色［6－9］。根據(jù)黃河水利委員會的西峰、天水、綏德3個水保站在典型小流域的調(diào)查，這些地區(qū)的重力侵蝕面積和侵蝕量均占很大比例，重力侵蝕十分嚴(yán)重［2，8－9］。1986年以后，由于溝道打壩抬高了侵蝕基準(zhǔn)面，穩(wěn)定了溝壁，加之溝谷大多種植了林草，據(jù)實地觀測沒有上述情況出現(xiàn)，溝壁的擴張得到制止［10］。在黃土塬區(qū)和丘陵區(qū)，溝頭前進多以土體崩塌形式進行，溝岸擴張是崩塌與滑 坡 共 同 作 用 的 結(jié) 果［2，8－9，11］。 淤 地 壩控制重力侵蝕的作用是在溝道建壩以后開始的，其減蝕量一般與溝壑密度、溝道比降以及溝谷侵蝕模數(shù)等因素有關(guān)［12］。溝道里修建淤地壩以后，隨著壩前泥沙的淤積，抬高了侵蝕基準(zhǔn)面，可防止溝道下切和溝岸坍塌，減少溝道侵蝕；壩內(nèi)淤積厚度的不同，對重力侵蝕的控制程度不同，且壩內(nèi)泥沙的淤積厚度越大，對重力侵蝕發(fā)生的控制越為有利［13］。

坡溝系統(tǒng)是黃土高原小流域特有的結(jié)構(gòu)。黃土高原的坡溝系統(tǒng)是從峁坡坡頂經(jīng)坡面、溝坡至溝道但不包括溝道的一個統(tǒng)一體，既是區(qū)域侵蝕產(chǎn)沙的主要源地，又是控制水土流失，減少下游河道淤積，恢復(fù)與重建生態(tài)環(huán)境的基本治理單元，歷來為黃土高原小流域治理中理論性與實踐性均很強的重要科學(xué)命題［14］。由于坡溝系統(tǒng)重力侵蝕的發(fā)生具有點多面廣和隨機性的特點，且基礎(chǔ)觀測資料相對薄弱，有關(guān)研究就更少［2，13］。

淤地壩（系）作為黃土高原地區(qū)溝道水土保持的重要工程措施，在減蝕、減沙、防洪、保收、發(fā)展地方經(jīng)濟等方面起到了重要作用［15－17］，而且有機結(jié)合了水土保持和致富的關(guān)系，深受黃土高原地區(qū)群眾的喜愛，同時又為水利水土保持部門所關(guān)注。修壩攔泥淤地、蓄水?dāng)r沙是流域水土保持綜合治理的一項重要措施，目前國內(nèi)外有關(guān)淤地壩的研究主要集中在壩系相對穩(wěn)定性，減蝕減沙，滯洪保收和生態(tài)環(huán)境效應(yīng)等方面，而對淤地壩控制坡溝系統(tǒng)重力滑坡侵蝕方面的研究則少有人提及。壩地淤高后抬高了溝底的侵蝕基準(zhǔn)面，增強了溝坡的穩(wěn)定性，減少了溝坡重力侵蝕發(fā)生的可能性。那么坡溝系統(tǒng)的穩(wěn)定性和重力滑坡侵蝕中峁坡和溝坡對溝道侵蝕產(chǎn)沙的貢獻隨著壩地的逐漸增高到底是如何變化的呢？本研究通過數(shù)值模擬，對隨著壩地的逐漸淤高而產(chǎn)生的坡溝系統(tǒng)的穩(wěn)定性、滑塌量以及坡溝系統(tǒng)重力滑坡侵蝕中峁坡和溝坡所占比例等方面進行研究，以期為小流域溝道重力滑坡侵蝕的預(yù)報和分析評價提供重要參數(shù)和科學(xué)依據(jù)，并為評價坡溝系統(tǒng)的穩(wěn)定性提供一定可靠度的依據(jù)。

1 研究區(qū)概況

本研究選取王茂溝流域左岸的典型支溝——關(guān)地溝作為研究對象。王茂溝是陜北綏德縣韭園溝中游左岸的一條試驗性治理支溝，為無定河的一條2級支溝，海拔高度940～1 188m，流域面積5.97km2，主溝長3.75km，溝道平均比降為2.7%，溝谷地面積2.97km2，占流域總面積的46.7%，溝壑密度4.3km／km2［2，17］。地形地貌主要由梁、峁以及分割梁峁的溝谷組成，溝內(nèi)地形破碎，溝壑縱橫，坡陡溝深。梁頂和峁頂面積不大，均略成穹形，坡度在8°～10°，梁、峁坡的坡度，一般在20°～35°；梁、峁坡以下為溝谷，溝谷橫斷面在上游及支溝均呈“V”字型，在下游略呈“U”字型。溝谷坡極陡，一般在35°以上。按土壤侵蝕區(qū)劃，屬黃土高原丘陵溝壑區(qū)第Ⅰ副區(qū)［2，17］。

流域?qū)俦睖貛Ц珊荡箨懶詺夂?，年平均氣?0.2℃，夏季多東南風(fēng)，春秋多西北風(fēng)。流域多年平均降水量513.1mm，降雨量年際變化大，年最大降雨量是年最小降水量的3.5倍；降雨年內(nèi)分配極不均勻，年內(nèi)降雨量主要集中在汛期7—9這3個月，汛期降水量占年降水量的73.1%，且多以暴雨形式出現(xiàn)，一次暴雨產(chǎn)沙量往往為全年產(chǎn)沙量的60%以上，土壤侵蝕以水力侵蝕及重力侵蝕為主，治理前屬劇烈侵蝕區(qū)，在黃土丘陵溝壑區(qū)具有一定的代表性［2］。

2 研究方法

2.1 重力滑坡侵蝕分析的極限平衡法

目前應(yīng)用于邊坡穩(wěn)定分析的方法主要有基于極限平衡的傳統(tǒng)方法和有限元法［18］。極限平衡法（如瑞典圓弧法、Bishop法、Janbu法）是邊坡穩(wěn)定分析中最常用的方法，它通過分析坡體在臨近破壞的狀況下，土體外力與內(nèi)部強度所提供抗力之間的平衡，計算土體在自身和外荷作用下的土坡穩(wěn)定性的程度。傳統(tǒng)的邊坡穩(wěn)定性分析方法中，為了便于分析計算做了許多假設(shè)［19－22］，如假設(shè)一個滑動面、不考慮土體內(nèi)部的應(yīng)力—應(yīng)變關(guān)系等。因此，傳統(tǒng)分析方法不能得到滑體內(nèi)的應(yīng)力、變形分布狀況，也不能求得巖土體本身的變形對邊坡變形及穩(wěn)定性的影響。在邊坡穩(wěn)定分析中引入有限元法始于20世紀(jì)70年代［20］，有限元法克服了傳統(tǒng)分析法的不足，不僅滿足力的平衡條件，而且還考慮了土體應(yīng)力—變形關(guān)系，能夠得到邊坡在荷載作用下的應(yīng)力、變形分布，模擬出邊坡的實際滑移面［19］。正因為有限元法的上述優(yōu)點，近年來已廣泛應(yīng)用于邊坡穩(wěn)定性分析。

2.2 重力滑坡侵蝕的有限元強度折減法原理

有限元強度系數(shù)折減法的基本原理［19－20］是將坡體強度參數(shù)黏聚力c和內(nèi)摩擦角φ，同時除以一個折減系數(shù)ω，得到一組新的c′和φ′：

式中：c——巖土體原本的黏聚力；φ——巖土體原本的內(nèi)摩擦角（°）；ω——巖土體強度折減系數(shù)；c′——巖土體折減后的黏聚力；φ′——巖土體折減后的內(nèi)摩擦角（°）。

然后作為新的參數(shù)輸入，再進行試算，利用相應(yīng)的穩(wěn)定判斷準(zhǔn)則，程序自動根據(jù)彈塑性計算結(jié)果得到破壞滑動面，確定相應(yīng)的ω值為坡體的最小穩(wěn)定安全系數(shù)。有關(guān)研究表明［18－22］，有限元強度折減法的安全系數(shù)在本質(zhì)上與傳統(tǒng)方法是一致的。趙尚毅等［19－20］通過多種比較計算說明有限元折系數(shù)法用于分析土坡穩(wěn)定問題是可行的，但必須合理地選用屈服條件以及嚴(yán)格地控制有限元法的計算精度。

本研究采用節(jié)點不平衡力的不收斂作為判據(jù)，同時根據(jù)塑性區(qū)的范圍及其連通狀態(tài)確定相應(yīng)的安全系數(shù)，屈服準(zhǔn)則采用不等角六邊形外接圓（DP1）屈服準(zhǔn)則計算安全系數(shù)然后按式（2）轉(zhuǎn)換為莫爾—庫侖等面積圓準(zhǔn)則下的安全系數(shù)［2］：

式中：η——外接圓（DP1）屈服準(zhǔn)則下安全系數(shù)轉(zhuǎn)換為莫爾—庫侖等面積圓屈服準(zhǔn)則下安全系數(shù)的轉(zhuǎn)換系數(shù)；φ——巖土體的內(nèi)摩擦角。

3 結(jié)果與分析

3.1 關(guān)地溝小流域溝坡系統(tǒng)概化模型及有限元模型

根據(jù)地質(zhì)資料［17］，關(guān)地溝地層構(gòu)造主要是馬蘭黃土（Qeol3），梁、峁頂、峁坡均有分布，厚度20～30m，其下層為離石黃土（Qeol2），厚度50～100m，多出露于溝坡上，再下層主要是三疊紀(jì)砂頁巖層，基本接近水平，多出露于干溝、支溝的下游及其兩側(cè)［17］。關(guān)地溝坡溝系統(tǒng)概化模型土層從上到下分別為馬蘭黃土（Qeol3）和離石黃土（Qeol2），厚度分別取為土層厚度的平均值25和75m。通過實測數(shù)據(jù)的統(tǒng)計分析，峁坡坡度在21°～25°所占比例最大，達到50%，而實測坡度在21°～25°范圍內(nèi)的平均值為23.1°，取峁坡坡度為24°；溝坡坡度在36°～45°所占比例最大，達到43.6%，而實測坡度在36°～45°范圍內(nèi)的平均值為40.8°，取溝坡坡度為41°［2］。

坡溝系統(tǒng)的計算采用平面應(yīng)變條件下摩爾—庫侖不等角六邊形外接圓屈服準(zhǔn)則（DP1）和非關(guān)聯(lián)流動的理想彈塑性模型［23－25］。坡溝系統(tǒng)基巖底邊界為固定約束，左邊界為水平約束。彈塑性有限元計算采用大變形靜態(tài)分析選項，最大迭代次數(shù)為1 000［2］。

3.2 計算參數(shù)取值

根據(jù)邊坡工程經(jīng)驗、現(xiàn)場查勘及室內(nèi)巖土物理力學(xué)試驗［17，26－29］和有限元計算的需要，坡溝概化模型各土層材料特征取值如表1所示。

表1 關(guān)地溝坡溝系統(tǒng)的計算參數(shù)

3.3 關(guān)地溝小流域隨壩地淤高溝坡系統(tǒng)穩(wěn)定性的變化

關(guān)地溝4號壩從1959年修建到1987年被洪水沖毀的28a運行過程中，累計於高4.86m，年均淤高0.18m?？鍓魏螽?dāng)?shù)厝罕妼误w進行了零星修復(fù)填筑及2005年壩體修復(fù)加高，壩地的自然淤積過程遭到破壞，壩地已不再是原來垮壩前的自然淤高，截至2007年壩地累計高10.8m。由于關(guān)地溝4號淤地壩垮壩前壩地的年均於高僅為0.18m，在垮壩前的淤高模擬計算中以5a為一個時間跨度，計算壩地每淤高0.9m時的坡溝系統(tǒng)穩(wěn)定特征參數(shù)；垮壩后為了與垮壩前保持一致性，壩地以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 的 淤 積 過 程［2］，每 淤積 （增 高）0.9m其相應(yīng)的穩(wěn)定系數(shù)和滑塌量的散點圖見圖1—2。

圖1 隨壩地淤高坡溝系統(tǒng)穩(wěn)定系數(shù)的變化規(guī)律

圖2 隨壩地淤高坡溝系統(tǒng)滑塌量的變化規(guī)律

由圖1—2分別擬合出相關(guān)系數(shù)最高的隨壩地淤高坡溝系統(tǒng)的穩(wěn)定系數(shù)和滑塌量變化規(guī)律的相關(guān)方程［2］（如表2所示）。

表2 隨壩地淤高坡溝系統(tǒng)的穩(wěn)定系數(shù)和滑塌量變化規(guī)律

由表2可以得出，隨壩地淤高坡溝系統(tǒng)的穩(wěn)定系數(shù)和滑塌量分別滿足二項式和直線分布規(guī)律，且相關(guān)性較高，可用于該溝道重力侵蝕的分析評價［2］。

由表3可得，在原始溝道建壩初期壩地還沒有淤積的時候，坡溝系統(tǒng)穩(wěn)定系數(shù)最小，由有限元強度折減法得此時坡溝系統(tǒng)在摩爾—庫侖不等角六邊形外接圓屈服準(zhǔn)則（DP1）下的安全系數(shù)為1.507，將其按式（2）轉(zhuǎn)換為莫爾—庫侖等面積圓屈服準(zhǔn)則下的安全系數(shù)為1.224，其與滑坡分析采用Spencer法求得此時的安全系數(shù)1.254 7的相對誤差為2.42%，在鄭穎等［19－20］的研究指出莫爾—庫侖等面積圓屈服準(zhǔn)則的安全系數(shù)與Spencer法的平均誤差5%范圍內(nèi)，精度滿足要求。限于篇幅，考慮到系統(tǒng)誤差的一致性，其余各壩地淤高的安全系數(shù)沒有兩兩比較，但精度滿足要求。

3.4 發(fā)生重力滑坡侵蝕時峁坡和溝坡各自的滑塌量及比例

在淤地壩減蝕減沙和坡面水土保持措施效益分析計算中，一直存在一個難點就是如何如把坡面侵蝕和溝坡侵蝕的量分別計算出來。由于重力滑坡侵蝕的總體積V和總重量M是已知的（已通過軟件模擬計算出來，見表3），通過式（3）及推導(dǎo)式（4）—（5）可分別計算得到峁坡和溝坡在重力滑坡侵蝕中各自滑塌土體的體積和重量及相應(yīng)的比例（表4）。其中，Vm和Vg為峁坡和溝坡重力滑坡侵蝕土體的體積；γm和γg為峁坡和溝坡重力滑坡侵蝕土體的容重。

通過變量代換和代入消元法可得到：

表3 隨壩地淤高坡溝系統(tǒng)穩(wěn)定系數(shù)和滑塌量計算結(jié)果

由表3—4可以看出，在坡溝系統(tǒng)重力滑坡侵蝕中，隨著淤地壩的淤高，峁坡和溝坡滑塌物的體積和重量之比基本為1∶50；峁坡滑塌物的體積和重量占總滑塌物體積和重量的2%，而溝坡滑塌物的體積和重量占總滑塌物體積和重量的98%。在重力侵蝕產(chǎn)沙中，溝坡侵蝕是溝道泥沙的主要來源。

4 結(jié)論

基于邊坡分析軟件和有限元技術(shù)，通過建立關(guān)地溝4號壩上游典型坡溝系統(tǒng)的概化模型，得出隨淤地壩壩地淤高，坡溝系統(tǒng)穩(wěn)定系數(shù)的增加滿足二項式分布規(guī)律，而滑塌量則滿足線性減少的分布規(guī)律，兩者精度均較高，可用于該溝道重力侵蝕的預(yù)報、侵蝕量的分析評價［2］。

表4 坡溝系統(tǒng)滑坡侵蝕中峁坡和溝坡侵蝕量及所占比例

通過公式推導(dǎo)和軟件計算分析，分別得到了峁坡和溝坡的滑塌量及各自在滑坡侵蝕總量中所占的比重。在坡溝系統(tǒng)重力侵蝕中，隨著淤地壩的淤高，峁坡和溝坡滑塌物的體積和重量之比基本為1∶50；峁坡滑塌物的體積和重量占總滑塌物的體積和重量的2%，而溝坡滑塌物的體積和重量占總滑塌物的體積和重量的98%。在重力侵蝕產(chǎn)沙中，溝坡重力侵蝕是溝道泥沙的主要來源。

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

［2］ 鄒兵華，袁潔，李占斌，等.淤地壩減輕坡溝系統(tǒng)滑坡侵蝕的數(shù)值模擬［J］.水土保持通報，2013，33（1）：265－270.

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

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

［5］ 韓鵬，倪晉仁，王興奎.黃土坡面細溝發(fā)育過程中的重力侵蝕實驗研究［J］.水利學(xué)報，2003（1）：51－55.

［6］ 焦菊英，王萬忠，李靖，等.黃土高原丘陵溝壑區(qū)淤地壩的淤地攔沙效益分析［J］.農(nóng)業(yè)工程學(xué)報，2003，19（6）：302－306.

［7］ 王光謙，薛海，李鐵鍵.黃土高原溝坡重力侵蝕的理論模型［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報，2005，13（4）：335－344.

［8］ 于國強，李占斌，李鵬，等.黃土高原小流域重力侵蝕數(shù)值模擬［J］.農(nóng)業(yè)工程學(xué)報，2009，25（12）：74－79.

［9］ 于國強，李占斌，張霞，等.黃土高原坡溝系統(tǒng)重力侵蝕數(shù)值模擬研究［J］.土壤學(xué)報，2010，47（5）：809－816.

［10］ 朱同新.黃土地區(qū)重力侵蝕發(fā)生的內(nèi)部條件及地貌臨界值分析［M］.北京：氣象出版社，1987.

［11］ 陳浩.黃土丘陵溝壑區(qū)流域系統(tǒng)侵蝕與產(chǎn)沙關(guān)系［J］.地理學(xué)報，2000，55（3）：354－363.

［12］ 付凌.黃土高原典型流域淤地壩減沙減蝕作用研究［D］.南京：河海大學(xué)，2007.

［13］ 魏霞，李勛貴，李占斌，等.淤地壩對黃土高原坡溝系統(tǒng)重力侵蝕調(diào)控研究［J］西安建筑科技大學(xué)學(xué)報，2009，41（6）：856－861.

［14］ 雷阿林.坡溝系統(tǒng)土壤侵蝕鏈動力機制模擬試驗研究［D］.陜西 楊凌：中國科學(xué)院水利部水土保持研究所，1996.

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

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

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

［18］ 黃正榮，梁精華.有限元強度折減法在邊坡三維穩(wěn)定分析中的應(yīng)用［J］.工業(yè)建筑，2006，36（6）：59－64.

［19］ 趙尚毅，鄭穎人，時衛(wèi)民，等.用有限元強度折減法求邊坡穩(wěn)定安全系數(shù)［J］.巖土工程學(xué)報，2002，24（3）：343－346.

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

［21］ Zienkiewicz O C，Taylor，R L.The Finite Element Method［M］.New York：McGraw－Hill，1989.

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

［23］ 吳海真，顧沖時.有限元強度折減法在土石壩邊坡穩(wěn)定分析中的應(yīng)用［J］.水電能源科學(xué)，2006，24（4）：54－56.

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

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

［26］ 王家臣.邊坡工程隨機分析原理［M］.北京：煤炭工業(yè)出版社，1996.

［27］ 祝玉學(xué).邊坡可靠性分析［M］.北京：冶金工業(yè)出版社，1993.

［28］ 龔曉南.土塑性力學(xué)［M］.2版.杭州：浙江大學(xué)出版社，1999.

［29］ 戚筱俊.工程地質(zhì)及水文地質(zhì) ［M］.北京：中國水利水電出版社，1996.