地理计算的本质特点是数据密集,多模型多尺度时空动态耦合,计算模型十分复杂。溃坝洪水等地理现象作为一个复杂巨系统,具有多主体、多因素、多尺度、多变性的特征,其时空动态发展过程往往比最终形成的空间格局更为重要^{[1, 2]}。因此,实时高效地进行溃坝洪水时空发展过程模拟分析及预测预警十分关键。

元胞自动机(cellular automata,CA)是一种时间、空间、状态都离散,空间相互作用和时间因果关系都为局部的网格动力学模型,具有模拟复杂系统时空演化过程的能力^{[2, 3]},在地理模拟方面,已广泛应用于城市演变^{[4, 5]}、土地利用^{[6, 7]}、交通模拟^[8]、森林火灾^[9]等领域中。在溃坝洪水模型研究方面,基于CA构建的溃坝洪水演进模型具有并行度高,几何边界容易处理,求解简单等特点,已经成为溃坝洪水时空过程模拟分析的一个研究热点^{[10, 11]}。由于洪水溃坝事件具有发生的突然性、扩张的快速性以及响应的紧迫性等特性,影响因素复杂多变,因而对其发展过程和影响区域分析的直观性、快速性和准确性的要求更高。而现有的基于元胞自动机的溃坝洪水模型^{[12, 13, 14, 15]}采用CPU (central processing unit)串行的计算方式进行溃坝洪水演进过程模拟与分析,模型的计算效率偏低。尤其是在使用高分辨率DEM数据进行精细模拟与分析时,系统响应和交互效率低下,难以达到"实时"地进行溃坝洪水演进模拟分析和风险评估的要求。

近几年,随着GPU (graphics processing unit)的性能提高,计算机系统的可编程性和功能都得到极大的扩展,能够支持更复杂的运算,并出现了可以将GPU作为数据并行计算设备的通用并行计算架构(compute unified device architecture,CUDA)^[16]。目前,基于GPU的并行算法已经被广泛用于土地利用模拟^[17]、图像处理^{[18, 19]}、水文模型^[20]、地球物理^[21]、天文计算^[22]、地形绘制^[23]等领域,并在大多数应用中获得了几倍、几十倍乃至上百倍的加速比^[24]。本文针对元胞自动机具有同步计算(并行)的特点,将GPU高性能计算技术和CA模型进行结合,提出了GPU-CA的溃坝洪水演进计算模型。将元胞自动机的局部迭代计算过程映射到GPU中进行,通过大量线程的并行计算提高溃坝洪水模拟的计算效率,实现溃坝洪水演进"实时"计算与模拟分析,支持快速的风险评估。

元胞自动机是定义在一个由具有离散、有限状态元胞组成的空间上,并按照一定局部规则,在离散时间维上演化的动力学系统^{[2, 3]}。它由元胞、元胞空间、邻域及规则等4部分组成,可视为一个元胞空间和定义于该空间的变换函数所组成^{[2, 3]}。

本文采用二维空间分布的元胞自动机模拟框架对洪水演进的时空动态过程进行模拟,将实际研究区按照一定分辨率划分成离散格网,格网单元为正方形单元,以固定的时间单元间隔Δt不断演化。邻域为冯\5诺依曼邻元类型,将4个相邻的格网作为一个元胞空间,其状态通过X、Y方向单宽流量和水面高程信息进行描述。

溃坝洪水演进计算采用离散的圣维南方程组^[2],整个演进过程分为两步进行计算,首先由t时刻的水深计算t+1时刻的单宽流量:

再由t+1时刻的单宽流量计算t+2时刻的水深:

式中,M_i,j^t和N_i,j^t为元胞(i,j)在t时刻X和Y方向上的单宽流量(m²/s);h_i,j^t为元胞(i,j)在t时刻的水深(m);u_i,j^t和v_i,j^t分别表示为元胞(i,j)在t时刻在X、Y方向上的流速(m/s);n_i,j²为元胞(i,j)的糙率(m^-1/3s);Δt表示时空比;ΔX和Δy分别表示X和Y方向上的尺度大小(m)。

元胞自动机建模采用典型"自下而上"的构模方法,每个元胞通过不断的迭代计算得到溃坝洪水动态演进过程,在计算过程中,每个元胞只与邻近元胞进行交互运算。GPU适合处理后能够转化为数据并行计算的问题。本文充分利用GPU与CA的优势,顾及多种边界约束条件及参数修正,修改溃坝模型计算算法,设计了如图 1所示的基于GPU-CA的溃坝洪水演进计算模型,将大量在CPU上串行计算的方法转换为在GPU上并行计算,以提高模型计算效率。计算过程中,首先在CPU端即host端进行任务的初始化,包括内存和GPU显存分配、数据初始化等;然后将原来基于CPU串行的元胞自动机模型进行任务划分,得到串行计算部分和并行计算部分的任务,其中并行计算部分主要是将计算密集和能够并行的部分利用元胞自动机的并行特性映射到GPU中,基于CUDA并行计算构架,通过采用功能和数据分解、线程映射等并行计算方法以及存储器访问优化、通信优化等优化策略,进行某一时刻下基于GPU的溃坝洪水演进模型并行计算;计算完成后,将GPU计算结果映射回CPU内存并进行结果统计,并对每个元状态进行相应CPU中的计算,并行计算则交给GPU负责,以提高溃坝洪水演进模型的计算效率,解决在利用高精度DEM进行溃坝洪水模拟的计算瓶颈问题。

图 1 GPU-CA的溃坝洪水演进计算模型 Fig. 1 Dam-break Flood Routing Computing Modelof GPU-CA

图选项

CUDA提供了细粒度数据并行度和线程并行度,可嵌套在粗粒度数据并行和任务并行中。并行粒度分为三个层次,分别为线程(thread)、线程块(block)和线程网格(grid)。线程是并行程序的基本构建块^[25],线程块可包含一定数量的线程,以一维、二维和三维的方式对线程进行索引。线程网格可包含一定数量的线程块,也可以一维、二维和三维的方式对线程块进行索引。本文采用二维元胞自动机模拟框架,采用二维的格网和线程索引。图 2为一个线程块包含(16,16)线程的示意图,一个线程网格包含多个线程块,并且以二维空间进行分布,同时一个线程块又包含多个线程,一个线程与元胞空间中的一个元胞相对应,负责该元胞的计算。线程和元胞索引的计算方法如式(3),其中i、j为元胞的索引位置,blockIdx、blockDim、threadIdx均为CUDA的内建变量。通过核函数调用kernel_function<<>>(param1,param2,…)来启动成千上万的并行线程进行计算。

图 2 元胞与线程的对应关系 Fig. 2 The Corresponding Relation Between Cellular and Threads

图选项

由于各个线程块的计算速度不同,若只采用一个元胞空间,可能造成读写冲突。为了解决这个问题,本文采用两个对应的元胞空间,第一个元胞空间负责数据的读取,第二个元胞空间则负责计算完成后结果的写入,将两个元胞空间状态信息交换进行下一时刻计算。

虽然在元胞自动机模型的扩展时应该尽量保持局部特征,但由于计算过程需要对一些全局信息(比如元胞修复数量、水量总体平衡等)进行统计与修正,而这些全局信息在所有线程中都有可能被修改,因此引发竞争冒险,导致不同的线程在对全局内存的修改是不规则的,从而导致统计错误。CUDA提供了原子操作,这些操作可以确保线程读取/修改/写入操作可以作为一个整体的串行执行^[25],如图 3所示。但是大量的原子操作将会对计算性能造成影响,因此只适用于少量线程发生竞争的情况。

图 3 线程的原子操作 Fig. 3 The Atomic Operations of Threads

图选项

对于另外一些统计值,比如统计所有元胞水深的总和,可以充分利用GPU并行计算的特点,进行并行规约计算,在GPU中新建立内存空间存储每个元胞的水深值,然后使用数据集数量一半的线程,每个线程将当前线程对应的水深与另外一个线程的水深相加。每增加一次,线程数减少一半,直到剩下最后一个线程为止,这个线程中的水深值就是所有元胞水深的总合。

如图 4所示,CPU和GPU协同进行溃坝洪水演进模拟与分析。CPU负责逻辑性强的事务处理和串行计算,首先在CPU端进行初始数据的读取,包括DEM数据、糙率、溃坝位置和溃坝相关信息等;然后进行内存和任务的分配,CPU端的串行计算任务主要包括过程流量线计算、溃口流量计算、剩余库容计算和演进时间计算等,串行计算得到的结果将和初始的数据信息一起映射到GPU端的显存中。

图 4 CPU/GPU协同的溃坝洪水演进模拟与分析 Fig. 4 GPU/CPU Collaborative Implementation for Dam-break Flood Routing Simulation and Analysis

图选项

GPU则负责执行高度线程化的并行处理任务,首先根据水深、总流量等数据计算溃口处的元胞流量;接着计算每个元胞的X方向单宽流量和Y方向的单宽流量;然后根据新的单宽流量计算新的水深,通过对水深进行统计来调整新的水深,各个GPU端的核函数之间通过CUDA内置的_syncthreads ()函数实现线程块内的线程同步;最后将计算结果映射回CPU端,得到水深、流速、淹没面积等数据,一方面可以进行下一时刻元胞状态值的更新,一方面也可以将其用于实时溃坝风险分析与可视化展示。

本文试验硬件环境为DELL T7610图形工作站,CPU的处理器为2×Intel Xeon (至强) E5-2670 2.5 GHz CPU,64 GB主机内存,显卡为NVIDIA公司的Quadro K4000,显存为3 GB,操作系统为Windows 7 64位专业版。GPU计算环境为CUDA5.5,开发环境为Microsoft Visual C⁺⁺ 2012,研发了用于GPU-CA模型下溃坝洪水演进实时模拟与分析的原型实验系统平台。

实验的案例如图 5所示,为肖家桥堰塞湖(104.26°E~104.294°E,31.635°N~31.667°N),距离晓坝镇7 km,流域面积154.81 km²,坝高57~67 m,坝长260 m,坝宽390 m,滑坡方量约242万m³,最大蓄水量3 000万m³。若发生溃坝,洪水将沿着茶坪河对晓坝镇、桑枣镇、安昌镇等近10个乡镇的200多个村及社区造成不同程度的影响。洪水影响区域内存在多个乡镇中心、一条河流及几条道路。

图 5 案例区域 Fig. 5 The Case Region

图选项

原始地形格网大小为10 m,为了进行计算效率的比较,将其进行重采样为10 m、20 m、40 m和60 m格网大小的栅格数据,并分别在CPU和GPU上进行溃坝洪水演进的模拟计算。并行计算的性能分析采用加速比的概念^{[17, 26]},即通过并行计算时间与CPU串行计算时间的比值来体现加速性能。对CPU-CA模型和GPU-CA模型的计算时间进行记录,每完成一次计算记录一次时间。从表 1可以看出,使用GPU-CA模型的计算效率得到了显著提升,尤其是元胞网格分辨率精细时更为明显。当元胞大小为10 m时,加速比为15.9倍,当元胞格网大小较大时,加速比小。这是因为数据传输时间与计算时间相比占了比较大的比例,因此,GPU-CA在大数据量情况下能获得更好的并行加速效果。

GPU-CA模型下溃坝洪水演进实时模拟与分析的原型实验系统平台,可以直观地展示溃坝洪水的演进过程,如图 6所示;并可以实时地查询不同时刻的水深、流速、淹没范围等信息,支持不同时刻下的溃坝风险评估(风险人口、风险城镇和重用公共设施等)。表 2为GPU-CA溃坝洪水模型和基于CPU-CA溃坝洪水模型进行风险评估的时间对比,计算范围为溃坝口到香溪村15 km的河段。当元胞格网大小为10 m时,基于CPU-CA的溃坝洪水模型完成风险评估需要40.89 min,而GPU-CA的溃坝洪水模型完成风险评估只需要2.59 min,将风险评估的效率提高了15.78倍,进而可在溃坝洪水演进模拟过程中实时地进行分析与风险评估。

图 6 溃坝洪水演进过程可视化模拟与分析 Fig. 6 Visualization Simulation and Analysis of Dam-break Flood Routing Process

图选项

与传统的基于CPU串行和并行的元胞自动机相比,GPU通用计算具有计算成本低,计算密度高,安装部署容易,计算加速比随计算规模提升等优点,更适合完成地理模拟和地理计算中的密集计算任务^[17]。本文基于元胞自动机的并行特性和GPU计算密集、高度并行的特点,设计了GPU-CA的溃坝洪水演进计算模型,并讨论了溃坝洪水演进元胞自动机模型,GPU模型映射,计算优化,CPU/GPU协同溃坝洪水演进模拟与分析等关键技术,并通过研发原型系统,开展案例测试与分析。实验结果说明,本文提出的方法能够有效地提高溃坝洪水演进模型计算效率,可实时地支持溃坝洪水发展过程的可视化模拟分析和风险评估。

由于基于元胞自动机的溃坝洪水演进计算模型需要考虑的因素复杂且相互制约,需要在计算过程中处理较多的特殊情况,导致并行计算程序中产生了较多分支,影响了并行线程的计算效率。因此,未来工作将重点对溃坝洪水演进模型涉及的多制约因素进行深入梳理,进一步减少循环判断分支,通过并行算法的改进与优化提高模型的计算效率。此外,结合当前较热的WebGL、HTML5等技术和基于GPU-CA的溃坝洪水演进模型实时计算算法,部署成网络服务,开展面向多用户终端的网络时空过程模拟与可视化分析研究,也是下一步工作的重点。