新型流域雨洪過程模擬方法研究

韓 浩，侯精明，金 釗

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

0 前 言

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

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

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

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是重力系數；Cf是槽底糙率系數，是由曼寧系數n和水深h計算得到：Cf=gn2/h1/3。

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

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

(5)

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

對方程進行離散，則每個網格單元(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)表示網格單元的索引；Δx、Δy分別是結構網格單元在x和y方向的尺寸大小；E、W、N、S分別代表網格單元邊界的4個方向；fE、fW、gN、gS分別表示網格單元的界面通量；Si ,j表示網格單元(i,j)在中心的源項。

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

1.2 GPU 并行計算方法

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

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

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

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

圖1 多GPU并行計算流程

1.3 模型計算性能評價

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

(1)計算精度分析

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

(7)

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

(2)計算效率分析

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

(8)

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

2 結果與分析

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

2.1.1“V”型流域的模擬結果分析

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

圖2 “V”型流域幾何構成

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

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

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

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

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

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

2.2.1葫蘆河流域概況

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

圖4 研究區域

2.2.2模型基礎參數

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

(9)

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

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

2.2.3模擬加速計算效率分析

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

表2 模擬計算效率結果

3 結 論

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

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

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