泥石流的數(shù)值模擬及危險(xiǎn)性評價(jià)

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

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

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

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

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

1 模型基礎(chǔ)

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

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

泥石流本質(zhì)是兩相流，但考慮到泥石流的運(yùn)移特征和流變模式，同時(shí)保證計(jì)算的簡潔高效，故而，不考慮泥石流流體密度的變化及溝道水沙侵蝕的情況下，通?？蓪⑵湟暈槿鯄嚎s流體，在三維N-S 方程的基礎(chǔ)上遵循嚴(yán)格的物理準(zhǔn)則，推導(dǎo)出二維控制方程[14]，矢量形式如下：

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

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

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

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

1.2 模型求解

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

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

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

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

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

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

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

1.3 GPU 原理及加速計(jì)算過程

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

2 驗(yàn)證算例

2.1 潰壩型泥石流模擬

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

圖1 潰壩泥漿不同時(shí)刻的運(yùn)動(dòng)結(jié)果Fig.1 Movement results of dam break mud at different times

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

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

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

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

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

2.2 泥石流危險(xiǎn)區(qū)預(yù)測

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

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

圖4 降雨量隨時(shí)間變化Fig.4 Rainfall over time

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

圖5 不同時(shí)刻滑坡點(diǎn)風(fēng)險(xiǎn)區(qū)劃Fig.5 Landslide risk zoning map at different times

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

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

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

2.3 GPU 加速時(shí)間效率比分析

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

3 結(jié)語

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

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

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

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

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