基于GPU加速的水文模型參數率定

闞 光 遠，洪 陽，梁 珂，何 曉 燕，丁 留 謙，張 大 偉

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

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

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

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

1 研究區域與數據

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

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

用于模型參數率定的水文氣象資料為日尺度，對應的日降雨和日蒸散發能力為2000～2010年共計11 a數據。為避免實測徑流資料誤差導致的參數率定結果不穩定，使用實測日降雨和蒸散發能力和一組人工設定的理想參數(這里采用新安江模型結構性約束KG+KI=0.7)驅動新安江模型，獲取一套人工理想徑流數據并將其作為參數率定的目標徑流過程。

人工設定的理想參數如下：

參數數理意義范圍理想值K蒸散發折算系數0.1～1.50.723B張力水容量分布曲線方次0.1～0.40.3C深層散發系數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地下水日出流系數0.1～0.60.4CG地下水日消退系數0.98～0.9980.994CI壤中水日消退系數0～0.90.8CS河網匯流日消退系數0～10.2L河網匯流日滯時0～5 d0 dXE馬法河道演算形狀系數-0.5～0.50.2

2 基于GPU加速的水文模型參數率定

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

(1) 硬件平臺搭建。基于4臺高性能HP工作站，搭建GPU并行計算Windows集群，集群中安裝的GPU見圖2。用于搭建硬件平臺的設備包括1臺HP Z820高性能工作站和3臺HP Z840高性能工作站，集群由TP-Link 8口千兆交換機通過4條山澤超五類網線實現高速互聯。為方便集群管理，使用邁拓維矩4口KVM DVI切換器實現多臺工作站共享一套鼠標、鍵盤和顯示器。本系統CPU核心總數為16 × 1(Xeon E5-2640v2) + 16 × 5(Xeon E5-2630v3) = 96個，GPU核心總數為2 880 × 6(Tesla K40c)=17 280個，峰值單精度浮點計算性能可達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位旗艦版操作系統。開發環境采用Microsoft Visual Studio 2010，采用VC++2010作為主力編程語言。采用NVIDIA CUDA C/C++ 6.5作為GPU端代碼的編譯器。

2.2 并行SCE-UA算法及其程序化實現

并行SCE-UA算法設計與GPU硬件架構的對應關系見圖3。并行SCE-UA算法由圖中多個并行子功能模塊組合而成，各子模塊介紹如下。

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

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

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

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

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

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

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

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

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

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

3 研究結果

3.1 率定結果一致性和穩定性分析

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

數值實驗結果表明，50次運行后，串行和并行SCE-UA算法均成功收斂至理想參數值。參數率定結果準確且穩定。串行和并行SCE-UA算法率定的模型表現一致，模擬精度高，納什效率系數均接近1。

以上結果證明，對于馬渡王流域11 a的日徑流模擬，1 000萬次目標函數評價次數足夠保證算法穩定一致的收斂至理想參數值。因此，本研究的計算性能比較均基于1 000萬次目標函數評價次數，排除由于評價次數不足導致的算法收斂不穩定影響。

3.2 運行時間比較

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

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

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

3.3 加速比比較

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

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

3.4 資料長度影響分析

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

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

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

圖6 資料長度與運行時間的關系Fig.6 Relationship between data length and execution time

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

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

4 結 論

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

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

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

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