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

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

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

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

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

1 模型基礎(chǔ)

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

1.1 控制方程及參數(shù)取值

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

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

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

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

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

1.2 模型求解

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

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

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

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

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

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

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

1.3 GPU 原理及加速計算過程

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

2 驗證算例

2.1 潰壩型泥石流模擬

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

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

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

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

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

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

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

2.2 泥石流危險區(qū)預測

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

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

圖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 個區(qū)域，即圖6 中的A、B、C，詳見圖6 細部放大圖。這可能是由于地形差異大，降水分布不均，弱巖層、滑坡和溝壑型產(chǎn)生了大量未固結(jié)的碎屑。對于巖石碎屑，雨水主要通過垂直張力裂縫穿透斜坡和滑動面，然后分別在后緣滑動面和底部滑動面的垂直裂縫中產(chǎn)生靜水壓力和升壓，從而影響斜坡的穩(wěn)定性，最終導致滑坡，進而在雨水的混合作用下，形成滑坡型泥石流。

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

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

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

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

2.3 GPU 加速時間效率比分析

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

3 結(jié)語

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

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

（2）利用所提出的模型對貴州省普灑村發(fā)生的泥石流進行反演，發(fā)現(xiàn)數(shù)值模擬的危險區(qū)范圍與調(diào)查結(jié)果基本一致，說明本模型可用于對泥石流易發(fā)區(qū)域的預測。

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

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