基于GPU加速的水文模型参数率定

2019-05-31阚光远，洪阳，梁珂，何晓燕，丁留谦，张大伟

人民长江 2019年5期

阚光远，洪阳，梁珂，何晓燕，丁留谦，张大伟

(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小节的结论相一致。