韓 浩，侯精明，金 釗

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

0 前 言

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

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

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

1 研究方法

1.1 二維水動力模型

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

(1)

S=Sb+Sf

(2)

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

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

(3)

(4)

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

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

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

(5)

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

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

(6)

其中

fE=fi+1/2, j

fW=fi-1/2, j

gN=gi, j+1/2

gS=gi, j-1/2

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

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

1.2 GPU 并行計算方法

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

(1) 對整個研究區(qū)域進行計算區(qū)域分解，將計算區(qū)域分解為與計算平臺的GPU個數(shù)相同的子計算區(qū)域，本研究區(qū)域分解為2個子計算區(qū)域。

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

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

圖1 多GPU并行計算流程

1.3 模型計算性能評價

從模型的計算精度和計算效率兩個方面對所提出的新型流域雨洪過程模擬方法進行評估。

(1)計算精度分析

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

(7)

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

(2)計算效率分析

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

(8)

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

2 結(jié)果與分析

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

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

“V”型流域常被國內(nèi)外學者開展雨洪過程模擬驗證分析[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。模擬是在均勻降雨條件下進行的，模擬的降雨歷時為5 400 s，模擬時長為12 000 s，均勻降雨強度為10.8 mm/h，網(wǎng)格分辨率為5 m。對于邊界條件的設置，流域溝道出口處是自由出流的開邊界，其他均為閉邊界。

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

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

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

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

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

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

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

2.2.1葫蘆河流域概況

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

圖4 研究區(qū)域

2.2.2模型基礎參數(shù)

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

(9)

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

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

2.2.3模擬加速計算效率分析

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

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

3 結(jié) 論

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

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

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