韓 浩，侯精明，金 釗

(1. 中國科學(xué)院地球環(huán)境研究所 黃土與第四紀(jì)地質(zhì)國家重點(diǎn)實(shí)驗(yàn)室,西安 710061；2. 西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利工程國家重點(diǎn)實(shí)驗(yàn)室，西安 710048；3. 陜西黃土高原地球關(guān)鍵帶國家野外科學(xué)觀測研究站, 西安 710061)

0 前 言

洪澇災(zāi)害是全球所有災(zāi)害中發(fā)生頻率最多的自然災(zāi)害[1]。近年來極端降水強(qiáng)度呈增加的趨勢變化，極端降水事件風(fēng)險(xiǎn)也呈上升趨勢，同時(shí)洪澇災(zāi)害與極端降水事件是緊密相關(guān)的[2]。極端降水事件突發(fā)性強(qiáng)、強(qiáng)度大且精確預(yù)報(bào)難，造成的損失更大。中國洪水災(zāi)害問題也十分嚴(yán)重，如近年來受到國內(nèi)外學(xué)者關(guān)注的河南省“7·20”極端暴雨事件，此次洪災(zāi)造成大規(guī)模農(nóng)作物損壞、房屋倒塌損壞，人員傷亡等[3]。山洪災(zāi)害是我國洪澇災(zāi)害中最嚴(yán)重的災(zāi)種，同時(shí)會(huì)伴隨著一些滑坡和泥石流問題，嚴(yán)重影響社會(huì)經(jīng)濟(jì)發(fā)展[4]。因此，暴雨事件致災(zāi)快且損失更加嚴(yán)重，通過研究洪水過程的規(guī)律和洪水預(yù)報(bào)的能力提升，可為洪水災(zāi)害風(fēng)險(xiǎn)評(píng)估和防洪應(yīng)對(duì)能力的提升有重要的科學(xué)和現(xiàn)實(shí)應(yīng)用價(jià)值。

目前關(guān)于暴雨洪澇災(zāi)害的成因分析、洪澇應(yīng)急管理、洪澇模擬仿真等方面的研究已取得了一定的研究進(jìn)展[5-7]。數(shù)值模擬技術(shù)是洪澇災(zāi)害模擬的強(qiáng)力技術(shù)手段，能系統(tǒng)地實(shí)現(xiàn)降雨徑流過程的模擬，也能精確地預(yù)測洪水的演變過程。目前關(guān)于雨洪過程的數(shù)值模擬方法，包含物理過程及水文過程機(jī)理的數(shù)值模型主要有水文模型、水動(dòng)力模型、水文水動(dòng)力耦合模型。水文模型不能全面地計(jì)算水力特征信息(流速、流量、水位等)，而且很難準(zhǔn)確地描繪出水流運(yùn)動(dòng)傳播的全過程，但水動(dòng)力模型在模擬具有復(fù)雜下墊面的地表產(chǎn)匯流方面具有很大優(yōu)勢[8]。水動(dòng)力模型主要包括一維、二維、三維水動(dòng)力模型。二維水動(dòng)力模型，相比一維模型計(jì)算時(shí)間長且需要更多的基礎(chǔ)資料，但能得到更精確的模擬結(jié)果[9]。三維水動(dòng)力模型需要計(jì)算量很大，更是需要進(jìn)行大量工作完成網(wǎng)格設(shè)計(jì)和模型校準(zhǔn)[10]。二維水動(dòng)力模型的理論已經(jīng)很成熟，國內(nèi)外學(xué)者也采用了科學(xué)的數(shù)值求解方法提高模擬計(jì)算精度[11-13]，通過實(shí)例驗(yàn)證，二維水動(dòng)力模型能夠滿足當(dāng)前雨洪過程模擬的精度需求。相比傳統(tǒng)的水動(dòng)力模型，耦合水文過程的二維水動(dòng)力模型具有較高精度和可靠性[4]能準(zhǔn)確地預(yù)測洪水發(fā)生過程的水力特征信息，包括水流速度、流量和水深變化。然而，大多數(shù)高效高分辨率雨洪過程數(shù)值模型面臨模擬計(jì)算耗時(shí)長的問題，應(yīng)用GPU并行計(jì)算技術(shù)為數(shù)值模型計(jì)算性能的顯著提升提供了技術(shù)支持。GPU技術(shù)主要被國內(nèi)外學(xué)者們用于開展地表徑流過程的數(shù)值模擬，利用GPU并行計(jì)算優(yōu)勢使得高分辨率大規(guī)模淺水流動(dòng)問題的水動(dòng)力模擬變得可行[14-16]。通過建立基于GPU加速技術(shù)的流域雨洪模型能實(shí)現(xiàn)雨洪過程的高效高精度模擬已被眾多學(xué)者所采用，驗(yàn)證結(jié)果均表明在模擬精度較高的基礎(chǔ)上，通過引入GPU加速技術(shù)能極大提升雨洪模型的計(jì)算性能[4，7]。然而，國內(nèi)學(xué)者開展GPU并行計(jì)算技術(shù)的洪水?dāng)?shù)值模型的研究工作較晚，目前現(xiàn)有單GPU數(shù)值模擬方法不能滿足大規(guī)模、精細(xì)化且具有大量計(jì)算單元的流域雨洪過程模擬。對(duì)于解決單GPU計(jì)算性能不足問題，除對(duì)數(shù)值模型算法進(jìn)行優(yōu)化外，可通過引入多GPU并行計(jì)算技術(shù)，提升模型的計(jì)算性能。

本文基于耦合水文過程的二維水動(dòng)力模型，采用多GPU并行計(jì)算技術(shù)，提出一種多GPU并行計(jì)算的流域雨洪過程模擬方法，量化分析單GPU和多GPU的雨洪過程模擬的計(jì)算效率，將多GPU的數(shù)值模擬方法實(shí)現(xiàn)大尺度流域與洪過程模擬，探討該方法模擬計(jì)算的性能優(yōu)勢，為洪水過程高性能模擬提供技術(shù)支撐。

1 研究方法

1.1 二維水動(dòng)力模型

二維水動(dòng)力模型的控制方程是二維淺水方程，控制方程的矢量形式可表示為：

(1)

S=Sb+Sf

(2)

公式(1)中:q為qx、qy、h的變量矢量，其中qx、qy分別表示x和y方向的單寬流量，m3/s;h是水深,m；f和g分別為x和y方向的通量矢量；S是源項(xiàng)矢量，包括底坡源項(xiàng)Sb和摩阻源項(xiàng)Sf，方程中不考慮科氏力、風(fēng)應(yīng)力等。

控制方程(1)的向量形式如下：

(3)

(4)

公式(3)～(4)中:zb為地形高程,m；h為水深,m，則水位η計(jì)算為η=zb+h；u、v分別為x和y方向的流速，m/s,則qx=uh和qy=vh；g是重力系數(shù)；Cf是槽底糙率系數(shù)，是由曼寧系數(shù)n和水深h計(jì)算得到：Cf=gn2/h1/3。

模型控制方程中考慮了水文過程，公式(4)中i表示降雨和下滲產(chǎn)生的源項(xiàng)。在計(jì)算區(qū)域的每個(gè)網(wǎng)格計(jì)算單元均涉及到降雨量、下滲的計(jì)算，開展暴雨雨洪過程的模擬時(shí)，可忽略蒸發(fā)。

采用Godunov型有限體積格式對(duì)控制方程(1)進(jìn)行離散求解，積分形式如下：

(5)

公式(5)中:Ω為控制單元的有限體積；q為qx、qy、h的變量矢量；f和g分別為x和y方向的通量矢量；S為源項(xiàng)矢量；t為時(shí)間變量。

對(duì)方程進(jìn)行離散，則每個(gè)網(wǎng)格單元(i,j)在新的時(shí)刻更新后的流量可表示為：

(6)

其中

fE=fi+1/2, j

fW=fi-1/2, j

gN=gi, j+1/2

gS=gi, j-1/2

公式(6)中:Δt為時(shí)間步長；n是時(shí)間步；(i,j)表示網(wǎng)格單元的索引；Δx、Δy分別是結(jié)構(gòu)網(wǎng)格單元在x和y方向的尺寸大?。籈、W、N、S分別代表網(wǎng)格單元邊界的4個(gè)方向；fE、fW、gN、gS分別表示網(wǎng)格單元的界面通量；Si ,j表示網(wǎng)格單元(i,j)在中心的源項(xiàng)。

此外， 采用二階MUSCL方法對(duì)對(duì)變量進(jìn)行空間二階重構(gòu)以提高計(jì)算精度[14]；基于HLLC近似黎曼求解器計(jì)算網(wǎng)格單元界面上的通量[4]；采用龍格-庫塔法對(duì)時(shí)間進(jìn)行推進(jìn)[14]；采用底坡通量法對(duì)底坡源項(xiàng)進(jìn)行處理[17]；利用半隱式法計(jì)算摩阻源項(xiàng)[7]；通過靜水重構(gòu)解決干濕邊界負(fù)水深問題[7]?；谝陨蠑?shù)值方法，本研究建立了耦合水文過程的二維水動(dòng)力模型，能實(shí)現(xiàn)雨洪過程的穩(wěn)健數(shù)值模擬。

1.2 GPU 并行計(jì)算方法

對(duì)于文中提出的基于多GPU并行的流域雨洪過程模擬方法，水文及水動(dòng)力過程的顯示格式計(jì)算更有利于采用GPU并行計(jì)算方法實(shí)現(xiàn)高性能模擬計(jì)算。對(duì)于本研究的計(jì)算平臺(tái)，包括兩個(gè)GPU，所采用的GPU均是NVIDIA公司的NVIDIA GeForce GTX 1080，CPU是Intel(R) Core (TM)i7-8700?；诙郍PU的雨洪過程模擬方法的實(shí)現(xiàn)流程如圖1所示。主要過程如下：

(1) 對(duì)整個(gè)研究區(qū)域進(jìn)行計(jì)算區(qū)域分解，將計(jì)算區(qū)域分解為與計(jì)算平臺(tái)的GPU個(gè)數(shù)相同的子計(jì)算區(qū)域，本研究區(qū)域分解為2個(gè)子計(jì)算區(qū)域。

(2) 采用CUDA流方法實(shí)現(xiàn)多個(gè)GPU進(jìn)行并行計(jì)算，解決多個(gè)GPU的數(shù)據(jù)通信。

(3) 最終將每個(gè)GPU中計(jì)算的結(jié)果數(shù)據(jù)合并，對(duì)結(jié)果進(jìn)行處理分析。

圖1 多GPU并行計(jì)算流程

1.3 模型計(jì)算性能評(píng)價(jià)

從模型的計(jì)算精度和計(jì)算效率兩個(gè)方面對(duì)所提出的新型流域雨洪過程模擬方法進(jìn)行評(píng)估。

(1)計(jì)算精度分析

本研究所采用的基于單GPU的二維水動(dòng)力模型已在發(fā)表的工作中得到驗(yàn)證[18]，模型具有較高的計(jì)算精度。故此處只需對(duì)單GPU和多GPU的數(shù)值模擬方法的計(jì)算精度進(jìn)行評(píng)估，采用相對(duì)誤差對(duì)模擬結(jié)果進(jìn)行比較分析，可采用如下的計(jì)算公式：

(7)

公式(7)中:e相對(duì)為采用單GPU和多GPU計(jì)算的模擬結(jié)果數(shù)值的相對(duì)誤差；ν2GPUs是基于多GPU的計(jì)算值，ν1GPU是采用單GPU的計(jì)算值。

(2)計(jì)算效率分析

開展流域雨洪過程模擬方法采用的計(jì)算平臺(tái)包括單GPU和多GPU，所提出的新型流域雨洪過程模擬方法是基于多GPU的，量化比較多GPU較單GPU并行計(jì)算的加速計(jì)算效率。采用相對(duì)加速比對(duì)計(jì)算效率進(jìn)行量化分析，公式如下：

(8)

公式(8)中:t1GPU、t1GPUs分別表示數(shù)值算法程序在單GPU、兩個(gè)GPU上運(yùn)行的時(shí)間，s;S相對(duì)是相對(duì)加速比。

2 結(jié)果與分析

2.1 “V”型流域的雨洪過程模擬加速計(jì)算效率分析

2.1.1“V”型流域的模擬結(jié)果分析

“V”型流域常被國內(nèi)外學(xué)者開展雨洪過程模擬驗(yàn)證分析[4,19]，“V”型流域幾何構(gòu)成如圖2所示。流域的坡面和溝道均是不透水的，坡面長、寬、坡度、曼寧系數(shù)分別為1 000 m、800 m、0.05、0.015 s/m1/3；溝道長、寬、坡度、曼寧系數(shù)分別為1 000 m、20 m、0.02、0.15 s/m1/3。模擬是在均勻降雨條件下進(jìn)行的，模擬的降雨歷時(shí)為5 400 s，模擬時(shí)長為12 000 s，均勻降雨強(qiáng)度為10.8 mm/h，網(wǎng)格分辨率為5 m。對(duì)于邊界條件的設(shè)置，流域溝道出口處是自由出流的開邊界，其他均為閉邊界。

圖2 “V”型流域幾何構(gòu)成

如前節(jié)所述采用單GPU的雨洪過程模擬方法具有較高的計(jì)算精度已在發(fā)表的工作中得到驗(yàn)證，故只需進(jìn)一步驗(yàn)證所提出的基于多GPU的雨洪過程模擬方法的計(jì)算精度。建立基于單GPU和多GPU的數(shù)值模型開展雨洪過程模擬，將多GPU和單GPU計(jì)算的水深結(jié)果進(jìn)行比較分析，單GPU和多GPU模擬的水深分布值(t=3 600 s)如圖3所示。由圖3可知，兩種情景下的水深分布值相同；通過計(jì)算單GPU和多GPU的“V”型流域雨洪過程模擬水深值的相對(duì)誤差，計(jì)算得到相對(duì)誤差為0，結(jié)果表明多GPU的模擬方法具有較高的計(jì)算精度。綜上分析，基于多GPU的雨洪過程模擬方法是可行且可靠的。

圖3 單GPU和多GPU的雨洪過程模擬計(jì)算結(jié)果對(duì)比

2.1.2雨洪過程模擬的加速計(jì)算效率分析

“V”型流域?yàn)? m地形分辨率時(shí)，網(wǎng)格數(shù)為6.48×104。由于網(wǎng)格較少，不能較好地反映所提出模擬方法的加速計(jì)算性能優(yōu)勢。在參數(shù)設(shè)置保持不變的基礎(chǔ)上，開展地形分辨率為1 m(網(wǎng)格數(shù)為1.62×106)、2 m(網(wǎng)格數(shù)為4.05×105)的“V”型流域的雨洪過程模擬。如表1所示，量化分析了兩種網(wǎng)格數(shù)工況下的多GPU計(jì)算較單GPU的相對(duì)加速比，加速比分別為1.57、1.27倍，表明所提出的方法具有明顯的加速計(jì)算效果。通過比較分析，計(jì)算單元的網(wǎng)格數(shù)越多，相對(duì)加速比越大，加速效率越高。在計(jì)算區(qū)域具有較多網(wǎng)格數(shù)的情形時(shí)，采用多GPU并行計(jì)算的雨洪過程模擬方法的加速計(jì)算效率越顯著。因此，所提出的流域雨洪過程模擬方法具有顯著的加速計(jì)算性能。

表1 “V”型流域雨洪過程模擬的計(jì)算效率結(jié)果

2.2 葫蘆河流域雨洪過程模擬加速計(jì)算效率分析

2.2.1葫蘆河流域概況

葫蘆河流域位于我國黃土高原地區(qū)，是渭河上游最大的支流，位于東經(jīng)105°05′～106°30′、北緯 34°30′～36°30′，面積約為1.07×104km2，區(qū)域地形如圖4所示。流域的北部和東部是六盤山山脈，南部是秦嶺山脈，西部是黃土高原[20]。葫蘆河流域年平均日照時(shí)數(shù)約2 300 h，年總太陽輻射量為544.284 kJ/cm2[21]。流域干流建有秦安和靜寧水文站，多年平均降水量為440.2 mm。

圖4 研究區(qū)域

2.2.2模型基礎(chǔ)參數(shù)

建立基于多GPU的葫蘆河流域雨洪過程數(shù)值模型，模型輸入的基礎(chǔ)參數(shù)包括地形數(shù)據(jù)、降雨數(shù)據(jù)、曼寧系數(shù)等。地形分辨率為30 m的地形數(shù)據(jù)來源地理空間數(shù)據(jù)云(http://www.gscloud.cn)，流域的計(jì)算網(wǎng)格數(shù)約為2.2×107。降雨數(shù)據(jù)采用設(shè)計(jì)暴雨，選用50年一遇的設(shè)計(jì)暴雨，降雨歷時(shí)為2 h，累積降雨量約為54 mm，采用平?jīng)鍪斜┯陱?qiáng)度公式，計(jì)算公式如下[22]：

(9)

公式(9)中:i為設(shè)計(jì)降雨強(qiáng)度，mm/min；P是重現(xiàn)期，a；t為降雨歷時(shí)，min。

為了量化所提出的模擬方法在開展大尺度流域雨洪過程模擬的加速計(jì)算性能，故對(duì)其他參數(shù)進(jìn)行概化處理，曼寧系數(shù)參考經(jīng)驗(yàn)值取值，河流區(qū)域的曼寧系數(shù)為0.033，其他區(qū)域的曼寧系數(shù)為0.06。由于是在場次暴雨條件下進(jìn)行模擬計(jì)算，忽略蒸發(fā)，設(shè)置河道以外區(qū)域的下滲為0 mm。

2.2.3模擬加速計(jì)算效率分析

采用多GPU并行計(jì)算的雨洪過程模擬方法，建立基于多GPU的葫蘆河流域雨洪過程數(shù)值模型，模擬重現(xiàn)期為50年的設(shè)計(jì)暴雨(2 h)條件下的葫蘆河流域雨洪過程。模擬總歷時(shí)為5 h，其中降水過程和退水過程為2 h和3 h。由表2所示，統(tǒng)計(jì)了在不同計(jì)算平臺(tái)上模型運(yùn)行時(shí)間，顯然在當(dāng)前設(shè)備條件下，由于網(wǎng)格數(shù)較多，內(nèi)存受限，基于單GPU的數(shù)值模型沒有完成計(jì)算過程；但所提出的多GPU模擬方法完成了葫蘆河流域的的雨洪過程模擬計(jì)算。由表2可見，基于多GPU的模擬計(jì)算時(shí)間為1.8 h，相對(duì)于設(shè)計(jì)模擬時(shí)間(5 h)的加速比約為2.8倍，表明多GPU并行計(jì)算能實(shí)現(xiàn)大尺度流域的雨洪過程模擬。因此，所提出的基于多GPU的數(shù)值模擬方法可滿足較大尺度流域雨洪過程模擬的需求。同時(shí)，該方法具有模擬精度高和計(jì)算效率快的優(yōu)勢，可拓展到城市洪澇過程的高性能模擬。

表2 模擬計(jì)算效率結(jié)果

3 結(jié) 論

本研究基于耦合水文過程的二維水動(dòng)力模型，引入多GPU并行計(jì)算技術(shù)，提出了一種基于多GPU并行的新型流域雨洪過程模擬方法，采用典型算例和實(shí)際算例對(duì)所采用的模擬方法進(jìn)行了驗(yàn)證和計(jì)算效率量化分析，形成結(jié)論如下：

(1) 模擬“V”型流域的雨洪過程中，通過比較分析基于單GPU和多GPU并行計(jì)算的模擬水深結(jié)果，表明多GPU雨洪過程模擬方法可靠可信；在量化分析地形分辨率為1(網(wǎng)格數(shù)為1.62×106)、2m(網(wǎng)格數(shù)為4.05×105)情景下，多GPU計(jì)算和單GPU的相對(duì)加速比分別為1.57、1.27倍，多GPU模擬方法計(jì)算優(yōu)勢明顯優(yōu)于單GPU。

(2) 建立基于多GPU的葫蘆河流域雨洪過程數(shù)值模型，量化分析葫蘆河流域雨洪過程模擬的加速計(jì)算效率，結(jié)果表明所提出的方法能實(shí)現(xiàn)大尺度葫蘆河流域雨洪過程的高性能模擬。