基于PFC3D位移判據的低滲堆積層滑坡臨界降雨閾值反演

左雙英，蒲 泉，史文兵，王 勇，吳道勇，王安禮

(1.貴州大學 資源與環(huán)境工程學院，貴州 貴陽 550025； 2.貴州大學 喀斯特地質資源與環(huán)境教育部重點實驗室，貴州 貴陽 550025；3.貴州省有色金屬和核工業(yè)地質勘查局物化探總隊，貴州 都勻 558099； 4.貴州省質安交通工程監(jiān)控檢測中心有限責任公司，貴州 貴陽 550081)

貴州省地處云貴高原的斜坡地帶，受河流的侵蝕切割，地勢起伏大、地表相當破碎，地質災害頻繁。隨著城鎮(zhèn)化建設的步伐加快，對山地斜坡的開發(fā)利用加快，加之貴州夏季多雨的氣候條件，各類崩滑事故愈加頻繁。近年來，關嶺崗烏[1]、福泉[2]、納雍張家灣[3]、水城雞場[4]等幾起重大地質災害事故對當地民生危害巨大。上述崩滑事故中，降雨誘發(fā)的堆積層滑坡是其中最主要的滑坡類型之一，其特點是堆積體表層土粒徑粗、孔隙大、吸水性強，下部為細顆粒粉土或黏土，滲透性低，在短時強降雨作用下迅速吸水入滲，在坡體內部大量聚集來不及排出，導致坡體孔隙水壓力增大，有效應力、抗剪強度降低[5]，破壞了原有的平衡狀態(tài)而產生滑坡。因此，降雨型堆積層滑坡坡體物質組成與降雨入滲地下水復雜作用的特殊性給斜坡傳統(tǒng)力學評價預測方法的有效應用帶來了極大的困難與不確定性[6-7]。為此，深入研究滑坡的先兆因素，得出易于監(jiān)測的指標，提高滑坡預警預報的可操作性和準確性是業(yè)內亟待解決的重點和關鍵問題。

滑坡的預警研究自日本學者齋藤率先提出的經驗公式至今已有近50余年的發(fā)展歷程[8]。期間，國內外眾多學者針對這一問題展開系列研究，對滑坡的預警判據研究做出了重要貢獻。例如，適用于長期預報的判據[9]：穩(wěn)定性系數(K)、聲發(fā)射參數以及可靠概率(Ps)；適用于臨滑預報的判據[10]：位移加速度、位移矢量角、蠕變曲線切線[11]等。以上判據的提出對于我們在滑坡的理論研究上有重要意義，但是對于實際的預警預報工作卻有諸多不便，例如指標不易觀測、經濟成本高等。許強等[12]強調應加強區(qū)域性氣象預警和單體滑坡預警工作，將預警的重心轉移到對歷史數據的統(tǒng)計分析和基于變形、地下水位、雨量等關鍵指標的預警模型和判據研究。Shruti Naidu等[13]利用事件進行回歸分析，建立降雨閾值方程，通過降雨閾值分析滑坡穩(wěn)定性的Fs值建立安博里滑坡預警系統(tǒng)，但是該法僅對淺層滑坡適用；陸研等[14]利用Liang Kleeman信息流方法研究得出重慶地區(qū)滑坡泥石流發(fā)生頻次與年均降雨量顯著相關，建議在暴雨發(fā)生時加強對滑坡災害的監(jiān)測與預報；吳益平等[15]基于信息—物元模型將傳統(tǒng)的滑坡災害危險性預警延伸到滑坡災害風險預警預報，為降雨誘發(fā)滑坡的預警工作提出了新的思考。綜上所述，目前國內外關于降雨型滑坡的預警研究大多屬于大面積區(qū)域上的預報，對于單體滑坡的降雨參數預警卻很少涉及。

文中依托貴州省地質災害自動化監(jiān)測預警項目對該坡體展開深入研究，并基于坡體表面堆積松散、內部兼具低滲透性的特性對其預警預報工作開展研究。通過現(xiàn)場調查和無人機圖像(UAV)結合數字高程模型(DEM)建立該滑坡的PFC3D數值模型，模擬該坡體在自重和不同地下水位靜水壓力作用下的變形和運動過程，監(jiān)測模擬過程中顆粒的位移和速度，將位移突變點(計算不收斂)對應的地下水位作為臨界參數；在Pradel & Raad提出的降雨入滲土坡深度的算法基礎上，反演臨界降雨時長以及降雨強度參數，并驗證降雨預警閾值的合理性。

1 研究區(qū)工程地質概況

研究區(qū)表層堆積體結構松散，部分降水透過松散土層入滲至薄層泥巖風化物。泥巖風化物級配曲線如圖1所示，土體顆粒主要為細粒組和砂粒組，粘粒含量約為21%，粉粒約占32%，砂粒組含量47%，基本符合含砂粘質粉土的粒徑組成。另外，土體曲率系數Cc=2.63，不均勻系數Cu=64.2為級配良好型土體，其顆粒級配較好在宏觀上表現(xiàn)為較小顆粒能夠很好的填充在空隙當中，表現(xiàn)出相對較低的滲透能力。

滑坡后緣、中部不斷有工程活動建設，各種填方堆載加大了斜坡下滑力(圖2)?；鲁醅F(xiàn)變形時間為2008年，于2013年4月開始，受強降雨作用影響，滑坡體上的房屋、道路、地面陸續(xù)出現(xiàn)開裂變形、沉降等變化。在多雨季節(jié)由于上部土體松散，降水入滲坡體內部；而坡體底部土體滲透系數低，造成水體易進難出，使得斜坡內部地下水位線持續(xù)升高，自重增加且孔隙水壓力上升，坡體的穩(wěn)定性朝著不利的方向發(fā)展(圖3)。從2017年10月25日第一次監(jiān)測開始，測點位移量較前期監(jiān)測有增大趨勢，滑體坡面上有多條新裂縫，據此推斷山體蠕動變形速度具有加劇趨勢，穩(wěn)定性較差，發(fā)生滑坡的可能性增大，嚴重威脅當地居民的生命財產安全。

圖3 納縫滑坡工程地質剖面圖(3-3剖面)Fig.3 Engineering geological profile of Nafeng landslide

2 離散元數值模擬

2.1 顆粒流計算原理

近年來，顆粒流(PFC3D)作為離散元軟件，已經成為邊坡工程、地下洞室開挖、路基穩(wěn)定等巖土工程領域中十分成熟的數值計算方法[16]。區(qū)別于有限元，顆粒流中的球體單元較為符合散體材料的特性，在涉及土質邊坡的計算中表現(xiàn)出一定的優(yōu)越性，因此被眾多學者沿用[17]。球體之間的接觸包括ball-ball接觸和ball-wall接觸。球體受到重力以及外力的作用失去原有的平衡，在顆粒的動態(tài)變化中通過牛頓第二定律計算其位移變化，通過不斷更新球體位移計算模型的運動情況。

2.2 微觀參數的獲取

在離散元的計算中，散體材料的宏觀力學特征通過微觀強度表現(xiàn)。但是，并沒有合適的函數關系可以合理地表示散體材料的宏觀參數與微觀強度之間的非線性關系。通常情況下，通過匹配數值試驗與室內試驗的應力-應變曲線獲得[18]。結合該研究區(qū)的勘查資料以及現(xiàn)場取樣，選擇室內大型剪切試驗進行力學測試(圖4)，不同法向力作用下的剪應力-剪切位移曲線見圖5，所得宏觀力學參數如表1，其中滲透系數是通過現(xiàn)場雙環(huán)滲透試驗測得。

圖4 試驗裝置 圖5 不同法向力下的剪應力-應變曲線

表1 土樣宏觀物理力學參數Table 1 Macro mechanical parameters of soil samples

為了更好的模擬土體中力的作用，選擇接觸粘結模型。試驗所得粘質粉土細顆粒粒徑小于0.05 mm，在模擬過程中，數值試驗的顆粒粒徑值選取很難完全與真實的顆粒粒徑保持一致，為了滿足計算效率需要，等比例擴大半徑生成顆粒。模型尺寸25×25×25 m，顆粒最大最小半徑比為1.66，通過內嵌fish語言伺服控制，實現(xiàn)大剪試驗三維數值模擬(圖6)。模擬過程中反復調整相應微觀參數，直到應力-應變曲線(圖7)逼近室內試驗結果。數值試驗所得微觀參數如表2所示。

圖6 剪切盒與數值試驗結果Fig.6 Shear box and shear results

圖7 大剪數值試驗剪應力-應變曲線Fig.7 Shear stress-strain curve of large shear numerical test

表2 滑坡模型微觀參數Table 2 Micro-strength parameters of the landslide PFC model

2.3 納縫滑坡(PFC3D)模型建立

納縫滑坡數值模型主要包括地形高程數據獲取以及幾何模型的生成。地形數字高程通過現(xiàn)場無人機勘測攝像(UAV)生成，利用該數字高程生成幾何模型(圖8)，模型的相關參數為已獲得的微觀強度參數。

圖8 納縫滑坡(PFC3D) 數值模型圖 圖9 滑體地下水靜水壓力計算示意圖

隨著降雨時長的增加，入滲到坡體內部的水也會隨之增多，而滑面以上的土層透水性能差致使這些水體不斷聚集抬高地下水位線，該處的靜水壓力值將不斷增大，在靜水壓力的作用下滑面滑動的趨勢增強，當地下水位線增加到一定程度時，土體內部的抗剪強度將無法承受變大的靜水壓力值而發(fā)生破壞(圖9)?；诖?，該模型在計算過程中，主要通過改變顆粒不同靜水壓力大小而全面分析滑坡各個位置的變形情況，模擬過程中分別對滑坡前緣、后緣以及中部9個區(qū)域進行全過程監(jiān)測。

2.4 模擬結果

模型在代入匹配的微觀參數后到達最初的平衡，此時對每個顆粒施加垂直向下的靜水壓力(p=ρhg)，每個顆粒在靜水壓力平行于坡面的分力作用下開始發(fā)生變形。通過反復調試靜水壓力的大小，找到模型位移不收斂的靜水壓力區(qū)間。隨著降雨入滲，土體顆粒的靜水壓力增大是一個動態(tài)增加的過程。當靜水壓力值在某個界限值以下時，土體顆粒最終會趨于平衡。在該模型的計算過程中，增加顆粒的靜水壓力值，顆粒的位移曲線會出現(xiàn)階梯狀態(tài)的增長，但是隨著靜水壓力值的增大，其平衡的程度會降低，具體表現(xiàn)為平衡階段曲線的水平程度降低。靜水壓力值增大到某一個界限值時，顆粒之間的力鏈聯(lián)結發(fā)生斷裂時，則認為坡體發(fā)生了破壞。

靜水位高度h=1.60 m時坡體監(jiān)測點的變形如圖10所示。在靜水壓力的作用下，坡體各個部位均發(fā)生不同程度的變形，變形速率先增加隨后減小最終趨于零，說明位移收斂，坡體趨于平穩(wěn)。由圖可知，監(jiān)測點1、監(jiān)測點4附近位移較大，其峰值分別為0.32 m和0.21 m；其余監(jiān)測點位置變形相對較小，位移變化在0.02 m-0.09 m之間。結合監(jiān)測點所對應位置分析，坡體左側的穩(wěn)定程度相對較低，當外部條件發(fā)生改變，其變形破壞跡象明顯；而監(jiān)測點3、6、9位移變化較小。

圖10 靜水壓力高度1.60m時監(jiān)測點位移-時步曲線Fig.10 Displacement-time curve of monitoring point at 1.60 m hydrostatic pressure

當靜水位高度增加到1.70 m時，隨著顆粒受到靜水壓力的增大，變形也呈現(xiàn)出增大的趨勢(圖11)，此時各個監(jiān)測點位置位移增長較為明顯。但是在整個坡面區(qū)域，位移變化卻存在較大的差異性，從各個監(jiān)測點的位移變化來看，監(jiān)測點1、4所在位置為整個坡面最不穩(wěn)定區(qū)域，監(jiān)測點3穩(wěn)定程度次之，監(jiān)測點2穩(wěn)定程度最高。整個坡面的后緣位移最大、中部以及前緣位移依次遞減，符合推移式滑坡特點。在該靜水壓力條件下，整個坡面的位移最終也是趨于穩(wěn)定，但是穩(wěn)定程度明顯降低，當外營力增加時極有可能誘發(fā)滑坡。

圖11 靜水壓力高度1.70 m時監(jiān)測點位移-時步曲線Fig.11 Displacement-time curve of monitoring point at 1.70 m hydrostatic pressure

當靜水位高度繼續(xù)增加到1.75 m時，各個監(jiān)測點的位移明顯變化，局部位置的位移表現(xiàn)出不收斂的變化特性，如圖12中監(jiān)測點1、3、4?？傮w來看，滑坡的失穩(wěn)發(fā)生在其后緣、中部位置，而滑坡前緣的位移變化均能趨于穩(wěn)定(圖13)。分析認為，這是由于上部顆粒在水壓力作用下向下滑落，在前緣位置堆積所致。在該靜水壓力條件下顆粒之間的力鏈聯(lián)結主要表現(xiàn)為拉應力和壓應力，數值模型破壞時的力鏈分布如圖14所示。圖中藍色代表壓應力，綠色代表拉應力，剖面3-3’監(jiān)測點1、4、7所對應位置顆粒之間的力鏈主要為壓應力，局部位置為拉應力，斜坡表現(xiàn)為整體下滑；而剖面1-1’監(jiān)測點3、6、9所在剖面在后緣位置有著明顯的張拉破壞，土體顆粒向下運動出現(xiàn)顆粒的堆積；在滑坡的中部位置以土體的剪切破壞為主，有明顯的剪切裂縫產生，并伴有部分張拉-剪切復合破壞，破壞沿著垂直于坡面方向逐漸發(fā)展貫通程度逐漸加大，最后裂縫貫通的區(qū)域出現(xiàn)整體的下滑并進一步擠壓前緣位置顆粒；監(jiān)測點2、5、8的位移變化情況可以看出斜坡2-2’剖面在該靜水壓力條件下的位移是收斂的，最終仍然趨于穩(wěn)定，在其中部顆粒最大位移是25.4 cm，分析認為可能是由于1-1’剖面的張拉-剪切復合破壞擴展發(fā)展所致，但是并未形成貫通性的破壞，總體表現(xiàn)出較穩(wěn)定的狀態(tài)。

圖12 靜水壓力高度1.75 m時監(jiān)測點位移-時步曲線Fig.12 Displacement-time curve of monitoring point at 1.75 m hydrostatic pressure

圖13 位移云圖Fig.13 Displacement cloud map

圖14 力鏈分布圖Fig.14 Force chain distribution

以上過程表明，納縫堆積體滑坡穩(wěn)定性與入滲坡體形成的靜水壓力大小關系緊密。經過計算，靜水位高度為1.60 m時，坡體位移最終會處于穩(wěn)定狀態(tài)；當靜水位高度增加到1.70 m時，坡體位移曲線平緩程度下降；繼續(xù)增加靜水位高度到達1.75 m時，部分位移不能收斂，其強變形區(qū)為N01、N03以及N04。通過現(xiàn)場無人機勘測攝像(UAV)生成數字高程模型(DEM),建立的納縫滑坡PFC3D數值模型能夠全方位、多角度地反應納縫滑坡的變形破壞過程，其位移軌跡能夠較好的吻合現(xiàn)場監(jiān)測影像。相比較PFC2D，PFC3D具有明顯的優(yōu)越性，更加真實地反應滑坡過程中各個部位變形的強弱，更有利于在預警工作中對強變形區(qū)的識別和判斷，做出合理的預警決策。

3 滑坡預警閾值

3.1 反演模型

納縫堆積體滑坡具有低滲透特性，在降雨過程中斜坡表面的土體含有不同粒徑的碎石，因此表現(xiàn)出一定的松散特性，加之在不斷的蠕動變形過程中斜坡各部位尤其是后緣發(fā)育有多條裂隙，這為降雨入滲提供了便利的滲流通道，更多的降水通過這些通道進入到坡體內部。蠕動變形當中再次松動的是泥巖風化土體，其細粒組含量達到53%，大小顆粒之間能夠互相填充表現(xiàn)出良好的級配特性，滲透系數為3.01×10-5cm/s，透水性能較差。因此，入滲到土體內部的水體最終會在該處富集提高地下水位線，抬高的水位線高度可近似與入滲的深度相同，水體的富集會使孔隙水壓力上升有效應力降低，同時水位線以下部分土體的靜水壓力會逐漸升高從而誘發(fā)滑坡。

降雨對滑坡穩(wěn)定性的影響可以分為以下2種情況：一是當降雨強度小于巖土體入滲強度時，坡面大部分降水通過表面徑流作用排泄，少部分降水通過表面入滲作用浸入巖土體內部。此時巖土體內部必須考慮非飽和滲流問題，根據相應水文地質、工程地質條件，綜合考慮其他相關因素，如滲透性、植被、地下水位變化等情況；二是當降雨強度大于巖土體入滲強度時，水通過已有的滲流通道以及表面的入滲作用能夠不斷浸潤坡體，增加土體自重從而誘發(fā)滑坡。上覆土層的飽和區(qū)高度受降雨強度影響，并隨著降雨時長增加而升高。本文僅針對后者，討論當降雨強度大于巖土體入滲強度時，通過離散元數值模擬所確定的臨界靜水位高度進而推算出降雨強度以及降雨時長，作為納縫滑坡的預警參數閾值，可為納縫滑坡自動化預警提供可用依據。

強降雨工況下，土體入滲深度的推算需要滿足的2個基本條件是降雨強度Imin大于土體入滲率v0并滿足一定降雨時長(T)。目前針對降雨土坡入滲深度的計算方法較為合理的是國外學者Pradel & Raad[19]提出的經驗算法。本文依托離散元計算出的靜水位高度，據此公式反演出相對應的臨界降雨強度以及降雨時長，并計算出降雨總量。計算方法如下所示：

(1)

(2)

(3)

式中，v0為入滲率；k為滲透系數(cm/s)；h為入滲深度(cm)；t為降雨時長(s)；S為基質吸力(cm)；θs為土體飽和含水率(%)；θ0為土體實測含水率(%)；ds為相對密度，e為孔隙比。

基質吸力通過粘質粉土的土水特征曲線[20]得到，計算方法如下:

θ0×100=-0.0002S2-0.0115S+46.04.

(4)

式(3)反映了土體體積含水率與基質吸引力的動態(tài)變化關系。本文研究土體在降雨過程中的臨界閾值問題，隨著降水增多土體含水量增加，因此在使用式(3)的時候需要代入土體破壞時的含水量。同時，在降雨過程中土體水量增加土體顆粒有一個從塑限到液限的轉變過程。當土體顆粒趨于液限時，土體就由可塑態(tài)變?yōu)榱鲃訝顟B(tài)，土體也就喪失了其應有的強度?；谝陨蟽牲c原因，在通過式(3)推導其基質吸力時，θ0的取值為土體液限值，即θ0=wL=39.5%。

將θ0=39.5%代入到式(4)求解方程組得：

基質吸力S=154.352 cm

降雨強度：

降雨時長：

累積總降雨量Q=I×T=5.73×10-5×8.87×24×3600=439.13mm

通過以上計算得出納縫滑坡臨界降雨強度I為5.73×10-5cm/s，降雨時長T為8.87 d，累計總降雨量為439.13 mm。

3.2 滑坡預警閾值指標擬定

根據滑坡預警指標劃分方式采用閾值的60%、70%、80%和90%進行預警指標劃分[21]，分別對應藍色預警、黃色預警、橙色預警和紅色預警，各等級判別的臨界位移見表3，降雨強度和降雨總量見表4。

表3 納縫滑坡位移預警指標Table 3 Early warning of displacement of Nafeng landslide cm

表4 納縫滑坡降雨預警指標Table 4 Rainfall early warning of Nafeng landslide

3.3 實測資料驗證

羅甸縣氣象局2008年至2017年10年降雨量統(tǒng)計，年最大降雨量1 326.2 mm(2014年)，年最小降雨量846.1 mm(2009年)，年平均降雨量1 098.68 mm。羅甸縣月平均降水量如表5所示，從降水量的時間分配來看，5月到10月是多雨的季節(jié)，集中了全年降水的89%。

表5 研究區(qū)2008-2017年每月平均降雨量Table 5 Average monthly rainfall in the study area, 2008-2017

研究區(qū)于2017年10月25日開始布設觀測點等基礎工作，并于10月25日進行第一次觀測點測量。截止2018年9月29日，對該滑坡體進行了20次監(jiān)測，根據各監(jiān)測點情況，最大累計水平位移量為310.3 mm(圖15)，對應剖面3-3’所在區(qū)域，其次是1-1’剖面區(qū)域，變化最小的是2-2’區(qū)域；從整個坡體來看，位移變化量最大的地方是其后緣公路所在位置，坡腳處位移相對較小。

通過數值計算的結果(圖12)與實測數據(圖15)以及近10年的降雨數據(表5)分析發(fā)現(xiàn)：斜坡在2017年10月到2018年5月的位移變化較小，而近十年的月平均降水最少的時間區(qū)間也是10月到次年5月，在降水充沛的月份里位移的變化則是相當明顯的，這說明了導致該研究區(qū)變形破壞的最主要因素是降雨。從表3擬定的滑坡位移預警指標分析來看，監(jiān)測期間最大位移為3-3’剖面附近的CJ10處310.3 mm，1-1’剖面最大位移為CJ6處190.5 mm，監(jiān)測值與計算值接近，并且達到藍色預警等級，證明了離散元數值計算得出的臨界位移變化量的合理性；由表5可知，羅甸縣每年5月、6月為全年降雨量最大季節(jié)，平均降雨量分別

圖15 位移監(jiān)測結果Fig.15 Displacement monitoring results

為213.27 mm，204.73 mm。而滑坡初現(xiàn)時間2008年5月的降雨量為288.9 mm達到了雨量的藍色預警值；研究區(qū)降雨充沛，為該區(qū)的低速緩慢活動提供了不利條件。這種低速緩慢運動使得內部巖土體的強度逐漸變弱，在持續(xù)的強降雨條件下都有可能誘發(fā)滑坡。2017年8月降雨量為246.4 mm，基本達到藍色預警雨量，而監(jiān)測位移結果表明其位移已經達到藍色預警值。

綜合理論計算結果和實際位移監(jiān)測以及氣象局所獲得降雨量結果來看，本文計算的降雨閾值具有一定的合理性。但是，工程地質條件的復雜多樣性和土體強度隨著含水率變化的非線性等使得滑坡預警工作在監(jiān)測過程中存在很多困難。因此，可加強宏觀跡象特征的群測群防來提高監(jiān)測預警工作保障。

4 結語

(1)納縫坡體裂縫發(fā)育為降雨提供便利的滲流通道，坡體內部土體粘粒含量約21%，粉粒約32%，砂粒組含量47%，級配良好，滲透特性較低。分析認為：坡體表面開裂，易于入滲水體不斷在坡體內部聚集，土體有效應力降低，并伴隨滑帶處不斷增大的靜水壓力是滑坡變形的內在原因。

(2)對斜坡的PFC3D數值計算，證明了靜水壓力的大小與坡體的穩(wěn)定性關系緊密：隨著靜水壓力增大坡體穩(wěn)定性逐漸降低，滑坡破壞臨界靜水壓力所對應的臨界靜水位高度在1.70-1.75m之間。

(3)通過降雨入滲深度反演靜水位高度所對應的降雨強度以及降雨時長，最終確定縫滑坡預警參數閾值：降雨強度I=5.73×10-5cm/s，降雨時長T=8.87d，累積降雨總量Q=439.13 mm，并進一步通過實測數據驗證了預警參數閾值的合理性。