泥石流的數值模擬及危險性評價

康永德，侯精明，石寶山，郭敏鵬，潘占鵬，于國強，洪增林

(1.西安理工大學，陜西 西安 710048；2.中國地質調查局干旱-半干旱區地下水與生態重點實驗室，陜西 西安 710054；3.陜西省地質調查院，陜西 西安 710054)

泥石流是一種破壞性極強的地表動力過程，對人類生命安全、財產及生產活動造成嚴重影響。此過程較為復雜，前人已在數值算法及模型構建方面取得了豐碩成果[1-5]，如離散元法、格子玻爾茲曼法、光滑粒子流體動力學法等數值模擬方法。上述方法各有優劣，離散元法和光滑粒子流體動力學法所建立的模型雖然能較好地模擬復雜環境中的泥石流運動，但是計算成本較高大。Savage 等[6-7]提出的方法可以很好地處理不連續問題且效果較佳、結果可靠。王純祥等[8]采用基于GIS 的深度積分二維數值模型模擬了泥石流的運動規律。Hungr 等[9]借助DAN3D 模型模擬了泥石流的運動過程。Peng 等[10-13]使用FLO-2D 軟件和PFC2D軟件對泥石流進行模擬，并進行了危險性評價。Ouyang 等[14-15]提出的基于深度積分的Massflow 數值分析法對泥石流運動規律進行了模擬，取得了理想結果。Han 等[16]對淺水方程與侵蝕模型進行耦合，模擬了泥石流的運動規律和流變參數。由于網格大小和邊界條件復雜等問題，泥石流模擬中尚存在諸多挑戰[17]。如果不簡化某些項，計算效率將受限制，特別是數值計算過程中的不收斂和發散現象[18]。

除了掌握泥石流的運動規律，提前預測也是非常重要的工作之一。前人已取得了大量的研究成果[9-11,13]，主要常用方法有層次分析法、信息量法、基于聯系期望的泥石流評價等[19-24]。其中，層次分析法基于背景要素的統計方法，根據地貌形態對影響因素的綜合評分進行分區，因分區結果主要依賴分區專家的經驗和專業背景的局限性，所得指標無明確的物理意義[25]；人工神經網絡法能處理復雜的線性關系，但是預測模型的學習樣本與訓練樣本和典型性密切相關，然而在泥石流危險性等級評價過程中選取訓練的典型樣本集是一個較為困難的問題；同理，支持向量機也存在上述問題及非線性高維模式識別問題[26-27]。

泥石流的形成是一個較為復雜的過程，數值模擬存在復雜的數據運算，計算時間長將會直接影響到災害預測和救援決策。雖然也有研究者利用多核多線程及分布式計算等方式，但效果提升不明顯，或者成本太高。因此，亟待從根本上解決泥石流數值模擬中的計算時間長的問題。本文在借助上述模型算法優勢、克服其缺點的基礎上，將GPU 通用計算機方法應用于泥石流數值模擬中，提出了一套顯卡加速的泥石流過程數值模型（即GPU Accelerated Surface Water Flow and Transport Model，簡稱GAST 模型），可提升數值模擬過程中的計算能力，有效縮短程序運行時間，提高模型整體效率。

1 模型基礎

西安理工大學水模擬及災害管理團隊開發的GAST 模型，借助Visual Studio 2013 和NVIDIA 推出的并行計算架構CUDA 平臺，利用GPU 求解復雜計算問題。本文的數值模擬計算均基于GAST 模型完成。

1.1 控制方程及參數取值

泥石流本質是兩相流，但考慮到泥石流的運移特征和流變模式，同時保證計算的簡潔高效，故而，不考慮泥石流流體密度的變化及溝道水沙侵蝕的情況下，通常可將其視為弱壓縮流體，在三維N-S 方程的基礎上遵循嚴格的物理準則，推導出二維控制方程[14]，矢量形式如下：

式中：x、y、t分別為笛卡爾空間坐標與時間坐標；q為變量矢量，含泥深h和單寬流量qx與qy；u、v為x、y方向上流速；F、G為x、y方向上通量矢量；S為源項矢量，包括降雨i、底坡源項及摩阻力源項；zb為河床底面高程；，δ為單位面積上的有效應力，φ為摩擦角。

泥石流發展過程較快，持續時間短，初始計算條件較為簡單，所有網格區域計算節點上認為無泥石流，以此為起始計算點，所在節點的水深和流速初始值均為零。計算時，泥石流入口處給定流量與泥深，運動過程可以通過流量與泥深的變化過程體現。只考慮泥石流運動停止的最簡單條件，計算時設定泥石流停止時的運動速度為0.000 1 m/s。

摩阻力計算采用O’Brien-Julien 提出的方法，如下所示：

式中：τ為屈服應力（賓漢應力）；η為動力黏滯系數；γm為泥石流單位重度；κ為泥石流阻滯系數；n為研究區曼寧系數。結合Woolhiser 建議[28]，取κ=228；η一般參照研究區歷史災害事件，進行流變參數的率定，參照王裕宜等[29]的研究取值；考慮到泥石流的運移特征和流變模式，且保證計算簡捷高效，由于泥石流在水槽中運動時邊壁將會對泥石流產生阻礙作用，而水槽比較順直，根據錢寧的研究結果取糙率n=0.05 s/m1/3；為粒徑小于0.05 mm 細顆粒的質量分數（小數表示）；P2為粒徑大于2 mm 粗顆粒質量分數（小數表示）；γv為黏性泥石流的最小重度，取17.6 kN/m3；γ0為泥石流的最小重度，文中取15.8 kN/m3。

1.2 模型求解

采用中心格式的有限體積法，將變量存于網格中心，即將式（1）在控制體Ω上積分，如式（5）所示。

進而，將其轉化為控制體邊界的面積分，如下式（6）：

此處，F(q)n為相應界面的通量，笛卡爾坐標中，式（6）繼續轉化為曲面積分如式（7）所示。

式中：FE、FW、FN和FS分別為網格單元（i，j）的東西南北4 個界面的界面通量。

對界面通量計算采用Godunov 格式的HLLC 近似Riemann 求解[30]，該格式在處理干單元時的功能要優于其他格式，可以較好處理復雜地形下干濕交替變化，具有較強的激波捕捉能力。為構造數值解的空間精度為二階精度，且為了避免在數值重構后水面梯度大的地方出現數值振蕩的現象，模型采用TVD-MUSCL方法進行數值重構[31]，表達式如式（8）所示。

式中：U為計算變量；ri=(Ui+1?Ui)/(Ui?Ui?1),ri+1=(Ui+1?Ui)/(Ui+2?Ui+1)；ψ為限制器函數，本模型采用min mod 限制器，其函數表達式如式（9）所示。

此外，時間步長采用兩階Runge-Kutta 格式保證其二階精度[32]。干濕動邊界處理方面，水深容差為0.000 001 m 以區別干濕網格單元[33]。摩擦阻力項的計算使用分裂點隱式法以提高穩定性[34]。計算過程中，Courant 為0.5，使用CFL 準則來預測下一個迭代新的時間步長。

1.3 GPU 原理及加速計算過程

當前，圖形處理器GPU 已成為一種多核多線程、可編程、并行化的處理器，通過統一計算機設備架構[35]（Compute Unified Device Architecture，CUDA）的開發平臺，為通用計算提供平臺。采用C++及CUDA 語言進行編程實現GPU 加速的并行計算過程，在此不再贅述。

2 驗證算例

2.1 潰壩型泥石流模擬

潰壩型泥石流一般過程為水流全面漫頂，緊接著壩體坦化，最后整體潰決形成泥石流。為驗證本文GAST 模型的合理性，選取Komatina 等[36]于1997 年的試驗作為算例。文獻中算例的水槽長4.5 m，寬1.5 m，為了模擬結果更優，在文獻[36]的基礎上適當擴大，將其長度按比例擴至135 m，寬擴至45 m 進行泥石流不同時刻的運動規律模擬。試驗漿體在水槽頂端由可動擋板阻隔。試驗開始時，抽出可動擋板使漿體自由下泄，模擬不同時刻的泥石流運動距離與速度。本文模擬時長為20 s，選取典型時刻t=0、1、5、7、9 和15 s共6 個時刻（圖1），模擬的泥漿運動規律基本與物理試驗運動規律吻合，向前移動的距離隨時間的增長關系與實測結果大致相同。潰壩型泥石流的破壞過程可分為3 個階段，t=0 時，處于靜止狀態；t=1 s 時，隨著水流的不斷滲入，發展形成剪切面，并且到達水槽的中部；t=5 s 時，泥石流快速移動，開始向水槽底端移動。一次泥石流事件中，水槽中的潛在物源在降雨條件下先失穩滑動，進入水槽后與原有物源一起再次被來流沖刷而啟動，進而將泥石流沖出水槽。另外，水槽內剛開始先為深紅色，逐漸向橘紅色和綠色及藍色分布過渡，說明水槽內物源在侵蝕一定時間后開始向下侵蝕，形成了一個完整的泥石流過程。

圖1 潰壩泥漿不同時刻的運動結果Fig.1 Movement results of dam break mud at different times

圖2 為泥石流數值模擬不同時刻的泥深分布。選取6 個典型時刻的斷面，t=0 s 時，處于靜止狀態，泥深在前緣漿體水平距離運動到18 m 時達到最大值5 m；t=5 s 時，隨著泥石流沿水平方向的移動，泥深呈現波動增厚的趨勢，即有堆積可能性，移動至45 m 時，泥深達到最大值1.2 m。

圖2 泥石流不同時刻泥深分布Fig.2 Mud depth distribution at different times of debris flow

圖3 為模擬不同時刻的速度分布。t=0 s 時，泥石流物源處于靜止狀態；t=5 s 時，隨著泥石流在水平方向的移動，速度逐漸呈波動增大的趨勢，水平距離為50 m 時速度達到最大值13 m/s，然后減小至0，此時泥石流滑動體積約為330 m3；t=7 和9 s 時，水平距離為82 m，速度分別達到12 和10 m/s；運動水平距離為110 m時速度減為0，泥石流移動到了水槽的末端，整個過程結束。

圖3 泥石流不同時刻速度分布Fig.3 Velocity distribution of debris flow at different moments

綜上所述，通過模擬泥石流漿體不同階段的運動現象，模擬的水平距離和泥石流的速度結果與實際物理過程基本一致，計算結果與實測值吻合良好，說明本文提出的GAST 模型可用于泥石流運動的模擬。

2.2 泥石流危險區預測

選取貴州省普灑村滑坡泥石流事件進行危險性模擬評價。該區海拔約為1 800～2 170 m，垂直高差約370 m，降雨豐沛，多年平均約為1 238 mm。由于強烈的地質構造運動、物理風化作用及外力侵蝕作用等內外動力共同作用，將基巖轉化為破碎巖體。另一方面，使原本密實的土體變為松散土體，導致原溝道本來的松散物源增加至一定數量后成為潛在的泥石流溝。極端氣候及人類活動等外動力作用也使得松散物源增加，增大泥石流發生概率。

降雨（圖4）是促使泥石流易發生的重要因素之一。通過GAST 模型輸入降雨數據和12 m 高精度地形模擬泥石流危險性易發區，結果如圖5 所示。根據資料顯示，此次易發生泥石流的危險區覆蓋有61 個煤礦點（圖5）。因此，必須考慮未來突發降雨階段易發泥石流危險區域對煤礦點的影響。

圖4 降雨量隨時間變化Fig.4 Rainfall over time

圖5 為0～1.90 h 降雨條件下煤礦點的風險等級。t=0 h 時，降雨事件開始，各個煤礦點的風險等級基本一樣，無風險。隨著降雨量的增加，從t=0.50 h到t=0.90 h，各個煤礦點的風險等級在逐漸增加；t=1.90 h 時，所有煤礦點中風險最大的有3 個區域，即圖6 中的A、B、C，詳見圖6 細部放大圖。這可能是由于地形差異大，降水分布不均，弱巖層、滑坡和溝壑型產生了大量未固結的碎屑。對于巖石碎屑，雨水主要通過垂直張力裂縫穿透斜坡和滑動面，然后分別在后緣滑動面和底部滑動面的垂直裂縫中產生靜水壓力和升壓，從而影響斜坡的穩定性，最終導致滑坡，進而在雨水的混合作用下，形成滑坡型泥石流。

圖5 不同時刻滑坡點風險區劃Fig.5 Landslide risk zoning map at different times

本文根據GAST 模型模擬結果，將其危險等級劃分為3 種，即低風險區（0～0.5），中風險區（0.5～1.5）和高風險區（1.5～2.0）。圖6（a）為中國地質信息網公布的實測結果、圖6（b）為GAST 模型模擬結果，不難發現，GAST 模型模擬的泥石流易發危險區和“中國地質信息網”調查公布的實測結果基本一致，從圖6（b）可以看出，A、B、C這3 處發生泥石流的風險等級最高，數值接近2，依次對應實測結果的主槽體、西支、東支，從其細部圖可知，風險等級均為中、高風險區，A處有2 個煤礦點比較危險，B處有4 個煤礦點風險較高，C處若發生泥石流時，下游有1 個煤礦點較危險。

圖6 模型模擬結果與實測結果Fig.6 Model simulation results and measured results

可見，基于GAST 模型的泥石流數值模擬為風險等級評價研究提供了一種新方法，更為諸如貴州、云南、甘肅等山區易發多發滑坡、泥石流地質災害等問題研究提供了新視角。

2.3 GPU 加速時間效率比分析

基于2.2 節計算可知，貴州省納雍縣普灑村泥石流易發區面積為2.68 km2，12 m 分辨率高精度地形應用GAST 模型進行模擬，時間為3 h。12 m 分辨率地形的計算網格數量為2 592 280 個，CPU 計算的時間為5 151.9 s，GPU 所用時間為352.21 s，本文將單核Intel CoreTMi7-7700 作為基準，加速比假設為1，通過對比CPU 和GPU 模擬時間發現，NVIDIA GeForce GTX 1080 型號GPU 的泥石流模擬計算效率相比較Intel CoreTM i7-7700 型單核CPU 的加速比約為14.62 倍。可見，GPU 計算能力遠超CPU。

3 結語

本文將GPU 計算能力強的優點應用于自主開發的GAST 模型，通過CPU+GPU 協同計算，最大效率提升模型的運算能力，在泥石流數值模擬及危險性預測等災害中做了初步分析，并取得了較好的效果，為災害預防和救災快速決策節省了大量時間。研究得出以下主要結論：

（1）基于GPU 加速技術的GAST 模型可以很好地模擬泥石流運動及預測危險區范圍。

（2）利用所提出的模型對貴州省普灑村發生的泥石流進行反演，發現數值模擬的危險區范圍與調查結果基本一致，說明本模型可用于對泥石流易發區域的預測。

（3）GPU 的加速性能優良，可以低成本、高效率地滿足大量運算的需求。

目前，泥石流試驗成本較高，試驗難度也很大，本文只是初步研究了泥石流理想水槽試驗的運動規律，泥石流的運動機理復雜，還需進行大量試驗數據的模型率定與驗證，這也是今后GAST 模型在泥石流數值模擬方面努力發展的方向。