闞 光 遠(yuǎn)，洪 陽，梁 珂，何 曉 燕，丁 留 謙，張 大 偉

(1.清華大學(xué) 水利系，北京 100084； 2.中國水利水電科學(xué)研究院 北京中水科工程總公司，北京 100048；3.中國水利水電科學(xué)研究院 水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038)

近年來，水文模型被廣泛應(yīng)用于水文模擬與預(yù)報、水資源預(yù)測規(guī)劃與管理、氣候變化對水文過程影響分析等領(lǐng)域，發(fā)揮著重要作用。水文模型模擬結(jié)果的好壞很大程度上取決于模型參數(shù)的優(yōu)劣。水文模型中的過程參數(shù)必須進(jìn)行率定。由于水文過程的高度復(fù)雜性[1]，模型參數(shù)率定一直是一個前沿和熱點(diǎn)問題[2-3]。模型參數(shù)率定方法有人工試錯法和自動率定法兩類。人工試錯法不易找到最優(yōu)的模型參數(shù)，調(diào)試時間長，結(jié)果依賴于調(diào)試人員的經(jīng)驗(yàn)，增加了模型的不確定性[4]。隨著計(jì)算機(jī)技術(shù)的發(fā)展和數(shù)學(xué)方法的進(jìn)一步應(yīng)用，模型參數(shù)自動率定已成為主流。這類方法具有尋優(yōu)操作簡便、尋優(yōu)結(jié)果客觀等優(yōu)點(diǎn)[5]。模型參數(shù)自動率定方法很多，常見的有單純形法[6]、Rosenbrock法[7]、模式搜索法[8]以及SCE-UA算法[9-10]等。目前在水文模型參數(shù)率定領(lǐng)域應(yīng)用較廣的方法是SCE-UA算法。該算法是Duan等在亞利桑那州大學(xué)發(fā)展的一種復(fù)雜的基于非線性單純形法和遺傳算法的混合單目標(biāo)優(yōu)化方法。

新安江模型[11-12]在我國已成功應(yīng)用幾十年，其參數(shù)率定一直是一個熱、難點(diǎn)問題，很多學(xué)者基于SCE-UA算法對新安江模型參數(shù)自動率定進(jìn)行過研究，取得了許多有益的成果。但在復(fù)雜的實(shí)際應(yīng)用中，新安江模型參數(shù)自動率定仍存在收斂速度緩慢的問題。隨著國民經(jīng)濟(jì)的發(fā)展，國家對水文預(yù)報的精度和響應(yīng)速度有了更高的要求，水文站網(wǎng)(尤其是雨量站)的布設(shè)密度顯著增加，實(shí)測降雨、蒸散發(fā)和流量序列長度逐年增長，給水文預(yù)報精度和可靠性的提升帶來了希望[13-15]。然而，數(shù)據(jù)量的增大導(dǎo)致了計(jì)算量的劇增，引發(fā)了預(yù)報精度和響應(yīng)速度之間的矛盾，增加了災(zāi)害損失的潛在風(fēng)險。對于預(yù)報斷面以上密集布設(shè)的雨量站和幾十年的長系列歷史資料，新安江模型參數(shù)自動率定的計(jì)算速度非常緩慢，往往需要數(shù)天時間才能完成一次率定。隨著氣候變化和人類活動的加劇，許多流域的下墊面狀況發(fā)生了顯著變化，以往率定好的模型參數(shù)不再適用，參數(shù)的重新率定成為當(dāng)務(wù)之急，過慢的參數(shù)率定速度嚴(yán)重降低了工作效率。傳統(tǒng)的基于中央處理器(CPU)的多核并行計(jì)算加速技術(shù)由于受到核心數(shù)目和功耗的限制，在大數(shù)據(jù)領(lǐng)域的表現(xiàn)不能令人滿意，急需基于圖形處理器(GPU)的新一代眾核并行計(jì)算技術(shù)的支持[15-20]。

綜上所述，基于GPU加速的新安江模型參數(shù)率定研究具有科學(xué)意義和實(shí)用價值。本文深入研究SCE-UA算法和新安江模型的原理，透徹挖掘算法的并行性，并與新一代高性能GPU并行計(jì)算軟硬件平臺緊密結(jié)合，研發(fā)并行參數(shù)率定軟件系統(tǒng)。以研發(fā)的并行參數(shù)率定軟硬件平臺為有力工具，深入分析并行率定方法在GPU上的加速效果,從而解決新安江模型參數(shù)自動率定加速這一關(guān)鍵問題。

1 研究區(qū)域與數(shù)據(jù)

本研究基于陜西省馬渡王流域。馬渡王流域是典型半濕潤半干旱流域，流域出口馬渡王水文站位于西安市馬渡王村，集水面積1 601 km2，距河口距離30 km。馬渡王河段順直，河中靠右岸有一河灘，河床為砂卵石。該站測驗(yàn)項(xiàng)目包括水位、流量、泥沙、降水、蒸發(fā)和水質(zhì)等水文參量；洪水過程峰型瘦，陡漲陡落，水位流量關(guān)系散亂。流域分布有10個雨量站，根據(jù)泰森多邊形法劃分為10個子流域，見圖1。

圖1 馬渡王流域Fig.1 Maduwang watershed map

用于模型參數(shù)率定的水文氣象資料為日尺度，對應(yīng)的日降雨和日蒸散發(fā)能力為2000～2010年共計(jì)11 a數(shù)據(jù)。為避免實(shí)測徑流資料誤差導(dǎo)致的參數(shù)率定結(jié)果不穩(wěn)定，使用實(shí)測日降雨和蒸散發(fā)能力和一組人工設(shè)定的理想?yún)?shù)(這里采用新安江模型結(jié)構(gòu)性約束KG+KI=0.7)驅(qū)動新安江模型，獲取一套人工理想徑流數(shù)據(jù)并將其作為參數(shù)率定的目標(biāo)徑流過程。

人工設(shè)定的理想?yún)?shù)如下：

參數(shù)數(shù)理意義范圍理想值K蒸散發(fā)折算系數(shù)0.1～1.50.723B張力水容量分布曲線方次0.1～0.40.3C深層散發(fā)系數(shù)0.08～0.20.15WUM上層張力水容量5～30 mm20 mmWLM下層張力水容量50～100 mm60 mmWDM深層張力水容量15～70 mm40 mmIM不透水面積0.01～0.020.01SM自由水容量10～50 mm12EX自由水容量分布曲線方次1～21.5KG地下水日出流系數(shù)0.1～0.60.4CG地下水日消退系數(shù)0.98～0.9980.994CI壤中水日消退系數(shù)0～0.90.8CS河網(wǎng)匯流日消退系數(shù)0～10.2L河網(wǎng)匯流日滯時0～5 d0 dXE馬法河道演算形狀系數(shù)-0.5～0.50.2

2 基于GPU加速的水文模型參數(shù)率定

2.1 基于CUDA和GPU的并行計(jì)算軟硬件平臺搭建

(1) 硬件平臺搭建?；?臺高性能HP工作站，搭建GPU并行計(jì)算Windows集群，集群中安裝的GPU見圖2。用于搭建硬件平臺的設(shè)備包括1臺HP Z820高性能工作站和3臺HP Z840高性能工作站，集群由TP-Link 8口千兆交換機(jī)通過4條山澤超五類網(wǎng)線實(shí)現(xiàn)高速互聯(lián)。為方便集群管理，使用邁拓維矩4口KVM DVI切換器實(shí)現(xiàn)多臺工作站共享一套鼠標(biāo)、鍵盤和顯示器。本系統(tǒng)CPU核心總數(shù)為16 × 1(Xeon E5-2640v2) + 16 × 5(Xeon E5-2630v3) = 96個，GPU核心總數(shù)為2 880 × 6(Tesla K40c)=17 280個，峰值單精度浮點(diǎn)計(jì)算性能可達(dá)0.032 × 1 (Xeon E5-2640v2)+ 0.0384 × 5(Xeon E5-2630v3) + 4.29 × 6 (Tesla K40c)= 25.964TFLOPS(萬億次/s)。

圖2 Tesla K40c GPUFig.2 Tesla K40c GPU

(2) 軟件平臺搭建。本研究采用Microsoft Windows 7 64位旗艦版操作系統(tǒng)。開發(fā)環(huán)境采用Microsoft Visual Studio 2010，采用VC++2010作為主力編程語言。采用NVIDIA CUDA C/C++ 6.5作為GPU端代碼的編譯器。

2.2 并行SCE-UA算法及其程序化實(shí)現(xiàn)

并行SCE-UA算法設(shè)計(jì)與GPU硬件架構(gòu)的對應(yīng)關(guān)系見圖3。并行SCE-UA算法由圖中多個并行子功能模塊組合而成，各子模塊介紹如下。

(1) 算法并行化-并行初始點(diǎn)群生成。SCE-UA算法初次運(yùn)行時需要根據(jù)算法參數(shù)、決策變量的個數(shù)來生成相應(yīng)數(shù)量的初始點(diǎn)作為進(jìn)化計(jì)算的父代樣本，包含兩個可并行的步驟：隨機(jī)點(diǎn)的生成和各點(diǎn)目標(biāo)函數(shù)值的計(jì)算。利用GPU生成多條并行線程，每條線程控制一個樣本點(diǎn)。在每一條線程內(nèi)，利用梅森旋轉(zhuǎn)算法生成一個隨機(jī)點(diǎn)，將隨機(jī)點(diǎn)的坐標(biāo)值作為模型參數(shù)，調(diào)用新安江模型，進(jìn)行模型模擬計(jì)算，并基于實(shí)測和模型模擬的流量過程進(jìn)行目標(biāo)函數(shù)值(納什效率系數(shù))的計(jì)算。

(2) 算法并行化-并行復(fù)合形競爭進(jìn)化。利用GPU生成多條并行線程，每條線程控制一個復(fù)合形的進(jìn)化過程。為每條線程建立一個隨機(jī)數(shù)種子，每個復(fù)合形基于自己的隨機(jī)數(shù)種子和梅森旋轉(zhuǎn)隨機(jī)數(shù)發(fā)生器各自獨(dú)立地進(jìn)行復(fù)合形進(jìn)化。當(dāng)規(guī)定的進(jìn)化次數(shù)完成后，多條線程的進(jìn)化結(jié)果合并入一條主線程，進(jìn)行復(fù)合形混合操作，之后再次拆分為多條線程，繼續(xù)進(jìn)行并行復(fù)合形競爭進(jìn)化，如此往復(fù)，直到算法收斂。

(3) 算法并行化-并行適應(yīng)度堆排序。復(fù)合形進(jìn)化每完成一批競爭進(jìn)化，需要將各復(fù)合形生成的新點(diǎn)群進(jìn)行混合、排序和重排序(又稱為“洗牌”操作)。排序基于目標(biāo)函數(shù)值進(jìn)行，采用基于GPU加速的并行堆排序算法，能夠在不改變原點(diǎn)群指針變量的基礎(chǔ)上對點(diǎn)群索引進(jìn)行排序，方便下一步基于索引進(jìn)行的點(diǎn)群及其目標(biāo)函數(shù)值重排序。

(4) 算法并行化-并行點(diǎn)群重排序。復(fù)合型混合、排序完成后，由GPU生成多條并行線程，每條線程控制一個樣本點(diǎn)，基于并行堆排序獲取的索引，并行地對原點(diǎn)群指針變量中的點(diǎn)群及其目標(biāo)函數(shù)值進(jìn)行重排序。

(5) 算法并行化-并行梅森旋轉(zhuǎn)隨機(jī)數(shù)發(fā)生器。在算法執(zhí)行的全過程中，很多子程序需要用到隨機(jī)數(shù)，為了避免各并行線程同時改寫同一個隨機(jī)數(shù)種子帶來的數(shù)據(jù)沖突和競寫問題(data racing)，由GPU生成多條并行線程，每條線程控制一個隨機(jī)數(shù)種子，調(diào)用屬于本線程的梅森旋轉(zhuǎn)隨機(jī)數(shù)發(fā)生器，并行地實(shí)現(xiàn)隨機(jī)數(shù)的生成。

(6) 程序化實(shí)現(xiàn)-CPU內(nèi)存與GPU顯存雙向高效拷貝。對于計(jì)算過程中需要長期駐留GPU顯存的標(biāo)量、指針等，在程序啟動后統(tǒng)一進(jìn)行一次顯存分配(調(diào)用cudaMalloc函數(shù))，程序運(yùn)行結(jié)束后，統(tǒng)一進(jìn)行一次清理(調(diào)用cudaFree函數(shù))，有助于提升代碼的執(zhí)行效率。對于需要在主機(jī)端(即CPU端，下同)和設(shè)備端(即GPU端，下同)頻繁傳遞的中間變量，通過調(diào)用cudaMemcpy函數(shù)實(shí)現(xiàn)CPU內(nèi)存與GPU顯存雙向高效拷貝。

(7) 程序化實(shí)現(xiàn)-GPU存儲器層次設(shè)計(jì)與性能調(diào)優(yōu)。為了提高代碼執(zhí)行效率，GPU將存儲器劃分為多個層次。因此，根據(jù)算法中各變量的生存周期、作用域等，將它們配置到不同的GPU存儲器層次中，達(dá)到最優(yōu)性能。需要全局可見且指向的數(shù)值需要頻繁修改的指針，分配為全局變量(如復(fù)合形點(diǎn)群坐標(biāo)等)。需要全局可見但數(shù)值在整個程序生命周期內(nèi)不變的標(biāo)量，分配為常量顯存(如復(fù)合型個數(shù)等)。需要全局可見但指向的數(shù)值在整個程序生命周期內(nèi)不變的指針，分配為紋理顯存(如人工理想徑流量等)。在線程塊內(nèi)線程間共享的變量，分配為共享變量。對于僅線程自身可見的變量，分配為線程私有變量(如臨時變量)。

(8) 程序化實(shí)現(xiàn)-多線程同步及原子操作。并行SCE-UA算法運(yùn)算流程中有幾個同步點(diǎn)需要進(jìn)行處理。主要包括：初始點(diǎn)群生成后的同步、復(fù)合形競爭進(jìn)化后的同步、適應(yīng)度排序和點(diǎn)群重排序后的同步等全局同步點(diǎn)(通過啟動和關(guān)閉核函數(shù)實(shí)現(xiàn))。復(fù)合形進(jìn)化過程中的計(jì)數(shù)器變量需要進(jìn)行自增運(yùn)算，這一運(yùn)算也需要進(jìn)行同步，由于為四則運(yùn)算，因此采用輕量級高效率的原子操作(atomicAdd)實(shí)現(xiàn)。

2.3 并行SCE-UA算法在GPU上的加速效果

通過調(diào)整SCE-UA算法的算法參數(shù)(復(fù)合形個數(shù))和水文資料規(guī)模(資料系列長度)，測試并行SCE-UA算法的性能特征。復(fù)合形個數(shù)調(diào)整范圍設(shè)置為4～1 024，遞增倍數(shù)為2。資料系列長度從1 a開始，每次增加1，直至最大年數(shù)(11 a)為止。與基于單核CPU的串行版程序進(jìn)行計(jì)算結(jié)果一致性、運(yùn)行時間、加速比、資料長度影響的比較和分析。

3 研究結(jié)果

3.1 率定結(jié)果一致性和穩(wěn)定性分析

為保證對比的公平性，將參數(shù)率定的總計(jì)算量(目標(biāo)函數(shù)評價次數(shù))固定為1 000萬次。本小節(jié)采用的水文氣象資料長度均為11 a。通過給定不同的初始模型參數(shù)，反復(fù)(50次)運(yùn)行串行、并行SCE-UA算法進(jìn)行一致性和穩(wěn)定性比較與分析。

數(shù)值實(shí)驗(yàn)結(jié)果表明，50次運(yùn)行后，串行和并行SCE-UA算法均成功收斂至理想?yún)?shù)值。參數(shù)率定結(jié)果準(zhǔn)確且穩(wěn)定。串行和并行SCE-UA算法率定的模型表現(xiàn)一致，模擬精度高，納什效率系數(shù)均接近1。

以上結(jié)果證明，對于馬渡王流域11 a的日徑流模擬，1 000萬次目標(biāo)函數(shù)評價次數(shù)足夠保證算法穩(wěn)定一致的收斂至理想?yún)?shù)值。因此，本研究的計(jì)算性能比較均基于1 000萬次目標(biāo)函數(shù)評價次數(shù)，排除由于評價次數(shù)不足導(dǎo)致的算法收斂不穩(wěn)定影響。

3.2 運(yùn)行時間比較

圖4為串行、并行SCE-UA算法運(yùn)行時間圖。從圖中可見，串行算法的運(yùn)行時間幾乎不變，這是由于總的計(jì)算量已固定為1 000萬次目標(biāo)函數(shù)評價。

對于并行SCE-UA算法，復(fù)合形個數(shù)代表了并行算法的并行度。通常，隨著并行度的增大，運(yùn)行時間呈下降趨勢。由圖4可見，隨著復(fù)合形個數(shù)從4增大至1 024，運(yùn)行時間相應(yīng)從774 196.6 s下降至4 403.4 s。此外，當(dāng)并行度較小時(如4～16個復(fù)合形)，并行SCE-UA算法比串行版的耗時要高，并沒有取得加速效果。這是因?yàn)樵诘筒⑿卸认?，程序的性能主要由處理器核心頻率決定，而Tesla K40c GPU的核心頻率(875 MHz)遠(yuǎn)小于Xeon CPU的核心頻率(2.0 GHz)，故并行版的運(yùn)行速度反而更慢。當(dāng)并行度增大到一定程度(如32～1 024個復(fù)合形)，并行版的性能大大超越了串行版。這是因?yàn)镚PU可以創(chuàng)建大量線程，并行地完成任務(wù)，大量線程并行帶來的加速效果已經(jīng)超越了低處理器核心頻率的不良影響，因此并行版能夠取得更好的加速效果。

圖4 串行、并行SCE-UA算法運(yùn)行時間Fig.4 Execution time of serial and parallel SCE-UA algorithm

3.3 加速比比較

串行、并行SCE-UA算法加速比見圖5。由圖中可見，隨著復(fù)合形個數(shù)的增加，加速比逐漸增大。當(dāng)并行度較小時(復(fù)合形個數(shù)介于4～16之間)，并行SCE-UA算法并未取得加速效果。這是由于并行度較小時，程序的運(yùn)行速度主要受核心頻率影響，此時GPU較低的核心頻率影響較大，導(dǎo)致加速效果不理想。當(dāng)并行度足夠大時(復(fù)合形個數(shù)大于等于32時)，并行SCE-UA算法取得了加速效果，且并行度越大，加速比越高，最高可達(dá)41.39倍，此時復(fù)合形個數(shù)為1 024。這是由于大并行度時，GPU能夠創(chuàng)建大量的并行線程，程序性能的決定性因素轉(zhuǎn)化為高并行度而非核心頻率，因此并行版能取得理想的加速效果。以上分析結(jié)果與3.2小節(jié)的結(jié)論相一致。

圖5 串行、并行SCE-UA算法加速比Fig.5 Speed-up ratio of serial and parallel SCE-UA algorithm

3.4 資料長度影響分析

目標(biāo)函數(shù)計(jì)算負(fù)荷是影響SCE-UA算法性能的重要因素之一。本小節(jié)通過改變資料長度(從1 a遞增到11 a)實(shí)現(xiàn)目標(biāo)函數(shù)計(jì)算負(fù)荷的改變，進(jìn)而進(jìn)行加速效果影響分析。為排除低并行度帶來的干擾，本小節(jié)中，復(fù)合形個數(shù)統(tǒng)一設(shè)定為1 024。

圖6為資料長度與運(yùn)行時間關(guān)系圖。由圖中可見，隨著資料長度的增加，串行、并行SCE-UA算法的運(yùn)行時間都有所增加，但并行SCE-UA算法的運(yùn)行時間增長幅度明顯小于串行版程序。

圖7為資料長度與加速比關(guān)系圖。由圖中可見，隨著資料長度的變化，并行SCE-UA算法的加速比變化并不十分顯著，介于35～42倍之間。

圖6 資料長度與運(yùn)行時間的關(guān)系Fig.6 Relationship between data length and execution time

圖7 資料長度與加速比的關(guān)系Fig.7 Relationship between data length and speed-up ratio

綜上，隨著資料長度的增加，串行SCE-UA算法運(yùn)行時間顯著增長，并行SCE-UA算法運(yùn)行時間呈增長趨勢，但增長幅度遠(yuǎn)小于串行SCE-UA算法。加速比隨資料長度的增長變化并不十分顯著。

4 結(jié) 論

針對水文模型參數(shù)率定速度緩慢的難題，采用基于GPU加速的并行計(jì)算技術(shù)建立了水文模型參數(shù)快速率定方法，在基于CUDA和GPU的軟硬件平臺上進(jìn)行了數(shù)值實(shí)驗(yàn)，通過性能分析，得出以下結(jié)論。

(1) 基于GPU加速的并行SCE-UA算法能夠穩(wěn)定收斂得到理想?yún)?shù)值。在算法結(jié)果一致性和穩(wěn)定性上取得了很好的效果，在新安江模型參數(shù)率定的實(shí)際應(yīng)用中具有實(shí)用價值。

(2) 并行SCE-UA算法收斂速度快，尤其是在大并行度(本文中復(fù)合形個數(shù)大于等于32)場合下，運(yùn)行時間遠(yuǎn)小于串行SCE-UA算法，能夠取得非常理想的加速比。

(3) 資料長度的增長會導(dǎo)致串行和并行SCE-UA算法的運(yùn)行時間變長，但并行SCE-UA算法的運(yùn)行時間增長幅度遠(yuǎn)小于串行SCE-UA算法。此外，盡管運(yùn)行時間受資料長度影響顯著，但資料長度對加速比的影響卻相對較小。