基于DEM方法的滑坡堰塞壩幾何特征分析

趙高文 喬建平 姜元俊 孟華君 王萌 楊宗佶

摘要：基于離散元法（ DEM）對滑坡堰塞壩的幾何特征進行了數(shù)值模擬，考慮溝道斷面形態(tài)、滑坡速度以及溝床坡度三種因素的影響，共27種組合。結(jié)果表明：溝道斷面形態(tài)、滑坡速度對堰塞壩橫向、縱向幾何形態(tài)都有影響，而溝床坡度只對堰塞壩縱向幾何特征產(chǎn)生影響;滑坡速度增大，堰塞壩縱剖面頂寬和底寬的比值減小，堰塞壩的有效厚度減小;溝床坡度增大，堰塞壩下游壩坡內(nèi)角減小，上游壩坡內(nèi)角增大，兩者的比值減小，堰塞壩壩坡傾角的變化規(guī)律與內(nèi)角的變化規(guī)律相反。

關(guān)鍵詞：滑坡;堰塞壩;泥石流;漫頂破壞;離散元

中圖分類號：P694;P642.22

文獻標(biāo)志碼：A

doi：10.3969/j.issn.1000- 1379.2019.05 .003

滑坡堰塞壩屬于堰塞壩的一種，是滑坡造成溝道或河道堵塞而形成的天然壩，該類壩體穩(wěn)定性很差。Costa等[1]統(tǒng)計分析了73座滑坡堰塞壩（堵江型）破壞案例，發(fā)現(xiàn)其中27%在形成后一天內(nèi)破壞、50%在形成后10 d內(nèi)破壞，破壞的主要原因是漫頂侵蝕。柴賀軍等[2]在Costa等的研究基礎(chǔ)上，增加24個破壞案例后進一步指出.21%的滑坡堰塞壩在形成后一天內(nèi)破壞.48%在形成后10 d內(nèi)破壞，78%在形成后半年內(nèi)破壞，88%在形成后一年內(nèi)破壞?；卵呷麎稳舳氯拥佬纬纱笮脱呷赡軙蜎]上游農(nóng)田與民宅，若潰決形成洪水或泥石流，往往給下游帶來巨大損失。例如：1786年6月四川瀘定摩崗嶺滑坡堰塞湖潰決，潰決洪水不僅摧毀了下游一切基礎(chǔ)設(shè)施，還造成了近10萬人死亡[3]：1933年疊溪滑坡堰塞湖潰決，造成下游近2 500人死亡[4]。

目前，對滑坡堰塞壩幾何特征的研究較少且多以經(jīng)驗性分析為主：Cui等[5]、Xu等[6]、Yin等[7]根據(jù)堰塞壩高度、庫容及壩體材料特性，對汶川地震災(zāi)區(qū)堰塞壩的危險性進行評估，并將其分為高危、危險、一般、低危等類型：王珊珊等[8]和Zhou等[9]基于GIS技術(shù)提取河谷橫斷面特征參數(shù)，評價滑坡堵江的可能性;Casagli和Ermini采用地貌指數(shù)法，建立了以堆積指數(shù)I_b、蓄水指數(shù)I_i[10]和無量綱堆積指數(shù)DBI[11]為控制因素的指數(shù)理論，用來評價堰塞壩的穩(wěn)定性;Komp[12]、Nash等[13]、Dong等[14]、Dal Sasso等[15]均采用地貌指數(shù)法進行了相應(yīng)的研究。除經(jīng)驗性分析方法外，數(shù)值模擬是較為常見的研究滑坡堰塞壩幾何形態(tài)的方法：Hun-gr[16]根據(jù)流變理論，建立了滑坡堰塞壩的動力學(xué)數(shù)值模型;Lo等[17]基于離散元法（ DEM），反演分析了臺灣小林村滑坡堰塞壩可能的堆積形態(tài)：Zhao等[18]采用CFD-DEM耦合的方法模擬滑坡堵江形成堰塞壩的過程。但現(xiàn)有關(guān)于滑坡堰塞壩幾何特征的研究仍然存在以下兩個問題：

（1）現(xiàn)有經(jīng)驗分析方法主要是根據(jù)已形成堰塞壩的幾何特征，分析堰塞壩未來的發(fā)展趨勢，并不能對堰塞壩的幾何特征作出預(yù)測。有的經(jīng)驗分析準(zhǔn)確性不夠，例如，KOILip[12]采用地貌指數(shù)法分析新西蘭境內(nèi)滑坡堰塞壩穩(wěn)定性時，結(jié)果并不理想。

（2）數(shù)值模擬研究不成熟。Hungr的研究并不適合滑坡體含水量較低、未發(fā)生液化的情況;Lo等的研究采用了小林村所在區(qū)域的真三維地形，這增加了模擬結(jié)果的真實性和可靠性，但由于地形過于特殊化和具體化，這種研究往往只能用于具體案例的反演分析，研究結(jié)果并非普遍適用，而且真三維地形模擬需要準(zhǔn)確的地形數(shù)據(jù)（包括滑坡前和滑坡后）;Dong等[19]指出，準(zhǔn)確的數(shù)字地形數(shù)據(jù)在絕大多數(shù)情況下是無法短時間內(nèi)獲得的：Zhao等的研究主要關(guān)注流固耦合的作用，未考慮地形條件等外在因素的影響，且該研究設(shè)置的離散元顆粒生成方式與野外真實的滑坡堰塞壩形成過程有較大區(qū)別。

1 堰塞壩幾何特征與潰口的關(guān)系

幾何特征對滑坡堰塞壩的潰決破壞有重要影響，在地貌指數(shù)法中，滑坡堰塞壩的高度、長度、寬度和體積等變量常常作為堰塞壩穩(wěn)定性分析的關(guān)鍵因素[19]。除了用于堰塞壩的穩(wěn)定性分析外，滑坡堰塞壩的幾何特征還是研究堰塞壩的潰決機理和過程以及開展模型試驗的基礎(chǔ)。由于滑坡堰塞壩的幾何形態(tài)很不規(guī)則，因此不同滑坡堰塞壩漫頂時的潰口位置不同，如著名的疊溪和紅石巖滑坡堰塞壩，如圖1所示。

絕大多數(shù)滑坡堰塞壩的破壞方式為漫頂侵蝕，由圖1可見，堰塞壩的幾何特征決定了堰塞壩漫頂時的潰口位置，潰口位置又影響漫頂時壩體的侵蝕過程。當(dāng)潰口位于壩體側(cè)面時，壩體受到的側(cè)向侵蝕為單側(cè)侵蝕，當(dāng)潰口位于壩體中間時，側(cè)向侵蝕為雙側(cè)侵蝕。這種因堰塞壩幾何特征引起的侵蝕差異對于滑坡堰塞壩的潰決過程和輸沙過程有重要意義。筆者曾做過驗證性試驗，模擬側(cè)向侵蝕差異對潰壩過程的影響，兩個壩體除初始潰口外，其他條件均相同，結(jié)果如圖2所示，可以看出，雙側(cè)侵蝕的效果明顯強于單側(cè)侵蝕，進而導(dǎo)致潰口的擴展速度和輸沙濃度等參數(shù)也不相同。

綜上所述，滑坡堰塞壩的幾何特征是堰塞壩的穩(wěn)定性分析、破壞侵蝕過程分析以及開展模型試驗的基礎(chǔ)，研究滑坡堰塞壩的幾何特征具有重要意義。

滑坡堰塞壩的形成過程是一個復(fù)雜且有待深入了解的山坡一溝谷耦合過程[12.20].很多因素都可能對滑坡堰塞壩的幾何特征產(chǎn)生影響，包括滑坡類型、滑坡的運動以及堆積過程等[19]。在滑坡堰塞壩的所有誘發(fā)因素中，地震是最主要的誘發(fā)因素之-[l]，汶川地震后，整個災(zāi)區(qū)地層出現(xiàn)松動，形成了上千個不穩(wěn)定斜坡[21]，很多不穩(wěn)定斜坡失穩(wěn)后都有可能形成新的滑坡堰塞壩，這給次生地質(zhì)災(zāi)害的發(fā)生提供了條件。如果能預(yù)測這些滑坡堰塞壩的幾何特征，將為災(zāi)害識別和應(yīng)急搶險工作提供有力參考。

2 研究方法

采用離散元法（ DEM）對滑坡堰塞壩的幾何特征進行初步分析，顆粒接觸采用Hertz-Mindlin with JKR模型。

離散元法通過力一位移法則計算顆粒間的接觸力，將顆粒運動當(dāng)作顆粒振動進行求解，并根據(jù)牛頓第二定律分析顆粒運動（考慮重力的影響），其控制方程為

2.2 Hertz-Mindlin with JKR接觸模型

Hertz-Mindlin with JKR模型是將Hertz-Mindlin無滑動模型與JKR（Johnson-Kendall-Roberts）黏聚力模型[22]相結(jié)合的一種顆粒接觸模型，該模型的優(yōu)點是可以考慮黏聚力的影響，模型中顆粒切向接觸力、法向阻尼力和切向阻尼力的計算方法與Hertz-Mindlin無滑動接觸模型的計算方法相同[23]。但由于考慮了顆粒表面能（y）的影響，因此顆粒的接觸半徑（a）比Hertz-Mindlin無滑動模型中的接觸半徑（ao）要大。法向力（ F_JKH）和重疊量（δ，可以為負值，表示顆粒未發(fā)生物理接觸）的計算公式為

當(dāng)表面能γ=0時，Hertz-Mindlin with JKR模型即為Hertz-Mindlin無滑動接觸模型。

在巖土工程數(shù)值計算領(lǐng)域，Hertz-Mindlin with JKR模型可以考慮土體含水量和黏性顆粒的影響，因此具有較高的適用性。

3 數(shù)值模擬

3.1 物理模型

一般地，滑坡堰塞壩的動力過程包括：坡體失穩(wěn)一加速運動一減速堆積等過程，如圖4（a）所示。由于數(shù)值模擬考慮所有過程的計算量過大，因此有必要進行物理模型概化，本研究不考慮坡體啟動和加速運動過程，通過API（Application Programming Interface）編程直接賦予坡體相應(yīng)的速度，模擬滑坡體以某個速度滑人溝道并堆積成壩，滑坡體的形態(tài)也設(shè)為較規(guī)則的塊體，如圖4（b）所示。

3.2 影響因素

滑坡堰塞壩的幾何特征受多種因素影響，筆者在此建議將各類影響因素按能空間尺度和能量尺度進行分類，以便今后開展滑坡堰塞壩動力過程關(guān)鍵因素的提取。

從空間尺度和能量尺度而言，溝道斷面形態(tài)、滑坡速度（ν）以及溝床坡度（θ）是決定滑坡堰塞壩幾何特征的基本參數(shù)，本研究主要考慮這三種因素對滑坡堰塞壩幾何特征的影響。

3.3 模擬條件

由于滑坡的規(guī)模變化非常大，從幾萬到幾百萬甚至上千萬立方米，因此本研究采用單位體積分析，滑坡體總體積0.5 m³.考慮到今后開展室內(nèi)模型試驗（滑槽模擬）的需要，此處將滑槽與溝道統(tǒng)稱為幾何體，其材質(zhì)按鋼材建模。滑坡體位于一矩形滑槽內(nèi)，尺寸為1.0mxl.0 mx0.6 m（長×寬×高），溝道斷面形態(tài)包括矩形、梯形和V形3種（見圖5），為保證滑坡物質(zhì)能夠完全堵斷溝道，3種溝道的斷面面積與滑槽斷面面積相等（ 0.6 m²）。同時，為保證不同溝道內(nèi)的計算結(jié)果具有可比性，溝道高度保持不變（0.6 m），即滑坡體的重力勢能保持不變，滑槽與溝道的夾角保持為30°，溝道縱向長度為5.0 m，滑坡體滑人溝道的角度與滑槽的底部平行。

滑坡堰塞壩的動力過程實質(zhì)上是能量的釋放和耗散過程，當(dāng)勢能保持不變時，坡體動能是影響滑坡堰塞壩動力過程的根本性因素。本研究設(shè)置的滑坡體滑人溝道的速度為單位速度（1 m/s）的整數(shù)倍，分別為1、2、3 m/s，代表低速、中速、高速（相對于溝道尺寸而言）三種條件，動能比值為1：4：9。溝床坡度設(shè)置為00、50和100三種，其中：00針對底床坡度較小的溝道或河道（統(tǒng)稱溝道），50和100針對底床坡度較大的普通溝道。

3.4 參數(shù)標(biāo)定

巖土工程數(shù)值計算對輸人參數(shù)較敏感，為保證模擬結(jié)果的可靠性，必須對輸人參數(shù)進行標(biāo)定。根據(jù)2.2節(jié)的介紹.Hertz-Mindlin with JKR模型中，影響顆粒接觸力的主要因素有表面能、靜摩擦因數(shù)、滾動摩擦因數(shù)、恢復(fù)系數(shù)等。本研究以驗證性試驗（圖2）中的土體材料作為標(biāo)定對象，該土體材料取自都江堰虹口鄉(xiāng)銀洞子溝滑坡堰塞壩，其力學(xué)參數(shù)已通過大型直接剪切試驗測出，黏聚力（c）和內(nèi)摩擦角（φ）分別為6.2kPa和32°。由于土體黏性較低，因此參數(shù)標(biāo)定方法以數(shù)值計算顆粒的自然堆積角（p）等于土體內(nèi)摩擦角時的參數(shù)組合作為數(shù)值計算的輸人參數(shù)。參數(shù)標(biāo)定結(jié)果見表1.所有輸入?yún)?shù)見表2.

4.1 溝道斷面形態(tài)的影響

因27種組合的幾何參數(shù)多達243個，此處僅以高速條件（相對于溝道尺寸而言）下，平直溝道內(nèi)滑坡堰塞壩的幾何特征為例進行分析。該條件下不同溝道內(nèi)滑坡堰塞壩的幾何數(shù)據(jù)見表3。

由表3可見，高速條件下.V形溝道內(nèi)的堰塞壩縱剖面底寬最大，梯形次之，矩形最小;與此規(guī)律相反，堰塞壩的壩坡傾角和內(nèi)角以矩形最大，梯形次之.V形最小。上述規(guī)律說明：V形溝道內(nèi)滑坡堰塞壩向溝道上、下游延伸堆積的范圍最大。同時，由表3還可看出，高速條件下，矩形溝道內(nèi)堰塞壩的壩頂寬度L_t=0，這說明堰塞壩的縱剖面為三角形，但梯形和V形溝道內(nèi)堰塞壩縱剖面均為梯形。低速條件下，三種溝道內(nèi)的堰塞壩縱剖面如圖7所示，全為梯形。由此可見，溝道斷面形態(tài)對堰塞壩的幾何形態(tài)有直接影響。

4.2 滑坡速度的影響

此處僅以溝床坡度θ=0°時梯形溝道內(nèi)堰塞壩的幾何特征為例進行分析，不同滑坡速度條件下堰塞壩的幾何數(shù)據(jù)見表4。

低速條件形成的堰塞壩，h_max與h_min的比值達到4：1，說明堰塞壩的橫斷面存在很大高差，這類堰塞壩比較常見，如Young河滑坡堰塞壩（圖8（a））和都江堰虹口鄉(xiāng)銀洞子溝滑坡堰塞壩（圖8（b））均屬該類型，其中銀洞子溝滑坡堰塞壩剛形成時的高差接近160m。橫斷面存在巨大高差的堰塞壩若發(fā)生漫頂潰決，則其潰口側(cè)壁失穩(wěn)滑動的規(guī)模會越來越大。

隨著滑坡速度的增大，堰塞壩橫斷面的高差h_max-h_min逐漸減小。由不同速度條件下梯形溝道內(nèi)堰塞壩的橫斷面圖（如圖9所示）可知，隨著滑坡速度增大，滑坡體可以向?qū)Π杜栏叨纬神R鞍狀的斷面形態(tài)，這與圖1中的兩個真實案例相吻合。由此可見，滑坡速度對堰塞壩漫頂潰決時的潰口位置有重要影響，低速條件下形成的堰塞壩，初始潰口位于壩體側(cè)面，側(cè)向侵蝕屬單側(cè)侵蝕（如圖l（a）、圖8所示），但當(dāng)滑坡速度增大后，堰塞壩的初始潰口向中部移動，側(cè)向侵蝕變?yōu)殡p側(cè)侵蝕（如圖l（b）所示）。

隨著滑坡速度的增大，堰塞壩頂寬（L_t）逐漸減小，堰塞壩頂寬與底寬（L_b）的比值亦減小。若以（L_t+L_b）/2表示堰塞壩的有效厚度，則隨著滑坡速度的增大，堰塞壩的有效厚度變小。

4.3 溝床坡度的影響

通過對試驗結(jié)果的分析可知，溝床坡度主要對堰塞壩的縱向堆積產(chǎn)生影響。27種組合對應(yīng)的壩坡傾角和內(nèi)角見表5。

4.3.1 壩坡內(nèi)角（φ'_u、φ'_d）

根據(jù)表5的統(tǒng)計結(jié)果，可得出以下認識：

（1）溝床坡度一定時，堰塞壩壩坡內(nèi)角的大小關(guān)系是矩形溝道≥梯形溝道≥V形溝道。V形溝道內(nèi)的堰塞壩壩坡內(nèi)角大小不受滑坡速度的影響。

（2）整體而言，溝床坡度增大，堰塞壩下游壩坡內(nèi)角（φ'_d）減小，上游壩坡內(nèi)角（φ'_u）增大，φ'_u與φ'_d的比值減小，如圖10所示。

5 結(jié)論

基于離散元法（ DEM）對滑坡堰塞壩的幾何特征進行模擬，采用可以考慮顆粒黏性的Hertz-Mindlinwith JKR模型進行了27種參數(shù)組合的模擬分析，得到以下幾點主要認識。

（1）滑坡速度和溝道形態(tài)對滑坡堰塞壩的橫向和縱向幾何形態(tài)都有影響，而溝床坡度主要影響滑坡堰塞壩的縱向幾何特征。

（2）低速條件下形成的堰塞壩在橫向存在較大高差，堰塞壩若發(fā)生漫頂潰決，初始潰口位于堰塞壩側(cè)面，側(cè)向侵蝕為單側(cè)侵蝕。當(dāng)滑坡速度達到高速條件后，滑坡體可以向?qū)Π杜栏叨纬神R鞍狀的橫斷面，這類堰塞壩若出現(xiàn)漫頂潰決，初始潰口位于壩體中部，側(cè)向侵蝕為雙側(cè)侵蝕。

（3）滑坡速度增大，堰塞壩縱剖面頂寬（L_t）和底寬（L_b）的比值減小，堰塞壩的有效厚度減小。

（4）溝床坡度增大，堰塞壩下游壩坡內(nèi)角減小，上游壩坡內(nèi)角增大，兩者的比值減小。但是堰塞壩壩坡傾角的變化規(guī)律和內(nèi)角的變化規(guī)律相反。

參考文獻：

[1]COSTA J E，SCHUSTER R L.The Formation and Failure ofNatural Dams[J].Ceological Society of America Bulletin，1988， 100（7）： 1054-1068.

[2] 柴賀軍，劉浩吾，劉漢超，等，天然土石壩穩(wěn)定性初步研究[J].地質(zhì)科技情報，2001， 20（1）：77-81.

[3]DAI F C，LEE C F，DENC J H，et al.The 1786 Earthquake-Triggered Landslide Dam and Subsequent Dam-Break Flood onthe Dadu River， Southwestem China [J]. Geomorphology，2005， 65（3）：205-221.

[4] 王蘭生，楊立錚，李天斌，等，四川岷江疊溪較場地震滑坡及環(huán)境保護[J].地質(zhì)災(zāi)害與環(huán)境保護，2000， 11（3）：195-199.

[5]CUI P，ZHUY Y，HAN Y S，et a1.The 12 May WenchuanEarthquake-Induced Landslide Lakes： Distribution and Prelim-inary Risk Evaluation[J]. Landslides， 2009，6（3）：209-223.

[6] XU Q，F(xiàn)AN X M， HUANC R Q，et al.Landslide DamsTriggered by the Wenchuan Earthquake， Sichuan Province，South West China[J].Bulletin of Engineering Ceology andthe Environment， 2009， 68（3）： 373-386.

[7]YIN Y P，WANC F W， SUN P.Landslide Hazards Triggeredby the 2008 Wenchuan Earthquake， Sichuan， China[J]. Land-slides， 2009， 6（2）： 139-152.

[8] 王珊珊，童立強，基于河谷橫剖面形態(tài)特征的滑坡體堵江易發(fā)性評價研究[J].地理與地理信息科學(xué)，2016， 32（5）：97-102.

[9]ZHOU L，TAKASHI 0.Longitudinal and Transverse Profilesof Hilly and Mountainous Watersheds in Japan[J].Geomor-phology， 2009， 111（1-2）： 17-26.

[10]CASACLI N，ERMINI L.Ceomorphic Analysis of LandslideDams in the Northem Apennine[J].Japanese Ceomorpho-logical Union， 1999， 20（3）：219-249.

[11]ERMINL，CASACLIN.Prediction of the Behavior of Land-slide Dams Using A Ceomorphologcal Dimensionless Index[J]. Earth Surface Processes and Landforms， 2003， 28（1）：31-47.

[12]KORUP 0.Geomorphometric Characteristics of New ZealandLandslide Dams[J].Engineering Geology， 2004， 73（1）：13- 35.

[13]NASH T， BELL D， DAVIES T， et al. Analysis of the For-mation and Failure of Ram Creek Landslide Dam， South Is-land， New Zealand [ J]. New Zealand Joumal of Ceologyand Ceophysics ， 2008， 51（ 3） ： 187-193.

[14]Dong J J， Tung Y H， Chen C C， et al. Logistic RegressionModel for Predicting the Failure Probability of A LandslideDam[J]. Engineering Ceology， 2011， 117（ 1） ： 52-61.

[15]DAL SASSO S F， SOLE A， PASCALE S， et al. AssessmentMethodology for the Prediction of Landslide Dam Hazard [J] .Natural Hazards and Earth System Sciences， 2014， 14（3） ：557-567.

[16]HUNCR 0. A Model for the Runout Analysis of Rapid FlowSlides， Debris Flows， and Avalanches [ J]. CanadianCeotechnical Journal， 1995 ， 32（4） ： 610-623.

[17]10 C M， LIN M L， TANC C L， et al. A Kinematic Modelof the Hsiaolin Landslide Calibrated to the Morphology ofthe Landslide Deposit [ J]. Engineering Geology， 2011，123（1） ： 22-39.

[18]ZHAO T， DAI F， XU N. Coupled DEM-CFD Investigationon the Formation of Landslide Dams in Narrow Rivers [J] .Landslides， 2017， 14（1） ： 189-201.

[19]DONC J J， LAI P J， CHANC C P， et al. Deriving Land-slide Dam Geometry From Remote Sensing Images for theRapid Assessment of Critical Parameters Related to Dam-Breach Hazards [ J] . Landslides， 2014， 11（ 1） ： 93-105.

[20]COSTA J E， SCHUSIER R L. Documented Historical LandslideDams from Around the World[ R] . Virginia ： US Ceologcal Sur-vey， 1991：91-239.

[21]殷躍 .汶川八級震地質(zhì)災(zāi)害研究 [J] .工程地質(zhì)學(xué)報，2008， 16（4） ： 433-444.

[22] JOHNSON K L， KENDALL K， ROBERTS A D. Surface Energyand the Contact of Ilastic Solids [ C] //Proceedings of the RoyalSociety of London A： Mathematical， Physical and Engneering Sciences. London ： The Royal Society， 1971： 301-313.

[23] SANTOS K G， CAMPOS A V P， OLIVEIRA O S， et al.Dem Simulations of Dynamic An～e of Repose of AcerolaResidue： A Parametric Study Using a Response SurfaceTechnique [J] . Blucher Chemical Engineering Proceedings ，2015 ， 1（ 2） ： 11326-11333.

[責(zé)任編輯張帥]