劉乾，張洋銘，2，萬定生*

網(wǎng)格化分布式新安江模型并行計算算法

劉乾1，張洋銘1，2，萬定生1*

（1.河海大學(xué) 計算機與信息學(xué)院，南京 211100； 2.南京銀行股份有限公司，南京 210019）（ ? 通信作者電子郵箱dshwan@hhu.edu.cn）

近年來，網(wǎng)格化分布式新安江模型（GXM）在洪水預(yù)報中發(fā)揮了重大作用，但在進行洪水過程模擬時，模型數(shù)據(jù)量與計算量巨大，GXM的計算時間隨著模型預(yù)熱期的增加呈指數(shù)增長，嚴(yán)重影響GXM的計算效率。因此，提出一種基于網(wǎng)格流向劃分與動態(tài)優(yōu)先級有向無環(huán)圖（DAG）調(diào)度的GXM并行算法。首先，對模型參數(shù)、模型構(gòu)件、模型計算過程進行分析；其次，從空間并行性的角度提出了基于網(wǎng)格流向劃分的GXM并行算法以提高模型的計算效率；最后，提出一種基于動態(tài)優(yōu)先級的DAG任務(wù)調(diào)度算法，通過構(gòu)建網(wǎng)格計算節(jié)點的DAG并動態(tài)更新計算節(jié)點的優(yōu)先級以實現(xiàn)GXM計算過程中的任務(wù)調(diào)度，減少模型計算中數(shù)據(jù)傾斜現(xiàn)象的產(chǎn)生。在陜西省大理河流域與安徽省屯溪流域?qū)μ岢龅乃惴ㄟM行實驗，在預(yù)熱期為30 d、數(shù)據(jù)分辨率為1 km的情況下，相較于傳統(tǒng)的串行算法，所提算法的最大加速比分別達到了4.03和4.11，有效提升了GXM的計算速度與資源利用率。

網(wǎng)格化分布式新安江模型；網(wǎng)格流向劃分；并行計算；有向無環(huán)圖；任務(wù)調(diào)度

0 引言

水文模型主要用于洪水預(yù)報，為防汛人員提供調(diào)度決策支持，在防洪減災(zāi)方面有著非常重要的意義。水文模型分為兩類：集總式水文模型和分布式水文模型［1］。集總式水文模型的計算過程中沒有考慮流域的空間差異性，采用一套統(tǒng)一、固定化的參數(shù)，在下墊面物理條件差異大、面積分布廣的流域，該模型的洪水預(yù)報精度并不理想；分布式水文模型則充分考慮了流域內(nèi)降雨空間分布不均以及下墊面物理參數(shù)差異對流域內(nèi)產(chǎn)匯流的影響［2］，預(yù)報精度較高。典型的分布式水文模型有：分布式新安江模型（Distributed Xin’anjiang Model）、TOPMODEL模型（TOPography-based hydrological MODEL）［3］和TOPKAPI模型［4］等。

分布式新安江模型按照流域劃分粒度又可分為：分塊型分布式新安江模型與網(wǎng)格化分布式新安江模型（Grid-based distributed Xin’anjiang hydrological Model， GXM）。GXM是基于數(shù)字高程模型（Digital Elevation Model， DEM）實現(xiàn)的、以新安江模型為基礎(chǔ)的分布式水文模型［5］。GXM將流域劃分為大小一致的高精度網(wǎng)格單元后進行水文預(yù)報計算，計算過程主要包括對每個網(wǎng)格單元計算蒸散發(fā)、產(chǎn)流、分水源以及根據(jù)上下游依賴關(guān)系逐網(wǎng)格進行坡地匯流和河道匯流演算［6］。因此GXM計算過程的耦合性很高，多采用傳統(tǒng)的串行計算方式，計算效率低，無法滿足快速預(yù)報的需求。

針對上述問題，本文提出了一種基于網(wǎng)格流向劃分與動態(tài)優(yōu)先級有向無環(huán)圖（Directed Acyclic Graph， DAG）調(diào)度的GXM并行算法，實現(xiàn)了GXM的并行計算與任務(wù)調(diào)度，提高了模型計算效率。本文的主要工作如下：

1）設(shè)計了GXM并行框架：首先從模型參數(shù)、模型構(gòu)件、模型計算過程三個角度分析GXM；然后對流域進行網(wǎng)格流向劃分，得到流域網(wǎng)格單元間的獨立或者依賴關(guān)系，并在此基礎(chǔ)上實現(xiàn)GXM的并行計算與任務(wù)調(diào)度。

2）提出了一種基于網(wǎng)格流向劃分的多層級可并行網(wǎng)格劃分算法，通過該算法得到可并行網(wǎng)格演算次序序列，從而實現(xiàn)GXM的并行計算，提高GXM的計算效率。

3）提出了一種基于動態(tài)優(yōu)先級的有向無環(huán)圖（DAG）調(diào)度算法，在計算過程中動態(tài)調(diào)整GXM資源分配，減少了GXM計算過程中數(shù)據(jù)傾斜現(xiàn)象的產(chǎn)生，提高了GXM計算資源利用率。

1 相關(guān)工作

分布式水文模型的并行化方式一般分為：空間上并行化、時間上并行化、時間和空間組合的并行化［7］。為了在計算過程中方便衡量模型并行化的性能，Liu等［8］提出了理論最大加速比（Theoretical Maximum Speedup Ratio， TMSR），即串行算法的計算時間與并行算法的計算時間的比值，并將TMSR作為分布式水文模型并行計算能力的評價指標(biāo)。Liu等［9］針對分布式水文模型FSDHM（Fully Sequential Dependent Hydrological Models），基于OpenMP（Open Multi-Processing）實現(xiàn)了一種依據(jù)流向?qū)⒛M單元分層的并行算法，并在河北省清水河流域?qū)υ摬⑿兴惴ㄟM行實驗驗證。在數(shù)據(jù)集分辨率為30 m，計算線程數(shù)為24的情況下，模型最大加速比達到12.49，并行算法有效提升了FSDHM的計算速度。秦澤寧等［10］針對WEP-L分布式水文模型匯流過程采用OpenMP設(shè)計了區(qū)域分解并行方法，有效提高了模型匯流過程的計算速度。目前對分布式水文模型并行化的研究大多數(shù)采用從空間上并行化的方式，通過將模型計算時相互獨立的模擬單元并行處理，取得了良好的應(yīng)用效果。筆者所在團隊的前期工作中提出了基于流向的分布式水文模型并行計算方法［11］，實現(xiàn)了分布式水文模型的標(biāo)準(zhǔn)構(gòu)建與并行化改造，本文在該研究的基礎(chǔ)上從空間并行角度對GXM進行并行化研究，提出相應(yīng)的GXM并行算法。

為了進一步提高分布式水文模型并行計算的效率，國內(nèi)外研究人員對于模型并行計算過程中的任務(wù)調(diào)度和資源分配方式進行了大量研究。Liu等［12］針對分布式水文模型提出一種可擴展的分布式水文模型的兩級并行化方法，該方法可以同時在子流域級和基本模擬單元級（例如網(wǎng)格單元）進行并行化。實驗結(jié)果表明，兩級并行化方法比僅在子流域級別并行化的方法具有更好的可擴展性，并且并行性能隨著數(shù)據(jù)量和子流域數(shù)量的增加而提高。秦澤寧等［10］為了實現(xiàn)分布式水文模型并行計算中任務(wù)分配的負(fù)載均衡，設(shè)計了基于貪心算法的優(yōu)化調(diào)度方法，有效提高了模型的并行效率。但是目前對GXM并行計算過程中的任務(wù)調(diào)度研究仍然不多。因此，本文提出一種基于動態(tài)優(yōu)先級的DAG任務(wù)調(diào)度算法，以實現(xiàn)GXM并行計算過程中的任務(wù)調(diào)度，提升GXM并行計算的資源利用率。

2 GXM并行算法

2.1　GXM并行框架

從模型參數(shù)、模型構(gòu)件、模型計算過程三個角度來對GXM進行分析，本文提出的GXM并行框架如圖1所示。

圖1　GXM并行框架

GXM參數(shù)主要包括集總式參數(shù)與以網(wǎng)格劃分的模型參數(shù)。集總式參數(shù)一般為水文流域模型系數(shù)，每個流域?qū)?yīng)一組參數(shù)；以網(wǎng)格劃分的模型參數(shù)包含了流域下墊面條件、流域狀態(tài)場與土壤含水量等信息，每個網(wǎng)格擁有獨立的參數(shù)且不同網(wǎng)格之間參數(shù)獨立。對上述參數(shù)采用NetCDF文件形式進行存儲，NetCDF具有很強的靈活性，特別適用于大量多維度數(shù)據(jù)的傳輸與存儲，具有讀寫效率高、靈活度高、壓縮性高、使用與平臺無關(guān)等特點［13］。

依據(jù)功能差異將GXM中的計算任務(wù)劃分為四種類型的構(gòu)件，分別為蒸散發(fā)構(gòu)件、產(chǎn)流構(gòu)件、分水源構(gòu)件和匯流構(gòu)件，以實現(xiàn)網(wǎng)格單元之間的關(guān)系解耦。每種構(gòu)件內(nèi)部又可以細(xì)化為不同的計算組件，構(gòu)件劃分如圖2所示。

圖2　GXM構(gòu)件分類

針對各個構(gòu)件在計算過程中的特點，將構(gòu)件劃分為兩種類型：獨立型構(gòu)件和依賴型構(gòu)件［14］。

1）獨立型構(gòu)件：計算過程中，獨立型構(gòu)件的網(wǎng)格計算單元之間的關(guān)系相互獨立，無須考慮流域上下游計算單元間的依賴關(guān)系。獨立性構(gòu)件包括蒸散發(fā)構(gòu)件、產(chǎn)流構(gòu)件和分水源構(gòu)件。

2）依賴型構(gòu)件：計算過程中，依賴型構(gòu)件需要考慮網(wǎng)格計算單元在時空間上的依賴性，根據(jù)計算單元上下游依賴關(guān)系計算。依賴型構(gòu)件包括匯流構(gòu)件。

2.2　構(gòu)建網(wǎng)格流向矩陣

GXM依賴高分辨率的流域下墊面信息，通過提取流域內(nèi)網(wǎng)格單元的流向信息，確定上下游網(wǎng)格單元間的流向關(guān)系，從而構(gòu)建流域內(nèi)的網(wǎng)格流向矩陣。

采用D8算法［15］確定網(wǎng)格單元水流方向，算法思路為：設(shè)定每個網(wǎng)格單元水流方向只有8種可能，即只流入與之相鄰的8個網(wǎng)格中的一個網(wǎng)格，網(wǎng)格單元的8個流向用不同的數(shù)字表示，網(wǎng)格單元流向標(biāo)記如圖3所示。

圖3　網(wǎng)格單元流向標(biāo)記

在3×3的網(wǎng)格窗口中，首先計算中心網(wǎng)格與各相鄰網(wǎng)格間的距離權(quán)落差，取距離權(quán)落差最大的網(wǎng)格作為中心網(wǎng)格的出流網(wǎng)格。網(wǎng)格間距離權(quán)落差計算公式如式（1）所示：

其中：表示需確定流向的網(wǎng)格單元高程，表示與之相鄰的網(wǎng)格單元高程；指代距離權(quán)重，對角線網(wǎng)格取，其他正向網(wǎng)格取1；為網(wǎng)格單元邊長。圖4為網(wǎng)格流向矩陣及對應(yīng)的水流流向矩陣示意圖。

2.3　基于網(wǎng)格流向劃分的GXM并行算法

2.3.1GXM過程分析

GXM以高精度網(wǎng)格作為模擬單元，有良好的預(yù)報精度。在GXM計算過程中，將每一個網(wǎng)格單元作為一個子流域處理，首先分別采用三層蒸散發(fā)模型、蓄滿產(chǎn)流模型及自由水蓄水庫結(jié)構(gòu)對每個網(wǎng)格單元進行蒸散發(fā)計算、產(chǎn)流量計算、分水源計算，得到每一個網(wǎng)格單元的產(chǎn)流量與三水源；其次，按照流向順序依次進行網(wǎng)格單元坡地匯流和河道匯流演算；最后，將流域內(nèi)網(wǎng)格的地表徑流、壤中流、地下徑流演算至流域出口［16］。

2.3.2多層級可并行網(wǎng)格劃分算法

GXM因為網(wǎng)格單元的計算耦合度較高，所以傳統(tǒng)上多采用串行計算方法，按照流向順序逐網(wǎng)格將流量演算至流域出口，導(dǎo)致GXM計算時間長、資源利用率低等問題。因此，本文提出一種多層級可并行網(wǎng)格劃分算法，通過構(gòu)建流域的網(wǎng)格流向矩陣、累積匯水面積矩陣和水系矩陣，分析挖掘出流域網(wǎng)格間的獨立與依賴關(guān)系，找出流域內(nèi)可并行網(wǎng)格的并行演算次序序列，實現(xiàn)多層級可并行網(wǎng)格的劃分。

通過網(wǎng)格流向矩陣可以確定流域內(nèi)所有網(wǎng)格單元的上游網(wǎng)格，按照網(wǎng)格流向矩陣中的流向，逐網(wǎng)格演算至流域出口。演算過程中所經(jīng)過的網(wǎng)格的匯水量均增加一個單位，從而可得到流域的累計匯水面積矩陣。累計匯水面積矩陣存儲了流域內(nèi)每個網(wǎng)格單元的累計匯水量，矩陣中每一個數(shù)值代表了從上游匯流區(qū)流入當(dāng)前網(wǎng)格的網(wǎng)格單元數(shù)量。當(dāng)網(wǎng)格的集水面積超過某一閾值時，將該網(wǎng)格單元設(shè)定為河道網(wǎng)格，標(biāo)記為1；低于此閾值的網(wǎng)格設(shè)定為非河道網(wǎng)格，標(biāo)記為0。最后用這些標(biāo)記好的網(wǎng)格單元構(gòu)成河網(wǎng)水系矩陣。通過流向矩陣生成累計匯水面積矩陣以及水系矩陣的過程如圖5所示。

圖5　累計匯水面積矩陣及水系矩陣生成過程

接著根據(jù)流域出流網(wǎng)格向上推演，對網(wǎng)格之間的演算關(guān)系進行解耦，通過多層級可并行網(wǎng)格劃分算法計算出流域內(nèi)網(wǎng)格的并行演算次序序列，從而實現(xiàn)GXM的并行計算。多層級可并行網(wǎng)格劃分算法描述如下：

輸入 網(wǎng)格流向，水系，累計匯水面積，累積匯水面積最大值，網(wǎng)格數(shù)，目標(biāo)網(wǎng)格；

輸出 返回每個網(wǎng)格的并行演算次序。

for←0 todo

所有網(wǎng)格賦初始值0；

for←0 todo

if

找出流域出口點，并賦值的演算次序為

（1）；

if1&&=

尋找流入該網(wǎng)格的上游網(wǎng)格；

返回的演算次序1；

；

break

for←0 todo

倒序排列；

while（?=0） do

if0&&=1

返回演算次序=（1）；

break

else if0&&≠1

尋找該網(wǎng)格匯入的下游網(wǎng)格；

返回的演算次序1；

break

return 返回每個網(wǎng)格的并行演算次序

設(shè)流域內(nèi)網(wǎng)格數(shù)量為，最大演算次序為，演算次序為，當(dāng)前次序重復(fù)的網(wǎng)格數(shù)量為C，它們之間的關(guān)系如式（2）所示：

河道網(wǎng)格的演算次序按照河道匯流方向依次遞增，并且非水系區(qū)域的網(wǎng)格單元在流域內(nèi)的占比很大，這些網(wǎng)格單元的演算次序往往較小并且分布較為密集。因為具有相同演算次序的網(wǎng)格單元計算過程中不需要進行數(shù)據(jù)交換，完全相互獨立，所以這些具有相同演算次序的網(wǎng)格單元可以并行計算，同時為不同層級演算次序的網(wǎng)格單元分配不同的計算線程，從而實現(xiàn)GXM的并行計算。各個網(wǎng)格單元間的數(shù)據(jù)通過內(nèi)存進行傳遞交互，調(diào)度方式基于先到先服務(wù)方式。在山西省大理河流域與安徽省屯溪流域?qū)ι鲜鎏岢龅幕诰W(wǎng)格流向劃分的GXM并行算法進行驗證，發(fā)現(xiàn)并行算法的計算效率相較于傳統(tǒng)的串行計算方法有著明顯的提升。

3 基于動態(tài)優(yōu)先級的DAG調(diào)度算法

3.1　數(shù)據(jù)傾斜

數(shù)據(jù)傾斜現(xiàn)象是指GXM在并行計算過程中，因為不同演算次序的任務(wù)量不同而造成的不同計算線程之間負(fù)載不平衡的問題。理想情況下，每個計算線程所分配的計算任務(wù)量接近，但在預(yù)熱期長、流域面積廣的情況下使用第2章提出的基于網(wǎng)格流向劃分的GXM并行算法進行模型計算時，極易出現(xiàn)計算任務(wù)分配不均衡、資源浪費的問題。GXM在并行計算時產(chǎn)生的數(shù)據(jù)傾斜現(xiàn)象如圖6所示。

圖6　數(shù)據(jù)傾斜現(xiàn)象

3.2　網(wǎng)格單元動態(tài)優(yōu)先級劃分算法

針對GXM計算中出現(xiàn)的數(shù)據(jù)傾斜現(xiàn)象，提出一種網(wǎng)格單元動態(tài)優(yōu)先級劃分算法。首先去除流域內(nèi)無效網(wǎng)格，其次依據(jù)并行演算次序在網(wǎng)格計算節(jié)點之間添加有向邊，從而形成由網(wǎng)格計算節(jié)點和有向邊所構(gòu)成的DAG，用于表達流域網(wǎng)格單元間的計算任務(wù)關(guān)系，如圖7所示，圖中的網(wǎng)格數(shù)值代表網(wǎng)格編號。DAG中的每個網(wǎng)格有若干個前置任務(wù)網(wǎng)格和一個后置任務(wù)網(wǎng)格。每個網(wǎng)格被認(rèn)為是一個計算節(jié)點，每個計算節(jié)點有自身相關(guān)參數(shù)、前置任務(wù)集與后置任務(wù)集。

圖7　流域有效網(wǎng)格單元的DAG

網(wǎng)格單元動態(tài)優(yōu)先級劃分算法首先會根據(jù)并行演算次序序列對計算節(jié)點的任務(wù)優(yōu)先級進行初始化，并在模型計算過程中動態(tài)調(diào)整網(wǎng)格計算節(jié)點的任務(wù)優(yōu)先級。通過定義式（3）、（4）計算節(jié)點的任務(wù)優(yōu)先級：

輸入 節(jié)點個數(shù)，計算節(jié)點優(yōu)先級，當(dāng)前節(jié)點，每個網(wǎng)格初始演算次序， 線程數(shù)當(dāng)前網(wǎng)格的下游網(wǎng)格坐標(biāo)，某節(jié)點位置；

輸出 剩余每個節(jié)點更新后的計算優(yōu)先級。

for←0 todo

if為葉子節(jié)點

［］1；

else if為非葉子節(jié)點

［］=；

forall （in （1）） do

repeat

until

更新后置任務(wù)網(wǎng)格優(yōu)先級［］1；

while（［］=［］） do

更新［］，節(jié)點在就緒隊列等待；

返回更新后的；

break

更新［］，節(jié)點在就緒隊列等待，節(jié)點更新為當(dāng)前計算節(jié)點；

返回更新后的；

break

return 剩余每個節(jié)點更新后的計算優(yōu)先級

3.3　基于關(guān)鍵路徑的DAG調(diào)度算法

根據(jù)上面提出的網(wǎng)格單元動態(tài)優(yōu)先級劃分算法，提出基于關(guān)鍵路徑的DAG任務(wù)調(diào)度算法來實現(xiàn)GXM并行計算中的任務(wù)調(diào)度。首先基于DAG計算任務(wù)的深度和動態(tài)優(yōu)先級，建立任務(wù)狀態(tài)隊列包括就緒隊列、運行隊列和已完成隊列。就緒隊列包含了當(dāng)前可執(zhí)行的任務(wù)，即所有不存在前置任務(wù)的節(jié)點和前置任務(wù)已經(jīng)執(zhí)行完畢的節(jié)點；運行隊列包含了正在運行的任務(wù)；任務(wù)執(zhí)行完畢后，進入完成隊列。其中，在任務(wù)進入就緒隊列前，賦予任務(wù)節(jié)點三種屬性類型：葉子節(jié)點、等待節(jié)點和關(guān)鍵節(jié)點。葉子節(jié)點為DAG中第一層節(jié)點，不包含前置任務(wù)，任務(wù)優(yōu)先級初始化為1，在任務(wù)開始執(zhí)行時就處于就緒狀態(tài)；等待節(jié)點是將要執(zhí)行的節(jié)點，當(dāng)前任務(wù)完成時，它的后置任務(wù)進入就緒隊列，成為等待節(jié)點；關(guān)鍵節(jié)點處于層數(shù)和深度之和最大的關(guān)鍵路徑上，任務(wù)總體執(zhí)行時間由關(guān)鍵節(jié)點決定，關(guān)鍵節(jié)點在所有等待節(jié)點中具有最高優(yōu)先級［17］。

基于關(guān)鍵路徑的DAG動態(tài)調(diào)度算法的基本思想為：將河道網(wǎng)格作為初始關(guān)鍵路徑，當(dāng)就緒節(jié)點進入運行隊列執(zhí)行完成后，從運行隊列中自行剔除，并按照動態(tài)優(yōu)先級劃分算法更新執(zhí)行結(jié)束的計算節(jié)點的后置任務(wù)優(yōu)先級，從而在后續(xù)計算過程中不斷更新關(guān)鍵路徑并將計算節(jié)點按照任務(wù)優(yōu)先級分配運行。GXM并行計算過程中的任務(wù)調(diào)度如圖8所示。

圖8　GXM并行計算過程中的任務(wù)調(diào)度

基于關(guān)鍵路徑的DAG調(diào)度算法描述如下：

輸入 計算節(jié)點個數(shù)，線程數(shù)，就緒隊列中葉子節(jié)點個數(shù)，就緒隊列中關(guān)鍵節(jié)點個數(shù)，就緒隊列中等待節(jié)點個數(shù)，任務(wù)優(yōu)先級，運行隊列任務(wù)數(shù)，進入運行隊列閾值；

輸出 每個網(wǎng)格節(jié)點任務(wù)的計算結(jié)果。

根據(jù)流向網(wǎng)格構(gòu)建計算任務(wù)DAG

for←0 todo

初始化任務(wù)優(yōu)先級；

forall （in （1）） do

分配葉子節(jié)點進入就緒隊列；

for←0 todo

while（<=） do

if［］=1 &&0

等待節(jié)點進入運行隊列，計算流量，后置任務(wù)成為等待節(jié)點；

-1；

else if［］=1 &&≠0

關(guān)鍵節(jié)點進入運行隊列，計算流量，更新后置任務(wù)；

-1；

執(zhí)行計算任務(wù)，結(jié)果進入已完成隊列；

if=0

返回計算結(jié)果；

break

return 每個網(wǎng)格節(jié)點任務(wù)的計算結(jié)果

4 實驗與結(jié)果分析

4.1　實驗數(shù)據(jù)準(zhǔn)備

實驗硬件環(huán)境是一臺內(nèi)存大小16 GB，CPU型號為Intel E5-2630，操作系統(tǒng)為Windows Server 2008的個人計算機。實驗選取陜西省榆林市大理河2020年6月1日后30 d、60 d、90 d的流域水文數(shù)據(jù)與安徽省屯溪2020年6月1日后30 d的流域水文數(shù)據(jù)作為模型輸入，對本文提出的基于網(wǎng)格流向劃分的GXM并行算法與基于動態(tài)優(yōu)先級的DAG調(diào)度算法進行實驗驗證。

4.2　實驗結(jié)果分析

首先采用流域網(wǎng)格分辨率為1 km的大理河數(shù)據(jù)，設(shè)置預(yù)熱期為30 d內(nèi)的不同天數(shù)，對比串行算法和基于網(wǎng)格流向劃分的GXM并行算法在大理河流域內(nèi)的模型計算時間，結(jié)果如圖9（a）所示，可以看出在使用并行算法后：在預(yù)熱期為3 d時，GXM加速比為1.25；當(dāng)預(yù)熱期增加到30 d，GXM并行加速比達到4.03；并且模型預(yù)熱期達到21 d以后，串行算法的計算時間明顯增加，而GXM并行算法的計算時間一直沒有超過80 s，與原有的串行算法相比，基于網(wǎng)格流向劃分的GXM并行算法的計算效率得到了明顯提升。

采用基于共享內(nèi)存模型的OpenMP實現(xiàn)GXM并行算法，在大理河流域進行對比實驗，結(jié)果發(fā)現(xiàn)GXM并行加速比在預(yù)熱期為30天的情況下達到了5.01，模型計算時間進一步減少，模型計算效率得到了提高。

隨后采用流域網(wǎng)格分辨率為1 km的屯溪數(shù)據(jù)，設(shè)置預(yù)熱期為30 d內(nèi)的不同天數(shù)，對比串行算法和基于網(wǎng)格流向劃分的GXM并行算法在屯溪流域內(nèi)的模型計算時間，結(jié)果如圖9（b）所示。在預(yù)熱期為3 d時，GXM加速比為1.1；當(dāng)預(yù)熱期增加到30 d，GXM并行加速比達到4.11。可以看出，基于網(wǎng)格流向劃分的GXM并行算法的計算效率比原有的串行算法得到了明顯提高。

Freitas等［18］針對MGB（Modelo de Grandes Bacias）水文模型提出了一種基于矢量化與線程并行的CPU并行算法。該算法利用AVX-512矢量指令集與基于共享內(nèi)存模型的OpenMP對MGB中的STE（timestep）、DIS（discharge）、CON（continuity）三個主要程序模塊進行并行化處理。在尼日爾河流域?qū)υ撍惴ㄟM行實驗驗證，實驗發(fā)現(xiàn)基于矢量化與線程并行的CPU并行算法：在CPU型號為Intel Core i7-7820X的硬件環(huán)境下，設(shè)置計算線程數(shù)為8，MGB模型的最大加速比達到了5.27；在CPU型號為Intel Core i9-7900X的硬件環(huán)境下，設(shè)置計算線程數(shù)為10，MGB模型的最大加速比達到了5.99。與本文提出的基于網(wǎng)格流向劃分的GXM并行算法相比，基于矢量化與線程并行的CPU并行算法同樣擁有較好的并行性能，這是因為該算法不僅利用OpenMP實現(xiàn)了線程并行化，還利用矢量指令集實現(xiàn)了數(shù)據(jù)并行化。但是基于矢量化與線程并行的CPU并行算法仍然存在以下局限性：首先，該并行算法對計算機的硬件環(huán)境要求較高，因為該并行算法需要在支持AVX512矢量指令集的硬件環(huán)境下運行，然而目前很多計算機的CPU并不支持AVX512矢量指令集；其次，該并行算法是通過對模型內(nèi)部中的矢量操作進行矢量化改造來實現(xiàn)模型的并行化，這就要求模型內(nèi)部中無法進行矢量化改造的線性操作不能過多，否則無法達到算法預(yù)期的并行效果；并且該并行算法沒有指出如何減少并行計算過程中線程間的通信開銷。而本文提出的基于網(wǎng)格流向劃分的并行算法是從模型自身結(jié)構(gòu)出發(fā)，利用流域下墊面的DEM數(shù)據(jù)對流域網(wǎng)格進行區(qū)域分解，獲得了流域網(wǎng)格的演算次序，從而實現(xiàn)了模型的并行計算，因此本文提出的并行算法對計算機的硬件環(huán)境沒有嚴(yán)格要求；而且本文提出的基于動態(tài)優(yōu)先級的DAG調(diào)度算法還實現(xiàn)了計算任務(wù)優(yōu)先級的動態(tài)劃分，從而實現(xiàn)了并行計算過程中的任務(wù)調(diào)度，能有效降低線程間的通信開銷。

接著采用流域網(wǎng)格分辨率為0.5 km的大理河流域數(shù)據(jù)，相較于分辨率為1 km網(wǎng)格數(shù)據(jù)，模型的計算任務(wù)節(jié)點增加了1倍以上，模擬單元多達11 448個。設(shè)置預(yù)熱期為30 d、60 d、90 d情況下，采用基于動態(tài)優(yōu)先級的DAG調(diào)度算法對GXM并行計算進行任務(wù)調(diào)度，驗證調(diào)度算法對于GXM并行效率的影響，實驗結(jié)果如表1所示。

表1　調(diào)度算法對并行計算的影響

從表1可以看出，當(dāng)預(yù)熱期達到30 d時，GXM并行計算時間縮短了11.31%；當(dāng)預(yù)熱期達到60 d時，GXM并行計算時間縮短了30.08%；當(dāng)預(yù)熱期增加到90 d時，GXM并行計算時間相較于未采用調(diào)度算法的計算時間縮短了35.63%。因此，采用基于動態(tài)優(yōu)先級的DAG調(diào)度算法對GXM并行計算進行任務(wù)調(diào)度能夠有效提升GXM并行計算的效率。

針對90 d的預(yù)熱期，設(shè)置不同的計算線程數(shù)，監(jiān)測GXM并行計算的并行效率［19］，即加速比與參與計算的線程數(shù)比值，結(jié)果如表2所示。

由表2可以看出：采用了調(diào)度算法后的并行效率明顯高于未采用調(diào)度算法的并行效率，在線程數(shù)為4時，采用調(diào)度算法后的并行效率達到了0.90。隨著計算線程數(shù)的增加，GXM的并行效率在逐漸降低，這是因為當(dāng)計算線程數(shù)增加時，用于任務(wù)調(diào)度的計算資源開銷在逐漸增大，制約了并行計算性能的提升；但采用調(diào)度算法后的并行效率一直優(yōu)于未采用調(diào)度算法的并行效率。這表明通過基于動態(tài)優(yōu)先級的DAG調(diào)度算法能夠有效緩解當(dāng)計算線程數(shù)增加時出現(xiàn)的線程間負(fù)載不均衡的問題，減少計算線程的閑置等待時間，從而能較好地提升計算機的資源利用率。

表2　并行效率對比

5 結(jié)語

針對GXM的串行計算效率低的問題，本文提出了基于網(wǎng)格流向劃分的GXM并行算法。首先構(gòu)建網(wǎng)格流向矩陣，利用多層級可并行網(wǎng)格劃分算法得出流域網(wǎng)格的并行演算次序序列，實現(xiàn)GXM的并行計算，提升了GXM的計算速度；其次，針對并行計算中的任務(wù)調(diào)度問題，提出基于動態(tài)優(yōu)先級的DAG調(diào)度算法，通過并行演算次序序列構(gòu)建流域計算節(jié)點DAG并使用動態(tài)優(yōu)先級劃分算法動態(tài)更新網(wǎng)格計算節(jié)點優(yōu)先級，實現(xiàn)GXM并行計算中的任務(wù)調(diào)度，有效提升了GXM并行計算的資源利用率。未來會考慮如何在并行計算中減少網(wǎng)格計算單元間的大量數(shù)據(jù)交互并利用矢量指令實現(xiàn)GXM的數(shù)據(jù)并行化，進一步提升GXM的并行計算效率。

[1] 芮孝芳. 對流域水文模型的再認(rèn)識［J］. 水利水電科技進展， 2018， 38（2）： 1-7.（RUI X F. More discussion of watershed hydrological model［J］. Advances in Science and Technology of Water Resources， 2018， 38（2）： 1-7.）

[2] 芮孝芳. 論流域水文模型［J］. 水利水電科技進展， 2017， 37（4）：1-7， 58.（RUI X F. Discussion of watershed hydrological model［J］. Advances in Science and Technology of Water Resources， 2017， 37（4）：1-7， 58.）

[3] BEVEN K J， KIRKBY M J， FREER J E， et al. A history of TOPMODEL［J］. Hydrology and Earth System Sciences， 2021， 25（2）： 527-549.

[4] LIU Z， MARTINA M L V， TODINI E. Flood forecasting using a fully distributed model： application of the TOPKAPI model to the Upper Xixian Catchment［J］. Hydrology and Earth System Sciences， 2005， 9（4）： 347-364.

[5] 姚成. 基于柵格的分布式新安江模型構(gòu)建與分析［D］. 南京：河海大學(xué)， 2007： 16-20.（YAO C. Development and application of grid-based distributed Xin’anjinag model［D］. Nanjing： Hohai University， 2007： 16-20.）

[6] 姚成，李致家，張珂，等. 基于柵格型新安江模型的中小河流精細(xì)化洪水預(yù)報［J］. 河海大學(xué)學(xué)報（自然科學(xué)版）， 2021， 49（1）：19-25.（YAO C， LI Z J， ZHANG K， et al. Fine-scale flood forecasting for small and medium-sized rivers based on Grid-Xin’anjiang model［J］. Journal of Hohai University （Natural Sciences）， 2021， 49（1）：19-25.）

[7] 葉翔宇，李強，郭禹含，等. 高性能并行分布式水文模型研究進展［J］. 地理科學(xué)進展， 2022， 41（4）：731-740.（YE X Y， LI Q， GUO Y H， et al. Progress of research on high-performance parallel distributed hydrological model［J］. Progress in Geography， 2022， 41（4）： 731-740.）

[8] LIU J， ZHU A X， QIN C Z. Estimation of theoretical maximum speedup ratio for parallel computing of grid-based distributed hydrological models［J］. Computers and Geosciences， 2013， 60： 58-62.

[9] LIU J， ZHU A X， LIU Y， et al. A layered approach to parallel computing for spatially distributed hydrological modeling［J］. Environmental Modelling and Software， 2014， 51： 221-227.

[10] 秦澤寧，黎曙，周祖昊，等. 分布式水文模型區(qū)域分解并行計算方法及其應(yīng)用［J］. 水電能源科學(xué)， 2020， 38（10）：1-4， 12.（QIN Z N， LI S， ZHOU Z H， et al. Domain decomposition parallel computing method of distributed hydrological model and its application［J］. Water Resources and Power， 2020， 38（10）：1-4， 12.）

[11] 河海大學(xué). 基于流向的分布式水文模型并行計算方法： 202210254598.1［P］. 2022-05-10.（Hohai University. Parallel computing method for distributed hydrological models based on flow direction： 202210254598.1［P］. 2022-05-10.）

[12] LIU J， ZHU A X， QIN C Z， et al. A two-level parallelization method for distributed hydrological models［J］. Environmental Modelling and Software， 2016， 80： 175-184.

[13] 王想紅，劉紀(jì)平，徐勝華，等. 基于NetCDF數(shù)據(jù)模型的海洋環(huán)境數(shù)據(jù)三維可視化研究［J］. 測繪科學(xué)， 2013， 38（2）： 59-61.（WANG X H， LIU J P， XU S H， et al. Visualization of marine environment data based on NetCDF data model［J］. Science of Surveying and Mapping， 2013， 38（2）： 59-61.）

[14] 劉軍志，朱阿興，秦承志，等. 分布式水文模型的并行計算研究進展［J］. 地理科學(xué)進展， 2013， 32（4）： 538-547. （LIU J Z， ZHU A X， QIN C Z， et al. Review on parallel computing of distributed hydrological models［J］. Progress in Geography， 2013， 32（4）： 538-547.）

[15] 鄔倫，汪大明，張毅.基于DEM的水流方向算法研究［J］. 中國圖象圖形學(xué)報， 2006， 11（7）： 998-1003. （WU L， WANG D M， ZHANG Y. Research on the algorithms of the flow direction determination in ditches extraction based on grid DEM［J］. Journal of Image and Graphics， 2006， 11（7）： 998-1003.）

[16] 李致家，姚成，汪中華. 基于柵格的新安江模型的構(gòu)建和應(yīng)用［J］. 河海大學(xué)學(xué)報（自然科學(xué)版）， 2007， 35（2）： 131-134.（LI Z J， YAO C， WANG Z H. Development and application of grid-based Xin’anjiang model［J］. Journal of Hohai University （Natural Sciences）， 2007， 35（2）： 131-134.）

[17] 王命全，于炯，田園，等.網(wǎng)格環(huán)境中基于負(fù)載均衡的工作流調(diào)度算法［J］.計算機應(yīng)用，2010，30（12）：3184-3186.（WANG M Q， YU J， TIAN Y， et al. Workflow scheduling algorithm based on load balance in grid［J］. Journal of Computer Applications， 2010， 30（12）：3184-3186.）

[18] FREITAS H R A， MENDES C L， ILIC A. Performance optimization and scalability analysis of the MGB hydrological model［C］// Proceedings of the IEEE 27th International Conference on High Performance Computing， Data， and Analytics. Piscataway： IEEE， 2020： 31-40

[19] ZHANG A， LI T， SI Y， et al. Double-layer parallelization for hydrological model calibration on HPC systems［J］. Journal of Hydrology， 2016， 535： 737-747.

Parallel computing algorithm of grid-based distributed Xin’anjiang hydrological model

LIU Qian1， ZHANG Yangming1，2， WAN Dingsheng1*

（1，，211100，；2，210019，）

In recent years， the Grid-based distributed Xin’anjiang hydrological Model （GXM） has played an important role in flood forecasting， but when simulating the flooding process， due to the vast amount of data and calculation of the model， the computing time of GXM increases exponentially with the increase of the model warm-up period， which seriously affects the computational efficiency of GXM. Therefore， a parallel computing algorithm of GXM based on grid flow direction division and dynamic priority Directed Acyclic Graph （DAG） scheduling was proposed. Firstly， the model parameters， model components， and model calculation process were analyzed. Secondly， a parallel algorithm of GXM based on grid flow direction division was proposed from the perspective of spatial parallelism to improve the computational efficiency of the model. Finally， a DAG task scheduling algorithm based on dynamic priority was proposed to reduce the occurrence of data skew in model calculation by constructing the DAG of grid computing nodes and dynamically updating the priorities of computing nodes to achieve task scheduling during GXM computation. Experimental results on Dali River basin of Shaanxi Province and Tunxi basin of Anhui Province show that compared with the traditional serial computing method， the maximum speedup ratio of the proposed algorithm reaches 4.03 and 4.11， respectively， the computing speed and resource utilization of GXM were effectively improved when the warm-up period is 30 days and the data resolution is 1 km.

Grid-based distributed Xin’anjiang hydrological Model (GXM); grid flow direction division; parallel computing; Directed Acyclic Graph (DAG); task scheduling

1001-9081（2023）11-3327-07

10.11772/j.issn.1001-9081.2022111760

2022?11?24；

2023?02?15；

國家重點研發(fā)計劃項目（2018YFC1508106）。

劉乾（1998—），男，江蘇南京人，碩士研究生，CCF會員，主要研究方向：分布式水文模型并行計算； 張洋銘（1997—），男，江蘇徐州人，碩士研究生，CCF會員，主要研究方向：分布式水文模型并行計算； 萬定生（1963—），男，江蘇溧陽人，教授，CCF會員，主要研究方向：數(shù)據(jù)管理、數(shù)據(jù)挖掘。

P333

A

2023?02?17。

This work is partially supported by National Key Research and Development Program of China （2018YFC1508106）.

LIU Qian， born in 1998， M. S. candidate. His research interests include parallel computing of distributed hydrological models.

ZHANG Yangming， born in 1997， M. S. candidate. His research interests include parallel computing of distributed hydrological models.

WAN Dingsheng， born in 1963， professor. His research interests include data management， data mining.