1、基于 GPU 栅格计算关键技术研究及实现 卢涵宇 陈军 卢天健 陈秀万 纪翔 北京大学遥感与地理信息系统研究所 贵州大学大数据与信息工程学院 成都信息工程大学资源环境学院 中南大学资源与安全工程学院 摘 要： 针对 CPU 处理栅格计算中的大量并行运算效率不高的瓶颈, 选择高并发计算能力的 GPU 来解决该问题。通过分析 GPU 通用计算的关键技术, 以 CUDA 为编程环境, 分别实现了 GPU 中简单栅格计算、基于窗口运算的邻域分析和以汇流计算为例的栅格迭代分析。实验结果表明:CPU 计算的时间消耗与输入数据长度和数据量线性正相关, GPU 并发能力容纳范围内随着输入数据长度幂数增长和数据
2、量增加, 时间消耗变化不明显。证明了 GPU 栅格计算的有效性, 具有一定的理论意义和实践价值。关键词： GPU; 并行计算; 栅格运算; CUDA; 作者简介：卢涵宇 (1978) , 男, 河南周口人, 贵州大学副教授, 博士后;E-mail:。收稿日期：2017-09-07基金：国家自然科学基金资助项目 (41671355) Research and implementation of key technologies for grid computing based on GPULU Han-yu CHEN Jun LU Tian-jian CHEN Xiu-wan JI Xiang
3、Institute of Remote Sensing and Geographical Information System, Peking University; Institute of Resources and Environment, Chengdu University of Information Technology; School of Resources and Safety Engineering, Central South University; College of Big Data and Information Engineering, Guizhou Uni
4、versity; Abstract： Aiming at the bottleneck of inefficiency of CPU in parallel processing grid computing, GPU which has high concurrent computing power is chosen to solve the problem in the paper. Based on the analysis of the key technology of GPU general-purpose computing, using CUDA as programming
5、 environment, simple grid computing, the neighborhood analysis based on window operation and grid iterative analysis of running water accumulation quantity computation are realizedrespectively. Experimental results show that the time consumption of CPU is linear related to the length of the input da
6、ta and the data quantity, and the change of time consumption is not obvious when the length of the input data grows exponentially and the data volume increases. It proves the validity of the GPU grid computing. The work has certain of theoretical significance and practical value.Keyword： GPU; parall
7、el computing; raster computing; CUDA; Received： 2017-09-07栅格计算模型众多, 支持大量的空间数据分析, 在地理信息系统 (geographical information system, GIS) 数据处理中应用广泛1-3。传统数学分析模型是针对中央处理器 (central processing unit, CPU) 环境设计的, CPU 具有强大的逻辑处理能力, 当大规模并行计算时执行效率往往成为应用的瓶颈4-6。图形处理器 (graphic processing unit, GPU) 是为了处理规模庞大又互不相干的相同模式的数学计算
8、。国内外学者对 GPU 计算及其应用开展了大量研究7-8。GPU 通用计算技术可实现更复杂的模拟, 逐渐与图形表达解耦, 由程序直接调用 GPU 硬件资源完成数学计算任务。栅格计算按照某种方法逐像元计算得到新的栅格数据, 其最小计算为栅格像元。一般情况下, 这种基本计算单元可设计为互不相干的并发单元, 因此栅格计算非常适合移植到GPU 中来提升计算性能。统一计算设备架构 (compute unified device architecture, CUDA) 是显卡厂商 NVIDIA 公司推出的 CPU 编程模型, 它采用了新的网格模型来抽象化 GPU 资源, 作为 GPU 通用计算工具广泛应用
9、到各领域的海量并发计算中9-13。在 GIS 空间分析领域, 成熟的主流商利用 GIS 解决方案中针对栅格计算模块尚未大规模转为硬件加速模式的问题。栅格计算的基础引擎还在大量使用 CPU 运算, 且多线程利用率并不高。以 CUDA 为编程工具的 GPU 栅格计算, 无疑是 GPU 应用于GIS 领域进行的一次有益尝试, 具有重要的理论意义和实践价值。1 GPU 通用计算技术与栅格计算1.1 GPU 通用计算关键技术1.1.1 GPU 通用计算方案相较于 ATIStream、NVIDIA、Open CL 和 Direct X, CUDA 硬件资源的逻辑化更加适配传统栅格计算数学模型, 有足够直观
10、的底层 API 调用提供, 已经具备完善的线程管理和调度机制。因此本文采用 CUDA 技术框架实现 GPU 栅格计算。1.1.2 CUDA 网格计算模型中 GPU 的逻辑模型CUDA 将 GPU 看作可以承担若干计算线程异步执行的硬件设备。GPU 执行的代码称为核函数。GPU 硬件在逻辑上分为若干层, 每一层有若干个处理器集群, 每一个集群中有若干个流多处理器 (streaming multiprocessor, SM) 。SM 是一个线程池, 在同一个 SM 中执行的线程共享数据。一个 SM 中可容纳 8 个 Block (线程块) 或 1 024 个 Thread (线程) , 每个 Bl
11、ock 可以容纳 512 个线程, 这些线程在 Block 中以 Warp (线程束) 的形式被管理。1.1.3 CUDA 设计模型中线程号索引机制激发核函数的调用参数 Dg, Db, Ns, S, 根据不同的参数配置可实现不同类型的索引机制。CUDA 模型中线程通过线程块分割, 每个线程块拥有一个块 id。如果一个线程块只分配一个线程, 使用 N 条线程块容纳 N 条线程, 线程索引号等于线程块索引号。可以将线程集合看做一个列向量, 向量的每一行表示一个块, 块中只有一个元素。此时线程块 id 的 x 索引值 (block Idx.x) 等于线程 id (Tid) , 如下式:如果使用 MN
12、 个线程块容纳 (MN) 条线程, 每个线程块依然只容纳一条线程, 则将线程块的分布展开到二维空间上, 线程 id (Tid) 等于线程块 id 的 x 索引值与线程块 id 的 y 索引值乘线程块定义 x 方向长度乘积的和。计算如下式:如果一个线程块容纳多条线程。即使用 MN 个线程块, 每个线程块容纳 PQ条线程, 即一个线程块索引 Tblockid, 假设 P=Q=3。Tblockid 内部分布如下:可继续分解 Tid 到 x 和 y 两个方向。则线程索引也具有两个维度, 各维度计算方式如下:这样就可以把二维矩阵上的像元映射到线程索引上。计算方式如下:其中, Grid Dim 表示网格矩
13、阵的定义;block Dim 表示线程块的定义。Grid Dim.x=N, 表示网格划分为 N 列;Grid Dim.y=M, 表示网格划分为 M 行;block Dim.x=Q, 表示区块划分为 Q 列;block Dim.y=P, 表示区块划分为 P 行。通过多个维度的线程索引号计算偏移量, 各个线程即可独立访问各自管理的内存空间。线程索引在 GPU 栅格计算过程中作为栅格像元的映射机制在每一个核函数中使用。核函数中所调用的参数 i 均是由此计算所得, 该值配合激发核函数时配置的激发参数。1.1.4 GPU 编程和 CPU 编程的主要区别二者的区别主要有两点: (1) GPU 编程本质上是
14、并行程序设计, GPU 拥有远强于CPU 的并发处理能力。 (2) 并行计算过程的逻辑无关性和跨线程性要求对传统CPU 算法逻辑重构, 以快速、正确的得到结果。针对上述区别, 对 GPU 栅格计算进行分类处理, 合理利用优势和规避劣势。1.2 GPU 栅格计算分类栅格计算按照计算过程的时间复杂度和单元逻辑相关性复杂度可分为三类:简单栅格单元计算、栅格邻域分析和栅格迭代计算。1.2.1 简单的栅格单元计算栅格计算器可完成诸如植被指数计算、栅格叠加分析等包含了两组或更多栅格数据之间的多元运算, 其过程本质是栅格数据中对应位置的各个栅格单元的重复代数计算。简单的栅格单元计算主要通过合并算子减少内存拷
15、贝的时间开销, 提高整体计算效率。1.2.2 栅格邻域分析栅格邻域的并发计算过程中, 不仅需要索引到计算位置的栅格, 还需要索引符合邻域表达式的该栅格的周边栅格。核心问题是将内存中数据区块的索引机制设计为支持二维栅格索引的方式, 并在一定程度上简化邻域栅格的重复求解过程。1.2.3 栅格迭代计算栅格迭代计算需要对同一栅格多次迭代计算且每次计算依赖于上一次的结果, 典型的如汇流计算等。每一次迭代计算是各个栅格数据并发计算, 栅格迭代分析计算量大, 需要利用 GPU 大规模的并行计算能力。2 简单的栅格单元计算的 GPU 实现简单的栅格单元计算中同时有加法和乘法时必须进行算子的合并操作。算子合并在
16、 CUDA 中通过 Thrust 开源库的 Transform 函数实现。2.1 栅格和计算CUDA 编程实现栅格 A 和栅格 B 相加, 需定义核函数 add Kernel () , 并由标识符_global_修饰。函数接受 3 个都是指针的参数。核函数内部通过该指针与内建的线程索引机制实现计算数据的索引。参照式 (1) , 做简单映射可得:调用该核函数需要定义与索引方式匹配的激发参数以完成线程分配。激发函数负责检索输入参数, 并完成设备端内存的创建和拷贝。激发核函数后等待线程同步, 并返回计算结果。transform () 函数利用参数完成动态的内存分配与核函数激发, 核函数的具体执行动作
17、由泛型函数对象 thrust:plus () 指定。这个函数对象即规约算子。规约算子的数学意义即为一个单步操作逻辑, 它在 CUDA 底层 API 中代表着一个核函数的实现。2.2 栅格的 SAXPY 表达式执行一次二元相加操作所使用的规约算子是 thrust:plus () 。通过适当运用 Thrust 可以提高计算效率。以 Thrust 实现 SAXPY 表达式合并计算为例:在 Thrust 中分布计算过程有 4N 次读取和 3N 次写入。使用单一规约算子能将计算过程简化为 2N 次读取和 N 次写入 (省略掉了数组 a 的填充和一次 transform函数调用) , 以及减少一次线程同步
18、。而合并算子后的计算过程跳过中间步骤的线程同步和内存拷贝操作, 节省大量的时间。3 栅格邻域分析的 GPU 实现栅格邻域分析过程如图 1 所示。图 1 中 33 的深色数据为邻域窗口, 浅色栅格为当前计算栅格。计算过程中邻域窗口将被用于匹配计算像元的周边参与计算的像元。图 1 窗口计算原理 Fig.1 Window calculation principle 下载原图栅格邻域分析在数学上是矩阵的窗口运算。并发过程中只会多次读取原始数据中计算栅格的邻域栅格数据, 不存在跨线程数据的写入操作。邻域分析过程中只涉及到跨线程数据读取操作, 核心问题在于高效实现邻域的索引。3.1 二维数据的空间邻域关系
19、求解基于复杂线程号映射机理, 可以实现相邻栅格的空间索引。有二维数据映射:则距离计算像元距离在 N 以内的像元为:3.2 邻域解的复用计算关键在于邻域窗口的复用, 在核函数激发前先计算出无偏移量的所有领域范围索引, 并将其作为参数传入核函数。核函数中仅通过基础邻域加上自身偏移量即可匹配到自身的领域, 以节省重复邻域计算带来的开销。该核函数不实现邻域求解, 其所指数据的邻域由给定邻域参数*b 和邻域长度*ns 加上自身偏移量 i 解得。4 GPU 支持下的栅格迭代计算4.1 迭代计算的 GPU 实现流程以流水累积量计算模型阐述栅格迭代计算的 GPU 实现流程: (1) 初始化计算数据, 实例化流
20、水方向矩阵和水源供应矩阵; (2) 按照帧计算模型进入下一帧 (第一帧) , 激发一次核函数; (3) 如果是初次进入 Device 端计算, 则先完成内存拷贝和线程分组; (4) 完成一次计算; (5) 校验结果, 校验未通过则迭代步骤 2。4.2 汇流计算的栅格迭代计算汇流计算的核心理念是反复计算每一栅格的流入水量和流出水量, 直到水源供应矩阵中所有栅格均被计算完全。定义这个迭代过程为函数 f () 。拆分为 f1 () fn () 。计算流程如图 2 所示, 浅黄色代表该栅格有新水源流入, 浅蓝色表示仅流出。图 2 汇流部分计算过程 Fig.2 Calculation process o
21、f flow accumulation 下载原图4.3 汇流计算测试汇流计算测试数据检测结果见图 3。其中原始数据由伪随机算法生成, 通过简单的取模运算创建流水方向矩阵。初始矩阵 (Base) 中的“1”表示该点自身能为目标点独立提供的汇入水流量, 结果矩阵 (Result) 中的数字表示该点在演算过程中汇入水流量的总和。方向矩阵 (Direct Matrix) 表示该点的目标方向, 方向为八连通区域中的一个位置。图 3 输出测试数据检测结果 Fig.3 Test results of output test data 下载原图5 模型测试及数据对比分析针对本文设计的三种栅格计算的 GPU 实
22、现, 拟 5 组不同长度的数据多次测试。每组每次随机生成, 记录时间消耗平均值, 绘制图表进行对比。表 2 是简单的栅格单元计算、栅格邻域分析两类栅格计算的测试数据评估。从表 2 中可以看出, CPU 计算过程的时间消耗与输入数据长度呈现显著的线性正相关, 在 GPU并发能力容纳范围内输入数据长度幂数增长, 但时间消耗增加不明显。当数据输入达到 E7 时 CPU 的时间消耗已经明显高于 GPU 一个数量级。表 2 GPU 算力消耗与输入值长度的关系 Tab.2 Relationship between GPU Calculate consumption and length of the in
23、put value 下载原表 图 4 是完成栅格迭代计算的 GPU 栅格计算测试。CPU 计算的时间消耗与输入数据量成线性相关。但在 GPU 计算测试中则有显著区别, 当 GPU 算力未达到饱和时, 输入数据量的增长对时间消耗增加不明显。受限于硬件资源, GPU 中规模为 E7 以上的迭代计算测试只能分组实施。GPU 计算的实现就等效于 10 次 E6, 时间消耗开始成线性相关。图 4 汇流计算测试 Fig.4 Calculation test of flow accumulation 下载原图6 结语GIS 中的栅格计算用途广泛, 采用并行功能强大的 GPU, 可有效解决 CPU 中大量并行
24、运算出现效率低的瓶颈问题。本文通过线程 Grid 模型的构造、栅格计算在GPU 并发线程中实施时栅格数据的映射逻辑等的设计, 分别实现了简单栅格计算、基于窗口运算的邻域分析和以汇流计算为例的栅格迭代分析。通过 CPU 和GPU 的计算性能测试, 结果显示:CPU 计算的时间消耗与输入数据长度和数据量呈线性正相关, 在 GPU 并发能力容纳范围内输入数据长度幂数增长和数据量增长两种情况下, 时间消耗增加均不明显。可见 GPU 硬件资源足够的情况下, GPU栅格计算在很大程度上能提升业务执行效率, 能够为 GIS 空间分析提供重要技术支持。参考文献1SHANE C.GPU 编程指南M.北京:机械工
25、业出版社, 2014. 2ROB F.GPU 高性能编程 CUDA 实战M.北京:机械工业出版社, 2013. 3郑方, 张昆, 邬贵明, 等.面向高性能计算的众核处理器结构级高能效技术J.计算机学报, 2014, 37, (10) :2177-2181. 4冷小鹏, 苗放, 杨文晖, 等.基于基于 G/S 模式的三维 GIS 监测系统研究J.广西大学学报 (自然科学版) , 2012, 37 (2) :317-322. 5WEN M, WU H.GPU Computing GemsM.Holland:Elsevier, 2011. 6崔树林, 张旭, 张树清, 等.基于 GPU 的大规模栅格
26、数据分块并行处理方法J.计算机技术与发展, 2015, 25 (3) :19-22. 7SANDERS J, KANDROT E.GPU 高性能编程 CUDA 实战M.北京:机械工业出版社, 2011. 8唐滔.CPU/GPU 异构混合并行的栅格数据空间分析研究D.武汉:武汉国防科技大学, 2012. 9卢敏, 王金茵, 卢刚, 等.CPU/GPU 异构混合并行的栅格数据空间分析研究以地形因子计算为例J.计算机工程与应用, 2017, 53 (1) :172-177. 10商哲然, 谭贤四, 曲智国, 等.基于 GPU 的 RFT 算法并行化J.雷达科学与技术, 2016 (5) :505-509. 11陈丁, 方刚, 王龙, 等.基于 GPU 的二维矢量场 LIC 算法研究J.测绘工程, 2015, 24 (3) :26-30. 12李清泉, 李德仁.大数据 GISJ.武汉大学学报 (信息科学版) , 2014 (6) :641-644. 13洪镇填, 龚根生.基于 GPU 并行计算的巨幅遥感影像坐标转换研究与实现J.测绘与空间地理信息, 2016, 39 (3) :139-142.
