CASC2D 分布式水文模型异构并行算法研究

左 翔 ，赵杏杏 ，丛小飞 ，刘修恒

（1. 南京河海智慧水利研究院，江苏 南京 210012；2. 南京中禹智慧水利研究院有限公司，江苏 南京 210012；3. 河海大学计算机与信息学院，江苏 南京 211100）

0 引言

分布式水文模型因能较好地描述水循环的时空变化过程而得到广泛的应用[1-5]。采用超渗产流机制的分布式水文模型 CASC2D 以栅格为空间离散单元，沿流向逐栅格进行汇流演算，因而可以对半干旱半湿润地区的水文过程进行较为准确的时空描述[6-8]。CASC2D 模型在现阶段的应用主要受计算效率限制，而计算效率主要取决于栅格单元的精度间隔[9]65，尺度划分越小，代表的地理特征越精细，模拟结果越准确，因此导致 CASC2D 模型计算量和计算时间呈几何倍数激增。提高水文模型计算效率，满足防洪“四预”实时预报的要求是亟待解决的关键问题[10]，采用并行计算技术是加速计算的有效手段。目前图形处理器GPU 计算能力提升明显，可以提供远超主流 CPU 的多线程并行计算能力[11-12]，科研人员将 GPU 逐步应用于水文预报领域，已取得一系列的成果[13-15]。

CUDA 是异构并行计算架构，充分利用 CPU 和GPU 各自的结构和性能优势，CPU 负责复杂的串行逻辑计算任务，GPU 负责处理数据密集型计算任务，两者协同提高计算效率[16]。在并行计算环境下，原CASC2D 模型程序涉及的数据流转、求解策略和算法结构需要重新考量，进行必要的结构调整和计算资源分配，以提升并行加速效果，目前采用异构并行编程技术改进 CASC2D 模型的研究鲜有报道。

本研究基于异构并行计算技术，构建了 CPU+GPU协同加速的 CASC2D 模型。原 CASC2D 模型的单指令多数据流计算转化为 CUDA 的内核程序（kernel）在 GPU 上的并发执行，逐栅格的串行计算转换为基于线程网格（Grid）、线程块（Block）的平行线程调度计算，栅格产流、坡面汇流和河道汇流等计算步骤采用不同的 kernel 程序代替，GPU 上执行的并行计算结果从显存传递回内存经 CPU 汇总后输出，实现CASC2D 模型的高效计算。将优化后的 CASC2D 模型应用于前毛庄子流域，模拟结果与原串行计算的模拟数据、流域的实测数据进行对比分析，以评估并行改进后 CASC2D 模型的模拟精度和计算效率。

1 CASC2D 模型

CASC2D 模型是基于物理机制和连续栅格单元演算的分布式模型[17]，以 DEM 为基础提取水系和划分计算栅格单元，采用基于距离平方倒数的降雨空间插值方法估算整个流域的降雨强度分布。降雨扣除植被截留等初始损失后采用 Green-Ampt 方程计算每个栅格单元的下渗率和累积下渗量。模型采用霍顿产流理论，当降雨强度大于下渗率时，栅格单元上产生的水流向较低的栅格流动，在水平面上分为 2 个流向。描述坡面汇流的控制方程采用基于圣维南的连续方程和动量方程，运用二维显式有限差分扩散波法描述坡面汇流[9]63，水流沿坡面汇集后，经河网汇流到达流域出口断面。河道汇流计算采用一维显式有限差分扩散波法，水流控制方程为一维明渠水流的连续方程，流量计算采用曼宁公式。模型演算流程如图1所示。

图1 分布式水文模型 CASC2D 演算流程

2 并行计算优化

分布式水文模型通常按地形或栅格划分子单元，各模拟单元按照时间序列进行迭代计算，具有空间可并行性强，时间可并行性弱的特点[18]，因此分布式水文模型主要采用空间离散化并行算法提升计算速度。刘军志等[19]通过分析各种水文模型的子过程，以及模拟单元的空间依赖关系，将水文过程的模拟方法分为计算独立型、邻域输入型、顺序依赖型和紧密耦合型4 类，其中计算独立型的水文子过程模拟单元的计算完全独立于其他单元，是最适合采用并行计算的类型，而紧密耦合型的水文子过程，是联合各个模拟单元进行的，并行计算难度较大。因此有必要对 CASC2D 模型的各个子过程进行分析，确定模拟方法类别，并评价并行优化的可行性。最终并行异构加速的模型程序仅引入多线程并行执行高密度计算部分，计算步骤需要与原 CASC2D 模型的步骤保持一致。

2.1 可并行性分析

为研究合适的并行计算方法，对图 1 中 CASC2D模型演算子过程进行逐一分析：

1） 数据读写。读取和写入栅格单元数据，计算迭代的步长计数，各栅格单元的分时变量，如下渗量、水位和流量数据需要传回主机端的内存，上述过程均通过 CPU 进行，效率主要受数据量大小和通信机制影响，无需进行并行改进。

2） 降雨计算。采用距离平方倒数法计算栅格单元面积上的降雨量。栅格雨量与邻近雨量站的实测雨强、雨量站栅格之间的距离有关，栅格单元之间不进行数据交换，具备并行优化可行性。

3） 产流计算。降雨在扣除植物截留等各种初始损失后开始下渗，栅格单元的植物截留与土地利用类型相关，模型计算每一个栅格单元的下渗率，涉及的计算参数包括饱和土壤水力传导度、毛管水头和土壤缺水量等，均由栅格的土壤类型确定。当栅格单元上的降雨强度大于下渗率时，剩余水量将会流向地势更低的栅格单元。同一时间步长内，栅格单元之间不进行水量交换，具备并行优化可行性。

4） 坡面汇流计算。坡面汇流的描述采用二维显式有限差分扩散波法，坡面汇流过程被分解为相等间隔的时间步长，以栅格单元为基本单位进行模拟计算，上一个时间步长的水深和水速作为下一个时间步长的初始水深和水速，各栅格单元的水量主要来源为邻域栅格单元的来水和降水补给。显式有限差分扩散波法要求各栅格单元仅计算当前时间步长内的输入、输出水量，不能影响当前时间步长内周围栅格单元的水量，以保证栅格单元的计算独立性。模型原计算程序的坡面汇流是逐栅格循环计算的，具备并行优化可行性。

5） 河道汇流。河道汇流采用一维显式有限差分扩散波法，与坡面汇流计算类似，原计算程序为逐河道节点更新，但同一时间步长内各单元计算独立，不影响周围单元，具备并行优化可行性。

通过对模型子过程的分析，确定可并行改进的主要内容，进一步分析得到影响 CASC2D 模型计算效率的因素有以下几点：

1） 数据量规模。基于 CUDA 的异构并行计算要在主机端（CPU）和设备端（GPU）进行数据传输，数据量的大小会在一定程度上影响计算效率，影响数据量的因素主要包括栅格数量和数据精度等。

2） 时间步长大小。通过增大模型计算的时间步长可以减少计算量，缩短模拟时间，但 CASC2D 模型采用的显式有限差分法，受计算稳定性限制，对时间步长要求非常严格，小步长模拟结果较优。

3） 计算复杂度。模型计算过程的复杂度越高，单次计算的时间就越长，并行计算后效率提升就越明显。

目前计算条件下单台计算机单次洪水 CASC2D模型模拟常需数小时以上。随着信息化水平的不断提高，可获取的流域信息不断丰富，模型计算的数据量不断提高，为保障计算效率，需要从硬件能力和模型算法 2 个角度配合，采用合适的并行计算机制实现CASC2D模型的优化加速，提升模型的实际应用效果。

2.2 并行加速方法

CASC2D 模型的并行加速需要结合 CPU+GPU的异构并行计算硬件架构对原程序重新设计，硬件架构如图2 所示。通过对原模型计算子过程的重新评估，根据计算特征分配给 CPU 或 GPU 执行，原程序中适合 CPU 计算或者必须由 CPU 顺序执行的子过程仍保留为 CPU 执行。

图2 异构并行计算硬件架构

并行加速后的 CASC2D 模型程序由串行和kernel程序 2 个部分组成，其中，串行程序由 CPU 运行计算，而 kernel 程序被 CPU 启动后交由 GPU 运行计算。模型执行某个并行计算任务的过程如图3 所示。

图3 CASC2D 模型并行计算执行流程

主要步骤如下：

1） 数据准备。文件读取、参数读取、动态内存分配、变量初始化由 CPU 依照原程序执行，将输入数据写入主机端内存，包括每个栅格单元的数字高程、土壤类型、降雨、河道断面等数据。

2） 分配显存空间。调用 CudaMalloc 函数在显存上为数字高程、土壤类型、土地利用等数据分配相应的区域。

3） 内存数据拷贝至显存。调用CudaMemcpyHost ToDevice 函数将步骤 2） 所示数据从内存拷贝至显存。

4） 启动并行计算。CPU 启动 kernel 程序，将栅格单元上的降雨、产流和汇流计算任务，按照合理的Grid和Block 层次划分，分别派发到 GPU 各个线程上进行并行计算。线程按行列坐标分配，对于栅格单元内的计算只需按行列号交由对应的线程执行，每个线程只负责执行与之对应的栅格单元计算，例如降雨并行计算的主要步骤：a. 获取当前线程的行列号，判断该线程是否在栅格坐标范围内；b. 计算当前线程所代表的栅格坐标；c. 计算当前栅格与不同雨量站之间的距离与雨量关系；d. 获取当前栅格处的降雨量。对于与单个栅格有关的子过程，如坡面和河道汇流等均可按上述类似方式完成计算任务。

5） 显存数据拷贝至内存。当需要进行结果输出时，调用CudaMemcpyDeviceToHost 函数，将 GPU 的计算数据从显存拷贝至内存，再进行结果输出。如果数据交换过于频繁，会导致计算效率下降，所以针对不同实例，按照需求调整输出时间步长尤为重要[20]。

6） 释放显存。调用 CudaFree 函数释放CudaMalloc函数分配的显存空间。

针对 CASC2D 模型模拟时间较长、应用性不高的特点，基于 CPU+GPU 的异构并行计算架构，利用几百甚至几千个流处理器同时调用 kernel 程序进行栅格单元的浮点计算，可以显著提升 CASC2D 模型的计算效率。

3 测试与分析

3.1 测试流域概况

为验证并行优化后的 CASC2D 模型计算效率，选取前毛庄流域为研究区域，流域水系如图4 所示。前毛庄流域位于天津于桥水库流域上游，流域面积为 435 km2，属于典型的闭合流域，降雨扣除蒸散发以外，下渗的水量小部分被植被吸收，其余大部分汇入河川径流，河川径流量及其时空分布特征主要受降雨影响。洪峰水量相对集中，洪水过程线呈陡涨陡落的趋势。土地利用简化为森林、耕地和草地 3 种类型，其中森林占比为 16.6%，耕地占比为 42.9%，草地占比为 40.5%。土壤类型中内潮褐土占比为 87%，本研究以潮褐土代表土壤类型。

图4 基于栅格单元的前毛庄流域水系图

3.2 测试条件

分别基于 CPU 串行计算和 CPU+GPU 异构并行计算进行测试，客观比较不同计算方法的计算效率，以加速比（串行计算执行时间与并行计算执行时间的比值）为主要评价指标。本研究测试机 CPU 为 Inter（R）Core（TM） i7-13700KF，显卡为 NVIDIA GeForce GTX 3080，核心数量为 8 704 个。数字高程资料源于分辨率为 30 m 的 DEM 数据，模拟时通过重采样将数据转换为不同精度进行对比测试。降雨资料的输入时段长为 1 h，且单位时段内雨强保持不变。采用线性内插的方法使实测流量资料的时段长与降雨保持一致，模型的输出时段长为输入时段长的倍数，为便于比较模拟与实测流量，模拟流量的输出时段长也取 1 h。

3.3 参数率定

CASC2D 模型参数率定按照先调整水量再调整流量，先调整峰值再调整峰现时差的原则，最后采用人工试错法进行参数率定。根据流域 DEM、土地利用类型、土壤类型及河道特征值等估算模型参数，再对一些敏感参数进行率定。前毛庄流域 CASC2D 模型参数值如表 1 和 2 所示。

表1 前毛庄流域 CASC2D 模型土地利用参数值

表2 前毛庄流域 CASC2D 模型土壤参数值

因 CASC2D 模型测试所需实测资料较多，不同参数设定会对模拟结果产生一定影响，故在不同测试条件下，模型参数应保持一致。

3.4 结果分析

模型参数率定结束后，选取前毛庄流域 2012 年的 20120727 号和 20120731 号 2 场洪水，对并行优化后的模型进行验证，综合评估影响模拟速度的因素和并行算法的加速性能。具体结果如下：

1） 计算步长的选择。计算步长是模型计算时的时间间隔，与栅格单元精度有关，对于采用显式有限差分法的 CASC2D 模型，当精度提高时，需要设置较短的计算步长，否则在一个时间步长内栅格单元容易出现水量疏干现象，导致计算不稳定和模拟偏差，计算量也因此成倍增大。针对测试的 30，90 和 300 m等 3 种 DEM 分辨率，需要选择合适的计算步长，一方面满足不同栅格单元数据精度下稳定计算的要求，另一方面方便比较在同等计算步长下并行加速的性能。选择精度最高、分辨率为 30 m 的 DEM，分别采用 1～6，8，10 s的计算步长进行试算，结果表明，仅在取 1～3 s 时水流计算稳定，在此范围内计算步长虽然能够在一定程度上影响洪峰流量的计算，但影响非常有限；取 4 s 时模拟流量出现不规则波动和放大发散现象，模拟效果较差；取 5，6，8 和 10 s 时程序会因栅格单元出现负值而终止运行。经过比较，最终决定采用 2 和 3 s 作为不同测试精度下的模拟计算步长，模拟结果如图5 所示。

图5 DEM 分辨率为 30 m 时不同计算步长的模拟结果

2） 数据精度的影响。在中小流域水文模型的实际应用中，DEM 分辨率数值常取 90 m，甚至 30 m 或更小，目的是更充分地反映流域产汇流特征的时空分布。以 20120731 号洪水为例，取计算步长为 2 s，DEM 分辨率分别为 30，90，200，300 和 400 m 进行径流模拟，模拟结果如图6 所示。在模型参数完全一致的情况下，随着 DEM 精度降低，模型模拟精度也会降低。洪峰洪量方面，随着 DEM 精度降低，模拟得到的洪峰流量值的总趋势是减小的，洪峰流量相对误差增大。峰现时间方面，不同精度模拟下的峰现时间与实测时间相差不大。确定性系数方面，总体趋势随着 DEM 精度的降低而减小：当 DEM 分辨率为 30 m 时，确定性系数达到0.90 以上；当 DEM分辨率为 200 m 时，确定性系数在0.85 左右；而当DEM 分辨率为 400 m 时，确定性系数低于 0.80。栅格单元数据精度越高，对地形概化越精细，虚假洼地蓄水效应越小，故单元蓄滞水量越小。栅格单元数据精度对洪水预报精度的影响较大，但当数据精度提高，计算时间也因此成倍增大，DEM 分辨率小于300 m 时，串行程序执行时间达到数小时以上。具体结果如表 3 所示。

图6 计算步长为 2 s 时不同 DEM 分辨率的模拟结果

表3 CASC2D 模型并行加速结果统计表

3） 模拟效率的提高。为评价并行加速效果，测试了 20120727 号和 20120731 号 2 场洪水 CASC2D模型串行和并行程序的计算耗时，统计结果见表 3。可以看出，CPU+GPU 异构并行计算技术可以大幅提高原程序的计算效率。参与计算的栅格单元越多，加速比越大，加速效率越高。在同样的数据精度条件下，计算步长对加速比影响不大。在 DEM 分辨率为200 m，计算步长为 3 s 的测试条件下，能够保证流量数值的稳定计算，洪水模拟耗时分别为 14 和 16 min，基本可以满足实时洪水预报预警的需求。因此本研究提出的异构并行计算方法具有显著的加速性能，为CASC2D 模型的应用提供了支撑。

4） 模拟精度的保持。为评估并行程序改造对CASC2D 模型模拟精度的影响，与原串行程序的模拟结果进行比较，对比结果如图7 所示，计算所得的流量过程与原串行程序完全一致，说明 CPU+GPU 异构并行技术不仅可以大幅提高原程序的运行效率，模拟精度也与原程序保持一致。

图7 原 CPU 串行计算与 CPU+CPU 异构并行计算的模拟流量对比图

4 结语

基于栅格单元的 CASC2D 模型，由于采用二维显式有限差分扩散波法描述坡面汇流过程，在实际水文模拟过程中存在计算量大、模拟时间较长等问题。采用 CPU+GPU 异构并行计算架构对原 CASC2D 模型程序进行了重构，为子过程计算分配合理的 GPU资源，测试结果表明：栅格单元数据精度越高，对地形概化越精细，模拟精度越高；栅格单元数量越多，GPU 并行计算的加速性能越好。当 DEM 分辨率为200 m，计算步长为 2 s 时，模拟前毛庄流域 20120727 号和 20120731 号洪水，计算耗时约为 20 min，基本满足中小河流洪水实时预警预报需求。并行优化后的模型得到了与原模型等效的模拟效果，显著提高了CASC2D 的应用潜力和科研价值。本研究中 CASC2D模型水动力过程采用的是显式有限差分法求解，为进一步提高模拟精度，未来可以采用隐式有限差分法求解，并探索其并行加速的方法。