深圳“12.20”渣土場遠程流化滑坡動力過程分析

高 楊，衛童瑤，李 濱，賀 凱，劉 錚，王學良

(1.中國地質科學院地質力學研究所,北京 100081；2.中國科學院地質與地球物理研究所/中國科學院頁巖氣與地質工程重點實驗室,北京 100029；3.長安大學地質工程與測繪學院，陜西 西安 710054)

據統計，自20世紀90年代以來，我國的城市固體廢棄物排量占世界總排量的30%[1]，填埋場堆填成為最為有效的處理方式。然而，這些人工堆填的高陡邊坡，體積方量，一旦失穩就會造成滑坡災害，對周圍居民的生命財產安全形成巨大威脅。

目前在全球范圍已經發生了多起人工堆填體滑坡事件，如：20世紀70年代波黑薩拉熱窩堆填體滑坡[2]，1988年美國Kettleman Hills堆渣場滑坡[3-4]，1993年土耳其的伊斯坦布爾Umraniye Hekimbasi填埋場滑坡事件[5]，1996年美國Rumpke填埋場滑坡[6]，1997年南非德班Bulbul滑坡[7]，2000年菲律賓Payatas堆填體滑坡[8-9]，2005年印度尼西亞的Leuwigajah滑坡[7,10]，2015年我國的深圳紅坳渣土場滑坡[11-13]，等等。一般情況下，人工堆填邊坡失穩后呈現了滑體在坡腳處堆積的狀態，運動等效摩擦系數高于0.3，滑坡沖擊力相對較小，如：Payatas滑坡(0.45)、Bulbul滑坡(0.32)、Rumpke滑坡(0.33)等。同時也有另一部分堆填坡體失穩后形成遠程運動，并具有明顯的流化特征，運動等效摩擦系數僅為0.1左右，如：Leuwigajah滑坡(0.11)、深圳滑坡(0.10)、山西襄汾“9.8”滑坡(0.09)，該類型滑坡運動距離遠、沖擊力大、破壞力更強。人工堆填坡體通常距離人口居住區更為接近，坡體一旦失穩后形成遠程運動，往往會形成損失慘重的災難性事件。

在探尋人工堆填體滑坡后的遠程流化破壞運動堆積規律、成災風險的研究，數值模擬成為被廣泛應用的技術方法。目前遠程滑坡運動堆積數值模擬研究方法較多，主要分為單流體連續介質算法和顆粒流離散元算法。單流體連續算法，主要是基于滑坡運動等效流體理論的CFD、LBM、SPH和MPM等計算方法[14-20]，具有可視化、計算效率高等優勢，并且重點考慮到了滑坡運動過程中的動力本構模型，適用于流化型滑坡中的泥石流、泥流和山洪等災害的數值反演分析；顆粒流離散元算法，是基于非連續力學方法的PFC、EDEM和MatDEM等計算方法[21-23]，該方法基于不同的接觸類型和摩擦系數，能夠較好模擬大變形大位移，適用于崩塌、碎屑流的數值反演分析。另外，還有一些研究采用DEM-CFD的多相流耦合數值算法，考慮固體顆粒和流體之間的相互作用力。

在地質災害防災減災領域，針對不同滑坡類型能夠快速、高效地對影響區域進行預測和評估，才能真正發揮數值模擬技術的實用價值。本文以2015年12月20日發生的深圳紅坳渣土場遠程滑坡為研究實例，采用基于SPH方法的DAN3D的可視化、高效計算的數值模擬方法，建立三維滑坡模型，通過運動動力模型和參數選取，對高含水量的人工堆填體滑坡-泥漿流的全運動過程進行模擬分析，并對其運動堆積厚度、堆積形態及最大速度分布等進行深入討論，并總結了采用該分析方法下的人工堆填體遠程流化滑坡的動力模型和參數的適用范圍，結果適用于堆填場科學選址和危險區的快速精準評估。

1 深圳滑坡基本情況及失穩機制

2015年12月20日，廣東省深圳市光明新區紅坳渣土場發生特大型滑坡災害，滑體體積約2.75×106m3，滑坡前后緣長達1 100 m，造成77人遇難、33間房屋被毀。該滑坡主要物質組成為城市固體廢棄物，滑坡覆蓋面積為0.38 km2，滑程長達1 100 m，垂直高差達到113 m，等效摩擦角度僅為6°，屬于遠程運動滑坡的范疇[24]。

研究區地貌演化過程主要分為三個階段：采石場階段、渣土場階段和滑坡階段。

(1)采石場階段：渣土場前身為采石場，以人類工程活動開采花崗巖為主，通過遙感圖可以看出采石場后期已經形成“凹槽狀”深坑，并存有大量積水，為后期滑坡的發生埋下了巨大隱患(圖1a)。

(2)渣土場階段：采石坑的地形使之成為建筑垃圾渣土堆放的有利場地，并形成了體積為5.83×106m3、坡率約為1∶2.5的人工堆填邊坡。深圳滑坡的巖土結構主要由底部花崗巖基巖和上覆建筑垃圾渣土兩部分組成(圖1b)。

(3)滑坡階段：深圳滑坡發生后，一些研究者在災后應急調查工作基礎上，對該滑坡的失穩機制進行了探討，認為滑坡的觸發機理與地下水密切相關(圖1c)。殷躍平等[25]認為地表水入滲和固結滲流這兩種效應的疊加導致坡體穩定性降低，最終導致整體滑動；劉傳正[26]提出深圳滑坡為泥化地基、承壓浮托、堆載推擠和臨空滑移等綜合作用下的“竹筏效應”導致滑坡失穩破壞遠程運動；Peng等[27]和Zhan等[28]認為松散渣土體材料性質和地下水水位的迅速抬升是導致滑坡失穩破壞的關鍵因素。滑坡下滑后大量房屋被滑坡沖毀掩埋，造成多人傷亡(圖1d)。

圖1 深圳人工堆填骨坡發生前后的三個階段及其破壞Fig.1 Three stages before and after the landslide and its destroy of Shenzhen landfill

2 深圳滑坡運動堆積情況

深圳滑坡失穩啟動后，約2.75×106m3的渣土體高位滑出。滑坡體沿N20°W方向以流態化形式在受納場下方開闊地帶散開，形成了近似于喇叭形狀的堆積區。滑坡區域總面積約0.38 km2，南北主軸長約1 100 m，東西向最大寬度為630 m，最小寬度為150 m。

根據堆積特征，紅坳受納場滑坡影響區主要分為滑源區、流通區和堆積區(圖2)。(1)滑源區：滑源區位于紅坳受納場區域，后緣及剪出口位置高程分別為155 m和65 m。根據滑坡縱剖面Ⅰ-Ⅰ′，滑源區沿主滑方向水平長度約為460 m，寬度150～400 m，面積約為0.15 km2。滑前受納場渣土坡體的最大厚度為 95 m，平均厚度為66 m；滑動體的最大厚度達50 m，平均厚度約為35 m；滑后滑源區的平均堆積厚度為15 m；滑帶角度為4°。(2)流通區：該區域位于紅坳收納場的收口處，起始位置是滑坡前緣剪出口，滑前高程為65 m。流通區東西兩側地形較高；中間較低，是滑坡向前流通運動的主要通道。該區長約100 m，寬度為150～260 m，面積為0.02 km2。巨大體積的滑體從該區域失穩剪出后，以流態化的運動狀態經過流通區。(3)堆積區：堆積區主要位于受納場下方的工業園區，沿主剖面方向長約540 m，最大寬度為620 m，影響面積為0.21 km2。堆積區的最大堆積厚度為22 m，平均堆積厚度為10 m。滑坡發生前堆積區位置處為工業生產園區，滑坡后緣距離臨時建筑物位置約有720 m，距離最近廠房約820 m。

圖2 滑坡后破壞Ⅰ-Ⅰ’堆積平面圖及剖面圖Fig.2 Plane graph and I-I’ accumulation profile of Shenzhen landfill landslide

3 計算方法及基底阻力模型

SPH算法是典型的拉格朗日方法，它的基本原理就是通過粒子模擬流體的運動規律，適合于求解高速碰撞等動態大變形問題。DAN3D滑坡模擬軟件采用SPH方法，將滑坡體總體積分為多個粒子，每個粒子都有一個有限的體積，材料的密度被認為是一致的。流體單元柱體的厚度和該區域的體積呈正比，由相鄰顆粒深度相加而得，所以流動體單元柱體的深度可用插值總和計算得到，即式(1)和(2)。SPH簡單示意圖見圖3。

圖3 SPH自由插值法概化示意圖(據McDougall[32])Fig.3 Schematic diagram of SPH free interpolation(from McDougall[32])

流體厚度梯度公式：

(1)

式中：V——每個顆粒的體積；

W——差之內核；

i和j——顆粒編號。

模型采用高斯插值內核：

(2)

式中：e——粒子平滑長度(所有粒子的e相等)，這是一個衡量內核寬度的值，即確定每個粒子的影響半徑[29]。

(3)

式中：N——總粒子數量;

B——光滑系數，無量綱量。

總之，SPH方法對于給定的每個顆粒體積和給定時間條件下，厚度和厚度梯度可以通過參考柱體區域計算出，并滿足連續性，根據動態方程在下個時步繼續計算。可用于遠程滑坡-碎屑流運動過程中堆積特征的動態分析。

該模擬方法基于結合流體深度的圣維南方程拉格朗日解，適用于表面層流、湍流、一般土體流動及碎屑流等。根據右手準則，在笛卡爾坐標系中(x，y，z)進行單元分析，z為滑動基床法線方向。

(1)質量守恒連續方程：

(4)

(2)動量守恒單元方程：

(5)

式中：ρ——滑體密度；

T——剪應力張量；

t——時間；

g——重力加速度；

▽——拉普拉斯算符；

v——速度；

?——張量積；

·——點積。

(6)

(7)

(8)

式中：ρ——等效流體密度；

h——流體深度；

vx、vy——流體x、y方向的速度；

τzx、τzy——基底剪切力；

kx、ky、kxy、kyx——切應力系數;

σz——等效流體對基底正應力。

式(6)和(7)中等號右邊第一項為滑體重力提供，第二和第三項為運動材料側向壓力，第四項為基底剪切力。式(8)在平行滑坡運動路徑x方向上考慮運動剪切摩阻力和側向壓力；在垂直滑坡運動路徑y方向上主要考慮土體側向壓力；在滑坡運動豎直方向上由于剪切力τxz和τyz同正應力σz相比較小，可忽略不計，因此不考慮。

4 模型建立及參數選取

4.1 滑坡模型建立

研究區地形網格文件主要為滑坡運動路徑文件和滑體網格，用來模擬滑體沿運動路徑的下滑運動，同時可以將運動路徑分為多個運動區域，分階段選用不同的基底阻力模型進行分析。在滑坡成災模式研究中，根據1∶2 000比例尺的數字高程模型，建立地形演化三個階段(采石場階段、填埋場階段和滑坡階段)，建立三維路徑地形和滑體數據網格文件(圖4)，隨后將兩個網格數據輸入到軟件中進行模擬分析。在DAN3D中選取組合基底阻力模型和模型參數在模擬過程中，滑體初始體積約為2.75 ×106m3，滑源區面積約為0.15 km2。

圖4 研究區DEM三維數字高程模型Fig.4 DEM three-dimensional digital elevation model of the research area

4.2 基底阻力模型選取

在分析中可以選用不同流體模型模擬不同類型和不同運動狀態下滑體的基底摩擦阻力對運動堆積過程的響應關系,在DAN3D中可采用等效流體運動阻力模型。目前針對滑坡運動模擬，在國際上最為常用的流體模型主要有三種：摩擦模型(Frictional model)、Voellmy模型和賓漢姆模型(Bingham model)。

(1)摩擦模型(Frictional model)，通常被選用于滑坡運動的初始啟動階段，且考慮了孔隙水壓力的重要影響作用。模型公式為：

τ=σ(1-ru)tanφ

(9)

式中：τ——滑體運動基底剪切阻力；

ru——孔隙水壓力系數；

φ——動力摩擦角；

σ——滑動路徑上正壓力。

(2)Voellmy模型(Voellmy model)，通常被選用于滑坡運動的流通和堆積階段，且考慮了運動過程中的湍流項和速度效應，剪切阻力和運動速度呈正比。模型公式為：

(10)

式中：σf——滑體運動摩擦基底剪切阻力，參數和摩擦模型一致；

γ——滑體重度；

v——滑體運動速度；

ξ——湍流系數。

(3)賓漢姆模型(Bingham model)，通常被選用于滑坡運動的流通和堆積的尾部階段，且考慮了流態化滑體的屈服應力和粘滯系數的影響作用。模型公式為：

(11)

式中：τ——滑體運動基底剪切阻力；

τy——滑體屈服應力；

v——滑體運動速度；

μB——賓漢姆粘滯系數；

h——滑體厚度。

摩擦模型通常用于開闊地形，地形起伏較平緩，速度影響較小的運動階段，更適用于滑坡滑源區失穩運動模擬，如：土質、巖質滑坡等；Voellmy流變模型更適合溝道中的湍流流體運動，模型針對溝渠地形，基底阻力與速度大小呈正比關系，更適用于滑坡鏟刮區域剪切，如：碎屑流、雪崩等；賓漢姆模型更適用于非牛頓泥漿流體的運動中剪切，如：泥石流、泥流等[18，30]。

4.3 模擬參數選取

在DAN3D的數值模擬參數設置中，主要包括SPH計算方法的控制參數和滑坡運動材料的模型參數。

(1)控制參數

SPH模擬方法控制參數主要包括顆粒數量(N)、粒子光滑系數(B)、速度平滑系數(C)、剛度系數(D)。連續體模擬是通過對控制方程進行離散化進行的，在滑動質量的每個重要位置，需要足夠多的計算單元(粒子)來捕捉其運動行為，因此計算分析中增加粒子數(N)，以提高SPH方法的計算精度；粒子平滑長度(B)影響粒子插值流深度的平滑度，并且可以調整，直到初始深度插值看起來平滑為止；速度平滑系數(C)決定相鄰粒子的速度對中心粒子的影響程度，速度平滑系數的選取控制了粒子擴散度，消除了運動過程中的大變形，增加運動過程中的穩定性；無量綱剛度系數(D)控制主動和被動型的滑體內應力狀態之間的應變率。針對深圳光明新區填埋場滑坡運動過程的分析，計算控制參數選取為：N=2 000，B=6，C=0.2，D=600。

(2)模型參數

為使滑坡后破壞的運動堆積模擬得到更精確的結果，通常采用試驗測試方法進行參數選取，包括模型中滑體密度、內摩擦角、黏聚力和動力摩擦角度等；其他模型參數可基于半經驗的試錯法進行確定，一些研究者已經在出版的著作中提出相關取值范圍，亦可作為參考值[31-33,18]。通常情況下，摩擦模型中孔隙水壓力系數(ru)參考范圍為0.1～0.9，基底動力摩擦角(φ)為10°～30°；Voellmy模型中孔隙水壓力系數(ru)參考范圍為0.5～0.8，動力摩擦系數(f)為0.05～0.1，湍流系數(ξ)為200～800 m/s2；賓漢姆模型中滑體屈服應力(τy)參考范圍為0.01～400 kPa，賓漢姆粘滯系數(μB)參考范圍為0.001～50 kPa·s，不同滑坡類型和滑體材料取值不同，具體參數參照實際的滑體屬性確定。

在深圳“12.20”滑坡實例模擬分析中，我們分別計算了三種組合模型的運動堆積過程，以此得到該類型滑坡最為適用的模型及參數。滑源區為滑體土失穩下滑，液化作用明顯，因此在滑源區均選用考慮高孔隙水壓力的庫倫摩擦模型；將流通區和堆積區選為同種模型，先后選用摩擦模型、Voellmy流體模型和賓漢姆流體模型，模型參數取值見表1～3，分界位置在距離滑坡后緣580 m處。三種情況下滑體容重和內摩擦角選用相同試驗參數：容重γ=16 kN/m3、內摩擦角φi=17°，再根據具體不同模型進行參數選取。

表2 摩擦-Voellmy組合(FV)模型參數Table 2 Parameters of Frictional-Voellmy combined model

表3 摩擦-Bingham組合(FB)模型參數Table 3 Parameters of Frictional-Bingham combined model

5 數值模擬結果

5.1 堆積狀態對比結果

圖5～6分別顯示了三種組合模型下的最終堆積等值線圖和堆積厚度分布剖面圖，通過對比分析可以看出FV組合模型運動距離最近、影響范圍和堆積區堆積厚度最小；FF組合模型運動距離、影響范圍居于三組模擬的中間，同實際滑坡堆積情況較為一致：影響范圍同實際情況相對較小，堆積區堆積厚度比實際情況小。FB組合模型運動距離、影響范圍居于三組模擬的最大值，堆積區厚度同實際滑坡堆積情況較為一致，但模擬結果影響范圍超出了實際堆積范圍，體現了滑體流動性的運動特征；因為考慮了運動路徑上房屋的阻擋作用，該模擬結果最適合于深圳滑坡堆積結果的反演分析。通過對比分析，摩擦模型和Voellmy模型雖有超孔隙水壓力的減阻作用，摩擦阻力仍然較大；摩擦模型在不考慮速度的情況下，剪切摩阻力小于帶有速度項的Voellmy模型。由于賓漢姆流體結模型剪切摩阻力小于以上兩者，更適合飽和渣土體泥漿的流態化運動情況。FB組合模型體現了滑體從滑源區剪出潰散解體后，呈泥流狀狀態運動堆積直到停止，流體的觸變和剪切變稀特性使滑體剪切屈服降低，以至于形成遠程滑坡。

三組模擬試驗都反映出，滑源區中具有高孔隙水壓力的摩擦模型能較好還原高含水量渣土體的失穩下滑過程；流通區和堆積區采用賓漢姆流體阻力模型能更好地模擬滑坡流化運動堆積狀態，并且模擬結果同實際滑坡堆積情況最為一致(圖6)。FB組合模型的選取體現了水與滑體的相互作用：一方面地下水導致滑體內孔隙水壓力的升高，有效應力下降，剪切摩阻力降低，使得滑體加速下滑，體現在摩擦模型中的孔壓系數(ru)；另一方面飽水渣土剪出后呈現流態化運動，與固體剪切采用不同的阻力模型，更能表現出滑體的流動狀態，剪切過程中摩擦阻力和運動速度呈正比關系，但流體屈服應力更小。因此，最終得出FB組合模型最適合深圳人工堆填體滑坡從失穩下滑到后破壞流化運動堆積的數值反演。

圖5 三組模型最終堆積厚度等值線圖Fig.5 Final stack thickness contour map of three groups of models

圖6 三組模型主滑方向堆積厚度變化剖面對比圖Fig.6 Comparison of the thickness variation profiles of the main sliding direction of the three groups of models

5.2 最大運動速度對比結果

圖7顯示了三種不同模型組合下深圳光明新區滑坡運動過程最大速度分布圖。FV組合模型最大速度為23 m/s，位于距離滑坡后緣590 m處，速度變化呈現上升快下降快的趨勢，體現了模型中速度同摩擦阻力成正相關關系的固體湍流剪切；FF組合模型最大速度為28 m/s，位于距離滑坡后緣620 m處；24～28 m/s的速度分布區域大于FV組合模型，但速度下降較快，滑體達到最大速度后，隨即很快下降為0 m/s，體現了模型中速度與剪切阻力不相關的高孔壓和低摩擦的固體滑體剪切。FB組合模型最大速度為30 m/s，位于距離滑坡后緣620 m處；24～30 m/s速度分布區域小于FF組合模型，但速度下降緩慢，滑體持速效應明顯，滑體堆積形態持續調整后速度降為0 m/s，體現了滑體呈現流態化運動。

圖7 三組模型最大速度分布等值線圖Fig.7 Maximum velocity distribution contour map of the three sets of models

組合模型計算結果(表4)顯示，FB組合模型更適用于深圳光明新區滑坡的后破壞運動模擬分析。深圳光明新區滑坡的速度經歷了“啟動-加速-持速-減速”的運動過程，水是影響滑動速度變化的關鍵因素。滑體破壞失穩之后，由于地形因素滑體經過勢動能轉換，滑動速度不斷增加。當滑體到達流通區時，滑動速度持續增加，最大速度均出現在流通區位置。隨后滑體呈流態化繼續向堆積區運移，在運動路徑上掩埋和沖毀33棟房屋建筑。

表4 組合模型計算結果Table 4 Combined model calculation results

5.3 FB模型運動過程結果

圖8顯示了DAN3D軟件中FB組合模型對深圳“12.20”填埋場滑坡下滑運動堆積過程的模擬結果，體現了滑體在滑源區固體剪切(F模型)、流通和堆積區流化剪切(B模型)的不同時段的運動堆積等值線圖，并采用紅黃藍顏色進行填充以表現滑體堆積厚度變化過程。在圖中滑體主滑方向為N 20° E，堆積距離為1 200 m，運動時間為70 s，紅色輪廓線表示了滑坡真實堆積范圍。FB組合模型下滑坡運動堆積動態過程如下：

圖8(a)顯示填埋場坡體失穩下滑的初始狀態，反映了滑坡失穩前渣土體堆填形態；

圖8(b)顯示填埋場滑坡的啟動狀態，反映了滑體開始高位剪出、啟動向下游運動；

圖8(c)顯示滑體移動下滑，顯示了在地形條件的影響下滑體經過流通區繼續向前運動，滑體呈流態化的運動特征向下游運動，且在流通區堆積厚度逐漸增加；

圖8(d)顯示滑體移動下滑到達堆積區，顯示了滑體主滑方向為N20°W，滑體前緣到達堆積區后呈流態化向周圍擴散運動，速度持續增加；

圖8(e)顯示滑體到達堆積區后再無側限阻擋的條件下逐漸散開，進入堆積區后滑體繼續向前運動，并向周圍擴散，影響范圍逐漸增大，速度逐漸下降；

圖8 FB組合模型流體基底剪切阻力模型運動堆積等值線圖Fig.8 Accumulation contour map of FB combination model fluid base shear resistance model motion

圖8(f)顯示滑體運動堆積的結束階段，滑體運動的最大范圍超過了真實滑坡堆積范圍。滑源區的最大堆積厚度位于滑源區的西北側，約為53 m；流通區最大厚度位于流通區的中部，約為18 m；堆積區最大厚度位于堆積區的中部，約為15 m。堆積范圍和堆積距離同實際情況較為一致。

6 結論

(1)滑坡在后破壞運動過程中體現固-流剪切轉換特征，主要可分為兩個階段：在滑源區內運動階段，為高孔隙水壓力下滑動剪切；在流通區和堆積區內運動階段，水動力作用使渣土體呈現較高的流態化，體現了高飽和度滑體的流動剪切。分階段的運動剪切方式增加了滑坡體的運動影響范圍。

(2)摩擦模型適合模擬堆渣土的孔隙水壓力作用下的失穩下滑剪切過程，適用于高含水量的人工堆填體滑坡的滑源區；賓漢姆模型適合模擬非牛頓流體飽和渣土體的流化剪切過程，適用于流通區和堆積區的流態化運動反演模擬分析。

(3)深圳滑坡后破壞運動速度變化主要經歷了“啟動-加速-持速-減速”的運動過程，固-流轉化導致滑坡造成巨大傷亡損失，其中水是影響滑動速度變化的關鍵因素。

(4)模擬結果顯示：堆積區平均堆積厚度為11 m，堆積范圍為0.4 km2，最大運動速度為30 m/s，最大速度發生于距滑坡后緣620 m處，堆積范圍、堆積厚度和運動速度同滑坡實際值基本一致。