基于淺水流模型的廣州市大坪泥石流預測與防治

鄭 偉，尹 培 杰

(1.深圳高速工程顧問有限公司，廣東 深圳 518000； 2.長安大學 公路學院，陜西 西安 710000)

0 引 言

中國國土面積廣闊，地貌類型多種多樣，但同時也使得各種地質災害頻發。泥石流作為一種山區常見的自然地質災害，對道路交通、村莊和其他基礎設施會造成極大危害[1-4]。夏季隨著雨水增多，使得原本松散的物源在水動力的驅使下快速運移，沿途侵蝕鏟刮溝道岸坡，形成泥石流并最終堆積在下游，對沿途造成不同程度的影響。及時準確預測泥石流的危害范圍和危害程度，為相關部門提供及時的災害預警具有重要的意義。

通過經驗公式計算法或現場物理尺度模型等傳統方法來定量評價泥石流的影響相對困難，且精確程度受人為判斷影響較大[1，5-6]。上述方法在預測泥石流的過程中往往缺少動態流動過程，模型的尺寸效應使得結果和現場吻合度較低。隨著計算機技術的高速發展，科研人員對復雜過程的求解能力顯著提高。數值模擬方法因成本低廉且結論可靠，成為泥石流精確計算的一種重要方法[7]。多種數值模擬方法，例如淺水流模型[4]、平滑粒子流體動力學模型[8]、離散元模型[9]、兩相流甚至多相流模擬[10]，都已應用于泥石流模擬。

淺水流模型因其高效性成為科研和工程領域應用最廣泛的一類計算模型。劉鐵冀等[1]利用Massflow模擬攔擋壩對泥石流溝的治理情況，認為其能較好地控制溝內不穩定物源的活動情況。潘青等[2]利用離散單元的方法，通過數值模擬手段分析了泥石流運動過程對攔擋結構的沖擊破壞作用，取得了良好的效果。張奮翔等[3]利用FlO-2D軟件進行數值模擬，發現在允許條件下，將泥石流流速和流深與其發生頻率相結合，得到危險評價圖可為災害預測提供良好的建議。Han 等[4]基于淺水流計算開發了SFLOW軟件對泥石流和滑坡進行數值模擬，對比FlO-2D其對沖擊范圍模擬有著較好的結果。丹羽論等[5]依據二維的淺水流方程，研究了泥石流溝谷中設置的攔擋壩(谷坊壩)對泥石流的控制效果，并開展了數值模擬研究。

不同降雨頻率下泥石流的沖出范圍是泥石流預測的重點，而攔擋壩不同高度和位置對泥石流的減災效應仍待研究。廣東省部分山區長期受到泥石流災害困擾[11-14]，由此造成的經濟損失達數十億元。本文選取廣州市中部山區大坪溝為對象，采用淺水流模型對流域內泥石流的沖出范圍展開數值模擬，分析其運動過程并提出相應的治理方案。

1 研究區概況

大坪溝位于廣州市從化區，屬于廣東省中部地區，溝口坐標為東經113°46′23.5″，北緯23°40′44.14″(見圖1)。研究區植被較為豐富，多灌木。根據走訪當地居民可知，溝內多次發生洪水和泥石流，最近一次為2014年5月份，從化區普降大雨，大坪村數座房屋倒塌，農田被沖毀，公路多處中斷且沿途多處塌方。據統計(http：∥roll.sohu.com/20140524/n399975660.shtml)，此次從化區溫泉鎮和良口鎮是暴雨的重災區，死亡4人，失蹤4人，約1.2萬人受災。派潭鎮8 h內錄得雨量291.6 mm，流溪河洪水已經超過50 a一遇標準。

圖1 廣州市大坪溝整體景貌及溝谷的坡度和坡向分布范圍Fig.1 The overall landscape of Daping gully in Guangzhou City and the distribution of slope and aspect

大坪溝溝谷流域面積約6.27 km2，主溝長度為4.85 km。左側發育兩條主要的支溝，主溝外側可見鴨洞水河道由此蜿蜒而過。整體上呈現西南高東北低的地勢形態，最高的海拔高程為1 168 m，最低海拔為98 m，整體高差為1 050 m，屬于低山山地地形。大坪溝的主溝整體走向為西北向，溝道陡緩相接，溝道平均坡度為10°左右，局部坡度可達20°。除溝道坡度外，大坪溝兩側岸坡的坡度多在10°～30°之間，兩側山地稍顯陡峭，坡度高于30°，局部山坡坡度甚至可達60°，這種岸坡坡度為泥石流匯聚水能提供了有利條件。此外，現場調查結果顯示大坪溝內松散物源主要包括坡面崩滑體、溝道往期泥石流堆積物和人工堆積體，這些物源呈松散堆積，密實度較低，在雨水的沖刷之下可能進入溝道參與泥石流活動。

研究區地處中國大陸的最南端，屬于亞熱帶季風氣候。降雨量較高且降水相對集中，每年4～9月的降水量占全年的85%以上，此時是泥石流的暴發高頻期。研究區附近觀測的多年的降水量平均值為1 696.5 mm，日最大降水量為423.9 mm (24 h)[12]。

2 淺水流模型

2.1 理論基礎

二維淺水流方程被應用于描述水平尺度遠大于垂直尺度的流體，例如對泥石流、洪水等現象的描述。由于三維不可壓Navier-Stokes方程是非線性方程，求解過程比較復雜，可以通過假設所計算的流體在其垂直方向上滿足靜水壓力并且速度統一，將上述三維方程降為二維深度上平均、質量及動量上平衡的方程形式，即二維淺水流方程。其運算過程相對簡單且仍保持著良好的準確性，因此在本研究中，將二維淺水流方程作為泥石流運動模型的控制方程。其方程的整體形式為[5-6]

(1)

其中：

q=(h，qx，qy)T= (h，uh，vh)T

(2)

f=(uh，u2h+1/2gh2，uvh)T

(3)

g=(vh，uvh，v2h+1/2gh2)T

(4)

(5)

式中：q為變量失量；f為x方向上的通量矢量；g為y方向上的通量矢量；S為源相矢量；q為平面x、y方向的單寬流量；u，v為平面x，y方向相應的流速；h為流深；zb為溝道海拔高程。

泥石流流體運動特征和普通水體相差較大。運動過程中與邊界的摩擦行為較為特殊。基于不同的側重點提出了多種泥石流的摩擦模型，比如Coulomb摩擦模型、Voellmy摩擦模型、O’Brien 等提出的流變摩擦模型，所得結果也存在一定差異。綜合比較幾種模型，流變摩擦模型包含了泥石流流體的摩擦性質、黏度性質及固體顆粒接觸能量損失等，已經被廣泛用于處理泥石流流體的摩阻行為[6]，其表達式為

(6)

(7)

式中：τ為屈服應力；ρm為泥石流中的固體物質密度；K為阻力系數；β為泥石流的黏度；ntd為等效曼寧阻力系數。

2.2 數據獲取和參數設置

以網絡下載的SRTM DEM為地形數據，通過重采樣將其精度轉化為8.75 m。利用ArcGIS軟件柵格裁剪工具獲得研究區范圍的柵格數據文件，并將其轉換為ASCII碼文件帶入程序。泥石流入流點對于數值模擬較為關鍵，通過實地野外調查和遙感解譯可以確定研究區泥石流啟動點[3，6]。本次模擬選擇起點為泥石流溝形成區的下游點。

在模擬過程中設置流量過程曲線對數值模擬進行流量和時間控制?，F實中泥石流的過程曲線是復雜的且多具有脈動性，與泥石流過程呈陣發性有關。泥石流按動力來源分為暴雨型泥石流和冰川融雪型泥石流，后者多分布于高山冰川地帶。本文研究區處于熱帶地區，屬于暴雨型泥石流，通常認為此類泥石流的暴發頻率和降水頻率是相互對應吻合的，可將洪峰流量簡化為三角形或五邊形進行計算，認為是單峰值過程[14]。此次研究將其概化為五邊形模型，通過計算一次泥石流暴發的幾個關鍵點，以一次泥石流時間的1/3段作為分界點，以峰值流量的1/4和1/3作為在這兩個時間節點下的泥石流流量[15]，繪制出泥石流暴發的流量過程線(見圖2)。

圖2 五邊形法的泥石流過程曲線[14]Fig.2 Debris flow process curve by Pentagon method[14]

在泥石流計算過程中，體積濃度決定著泥石流的基本性質和參數，也影響著泥石流運動狀態。體積濃度相對大小直接決定了泥石流的密度，從而影響了流體的黏性，控制著屈服應力和流變參數等水文、力學參數。對于泥石流的運動狀態，一方面，泥石流流體中含有著大顆粒的固體，彼此之間的鑲嵌作用形成骨架結構對泥石流運動速度有著重大影響，而且在運移過程中會呈現出不均勻的分布狀態；另一方面，泥石流的破壞力、侵蝕性等也受到固體物質含量的影響。體積濃度是指固體物質體積與泥石流總體積的比值。O′Brien等[16]通過總結不同條件下的泥石流體積濃度，給出一個經驗值為0.55～0.65。根據參考文獻[3-4]，本文選取Cv=0.55為泥石流的體積濃度。

(8)

式中：Vs是所含固體體積，Vw是所含液體體積。

2.3 泥石流洪峰計算

數值模擬過程中，泥石流洪峰流量的確定將會影響最終結果。不同的流域特征觸發泥石流的降雨閾值呈現出較大差異性，導致泥石流洪峰計算結果出現差異。因此，為了研究不同降雨條件下大坪溝泥石流的沖出規模及危害范圍，本文設計了6種降水頻率，分別為P=20%，10%，5%，2%，1%，0.5%。根據水文手冊計算洪峰流量，其計算公式[5，16]為

(9)

本文依據暴雨洪峰流量，利用雨洪法求取泥石流相應的峰值流量，并考慮泥石流攜帶一定比例的泥沙物質，對暴雨洪峰進行修正。針對溝道形狀對其修正，不同溝道形狀對應泥石流的洪峰相差甚大，對其運移的速度和侵蝕速率都有影響。引入清水洪峰修正和溝道堵塞修正系數，最終的公式為[5，16]

Qc=(1+φ)QpDc

(10)

式中：Qc為在降水頻率為P時，相應的泥石流峰值流量，m3/s；(1+φ)為清水洪峰修正系數，可以由體積濃度Cv計算得到；Dc為溝道的堵塞修正系數。

最終設計的6種不同降水頻率下的泥石流水文曲線如圖3所示。

圖3 在不同降雨頻率條件下的水文曲線Fig.3 Hydrological curves under different rainfall frequency conditions

3 計算結果

3.1 數值模擬結果

對泥石流的數值模擬結果用ArcGIS軟件進行處理，可以更加直觀地看到泥石流沖出范圍和危害情況。圖4為大坪溝泥石流流深的數值模擬最終結果。在降雨頻率為20%，10%，5%時，泥石流運動路徑主要為沿大坪溝溝道運移，堆積厚度變化不大，最大厚度分別為2.1，2.2，2.7 m，一般沿溝道的泥石流堆積僅在1.5 m左右。而在降雨頻率小于5%時，也就是降雨重現期超過50 a一遇時，泥石流運移至外側鴨洞水河道，在溝口形成堆積扇，堆積厚度明顯增加，最大厚度增至5.3 m。當降雨頻率降低至0.5%時，泥石流在溝口堆積長度近800 m，而且向下游方向堆積范圍更廣，符合泥石流的運動特征?？紤]到大坪溝內暴發特大型泥石流時可能會堵塞鴨洞水河道，導致河水暫時性斷流，該情況下下游沿途村莊和道路都受到巨大的威脅，建議相關部門建立泥石流的監測和預警預報系統。但由于泥石流的堆積厚度不大，數值模擬結果為5 m左右，在河水作用下應該可以自然下泄。

在流速分布圖5中，可見泥石流沿溝道保持高速運移，最大流速為3.3～6.1 m/s，說明泥石流在溝道內動能較大，表現為侵蝕搬運，尤其是對溝岸兩側和溝道的侵蝕。但是運移至外側鴨洞水河道流速明顯降低，泥石流在此則以堆積作用為主。在降水頻率為20%和10%時，泥石流的流速均為較低速，最大不超過4 m/s，與中國著名的泥石流案例相比，例如汶川震后泥石流[17]和蔣家溝泥石流[18]分別為11.9 m/s和8 m/s，相差較大，這主要和泥石流的性質和洪峰流量相關。該泥石流最大流速可達6.1 m/s，雖仍低于上述兩值，但此時泥石流已經具備較高的破壞力和沖擊力[19-21]。

泥石流的危險對象評價是泥石流災害模擬的最終目的[22-24]。在降雨頻率大于5%時，泥石流對大坪村的影響較小，泥石流主要沿溝道運移，運移2 km左右。由于大坪村房屋建設主要位于主溝右岸地勢較高之處，此時泥石流的危害性不大，在適當防護的條件下，對經濟建設和生命財產安全影響不大。當降雨頻率取值降低，也就是高于50 a一遇的暴雨時，泥石流的危險性明顯增加。泥石流沿途可能會沖毀居民的房屋，最終堆積在溝口對岸的公路上，導致當地的道路癱瘓交通受阻。而且泥石流堆積范圍也明顯增大，僅溝口的堆積范圍就達到0.3 km2，直接影響到下游村莊。

在數值模擬過程中，由于設置的集水點位于大坪溝左側小支溝上游，泥石流運移在支溝溝口出現明顯的堆積現象，也就是圖4和圖5中，物質剛運移時形成的“鼓包狀”堆積。由于支溝溝口地形開闊，適宜泥石流堆積。在實際泥石流發生過程中，降水在主溝形成區匯聚開始攜帶溝道內的物質進行運移。當其運動到下游溝道時，來自側面支溝的匯流會增大泥石流動力，在其溝口堆積現象不明顯，或呈小型的舌狀堆積。

圖5 6種不同降水頻率下泥石流的流速分布Fig.5 Velocity distribution of debris flow under six different precipitation frequencies

3.2 危險性區劃

以往的研究表明可以按照泥石流流深對評價結果進行危險級別劃分[3]，本文按照0～1.0 m劃分為低危險區，1.0～2.0 m為中危險區，2.0～4.0 m為高危險區，大于4.0 m為極高危險區。對處理的結果在ArcGIS中進行柵格統計分析，可以得出條形統計圖(見圖6)。

圖6 泥石流危險區評價和危害指數(HI)Fig.6 Debris flow risk area assessment and hazardous index (HI)

由圖6可知，在評價結果中，泥石流的大多數影響范圍均為低危險區。當設計的降雨頻率高于5%時，低危險區的占比都接近或超過80%，其余的中危險和高危險區均位于支溝溝口和部分溝道。這3種情況的模擬結果均未出現極高危險區，說明此時泥石流危險性較小。但當降雨頻率低于5%時，即降雨重復期分別為50 a，100 a和200 a一遇時，評價結果中出現極高危險區，所占比例分別為2.54%，7.59%，14.81%，對應的低危險區的面積明顯下降，由原來的80%左右下降至77.65%，50.63%，47.23%，并且中危險區和高危險區也是呈現增大的趨勢。在本文中，將高危險區和極高危險區占總體泥石流危害范圍的比例定義為泥石流危害指數(HI)，該指數大小可以用來表征泥石流對沿岸道路和居民危害的程度，即

(11)

對于模擬得到的泥石流影響范圍，在設計的6個降水頻率情況下，泥石流危害指數HI分別為3.71，3.10，3.94，6.57，26.38，30.68，對應圖6中的黑色折線顯示，呈逐漸遞增趨勢，說明泥石流的危害性越來越大。泥石流最終的分布范圍主要為外側的鴨洞水河道，尤其是降雨頻率為0.5%時。

3.3 攔擋壩治理

通過上述數值模擬的結果可見，泥石流對溝口村莊和道路均有著不同程度的影響。在降雨頻率為0.5%時，泥石流運移過程出現范圍較廣的極高危險區。以往研究表明，攔擋壩對泥石流的治理效果是明顯的，因此，建議在大坪溝溝內修建攔擋壩以降低泥石流對沿岸道路及居民的威脅。

本文采用地形重構技術以研究攔擋壩對泥石流的阻攔作用。將大坪溝溝道地形進行修改，具體操作是將地形柵格高程數據轉換為ASCII碼，以文本文件形式即可在溝道內提高相應的高程數據來模擬修建的攔擋壩。當攔擋壩壩高設置為10 m，即將原溝道高程最小ASCII碼增加10，其他值相應增大以保持“攔擋壩”頂部水平。然后進行數值模擬，所用的水文曲線完全根據降雨概率P=0.5%時所設置。攔擋壩設計位置和高度的不同會影響到最終的治理效果，本文在大坪溝內設置兩個攔擋壩的位置(A點和B點)，此外設計3種攔擋壩高度(5，10 m和15 m)，以探究其對泥石流的治理效果。

如圖7所示，最終的數值模擬結果表明，設置攔擋壩之后，泥石流最大流深主要分布在攔擋壩之上，形成一個小型的“水庫”。由圖7可見：A，B兩點設置5 m的攔擋壩治理效果均不明顯(圖中紅色曲線代表200 a一遇的泥石流范圍)，基本和原始情況一致，如圖4(f)和7(a)、7(d)所示。在A點設置15 m高的攔擋壩最終的治理效果最為明顯，泥石流在溝口外側的主要影響范圍被減小至不足原先的50%，且攔擋壩以下的泥石流流深全部為1.3 m以下，屬于中低危險范圍內，對沿岸的基礎設施和居民影響基本不大，達到預期治理效果。泥石流的最大流深受到攔擋壩高度影響，這主要是最深處集中于攔擋壩上游的位置。

利用3.2節提出的泥石流危害指數HI評價攔擋壩治理效果。在A點設置攔擋壩為5 m時，HI指數為27.72，設置10 m和15 m攔擋壩時，HI指數分別為15.81和11.39，較原本200 a一遇的HI指數(30.68)有明顯下降；B點設置攔擋壩時，HI指數對應為24.72，14.25和8.01，因此，據HI指數可見B攔擋壩優于A攔擋壩效果。但是由圖7可見，設置A攔擋壩后，其上游出現了大范圍的極高危險區和高危險區，導致HI指數明顯增大。在實際情況中，依據攔擋壩下游的危險性區劃求得HI指數才是正確評價的指標，因此，提出“有效”危險區，公式(11)可改寫為以下形式：

(12)

由上式計算得到在A點分別設置壩高為5，10 m和15 m攔擋壩時的有效HI指數為7.33，0，0。

對于大坪溝的設置攔擋壩的位置，本文在主溝A，B兩點進行研究發現，3種不同高度的攔擋壩在A點的治理效果均優于B點的治理效果，推斷這主要是地形因素所導致的。A點位于一支溝溝口，適合攔擋壩截流由主溝而來的大部分洪水和泥石流，如圖7(a)～(c)所示。但B點地形開闊，修建攔擋壩的治理效果較差。因此，增加設計一組數值模擬以驗證上述結論，如圖8所示，在A、B兩點設置兩級攔擋壩對溝內泥石流進行治理，攔擋壩高度分別設置為10 m和15 m。最終結果顯示泥石流經A阻攔之后，B攔擋壩能發揮的作用較為有限，基本與單獨設置A攔擋壩相近，說明A處設置攔擋壩能夠有效治理大坪溝的泥石流。

圖7 對大坪溝進行不同方案攔擋壩治理效果對比Fig.7 Comparison of dam treatment effects of different schemes for Daping gully

圖8 對大坪溝進行兩級攔擋壩治理效果對比Fig.8 Comparison of the effect of two-stage barrage treatment on Daping gully

綜上所述，根據泥石流沖出危害范圍、HI指數和攔擋壩合理選址，對于大坪溝泥石流可選擇在A點出建立15 m高的攔擋壩進行治理。

4 結論與分析

本文采用淺水流的數值計算方法對廣州市大坪溝泥石流進行數值模擬計算，基于6個不同的降水頻率特征，計算了相應的水文過程曲線，并提出8種攔擋壩方案進行災害治理的對比，分析表明：

(1) 泥石流的流深分布圖顯示，泥石流對大坪溝溝口居民和道路設施的影響最為顯著。設計降雨頻率為0.5%時，最大流深為5.3 m，可能沿外側鴨洞水河道堆積長約800 m的范圍。但在降雨頻率高于5%時，泥石流主要沿溝道運移堆積，對沿途影響甚微。

(2) 泥石流的最大流速為3.3～6.1 m/s，主要集中在溝道，此時泥石流表現為侵蝕作用。當泥石流運移至外側河道時，流速明顯降低，此時，泥石流表現為堆積作用。

(3) 按照最終的泥石流堆積厚度進行危險區劃分，表明當降水頻率設計值大于5%時，泥石流影響區域80%以上均為低危險區。但是當設計值高于5%時，出現極高危險區。

(4) 對泥石流溝設置8個攔擋壩的方案進行數值模擬，結果顯示，當在A點設置15 m攔擋壩時，主要的泥石流堆積量被控制在攔擋壩之上，此時下游均為中低危險區，有效HI指數為0。