一种基于GPGPU的指控系统空间查询优化方法

赵艳伟 于璐 张宁 刘凯 徐海旭

指控系统作为武器系统装备的一个重要组成部分,是整个武器系统的信息处理中心,也是指挥员了解战场态势、进行指挥决策的重要手段.而要在复杂的多维空间战场进行指挥控制,必须处理极为庞大和繁杂的战术信息.其中60%与地理背景密切相关[1],尤其在野战过程中,指挥决策对于地理信息及其相关地理属性的依赖性更为突出,地理信息直接影响作战方案的制定和战术思想的落实.通过指控系统的GIS功能,可实现对战场地理空间的实时综合处理和直观可视化分析,形象表示敌我态势,从而辅助指挥员进行战场指挥决策.

在指控系统中,最为常用的GIS功能就是战场电子地图的浏览与战场空间查询,根据查询条件,指挥员可以找到一定地理范围空间内的目标数据,例如战场军事标注点、远程分队的覆盖区域等.

随着网络化、一体化作战模式的发展,指控系统需要处理的各类数据越来越复杂,数据量越来越庞大,仅地理信息数据就能达到TB级以上,现有的指控系统在面对如此密集的计算任务时,传统的处理方式已无法满足应用对大规模计算的实时性需求[2].当面对海量战场空间数据的查询与分析等复杂处理任务时,数据处理周期通常需要控制在几百毫秒乃至数毫秒内,而传统基于树型索引的空间查询方法,当数据越来越集中时,空间要素的重叠度也将越来越高,过高的重叠必然导致空间检索过程中出现冗余路径,严重影响了检索效率,现有的GIS空间查询功能已无法满足指控实时性的需求.

从当今计算机科学发展的前景来看,多机集群、分布式、并行处理将成为应对大规模复杂计算模式的主流趋势,鉴于指挥控制系统的军事应用需求及应用特点,采用并行处理架构将成为必然.GPGPU(General Purpose Graphic Processing Unit)凭借其大规模并行计算能力和易编程性已成为高性能计算和科学计算应用的重要加速手段,并广泛应用于其他通用计算领域.GPU通用计算为复杂性超出CPU解决能力和计算需求的应用,提供了一个使用较少能耗却能发挥更高性能的解决方案[3].

在指控系统对密集型计算任务的实时处理需求下以及GPU通用计算迅猛发展的背景下,本文提出了一种基于GPGPU的指控系统空间查询优化方法.

本文首先介绍了GPU通用计算的相关理念,在此基础上重点设计了一种可应用于指控系统GIS分系统的CPU+GPU并行空间查询方法,并通过实验数据证明了本方法的高效性,为未来实现GPU通用计算在军事领域上的应用提供可行的技术参考.

1 GPU通用计算

图形处理器(Graphic Processing Unit,GPU)从诞生之日起就以超越摩尔定律的速度发展,运算能力不断提升.业界很多研究者注意到GPU进行计算的潜力,于2003年SIGGRAPH大会上提出了GPGPU的概念.GPU逐渐从由若干专用的固定功能单元(Fixed Function Unit)组成的专用并行处理器向以通用计算资源为主,固定功能单元为辅的架构转变.GPU无论在浮点运算能力上还是存储带宽上都远远超过了同时期的主流CPU,GPU通用计算是一片正在被打开的潜力巨大的市场[4].

近几年,越来越多的学者开始把对并行计算计算性能的提升关注重点从CPU转移到GPU上来,GPU内部集成了数百个并行流处理器,其并行计算效率比CPU要高很多,而且已经有很多学者和厂商将GPU应用于科学计算领域,取得了非常好的效果[5].例如犹他州大学在医学成像领域使用GPU可获得146倍加速效果;马里兰大学在基因排序领域使用GPU可获得30倍加速比,除此之外,GPU在分子运动、视频转换、量子化学、线性代数等领域都有不同程度的性能提升[6].

目前,在国内,GPU应用在军事指挥控制系统上主要用来处理图形图像等传统问题,针对计算密集型的指控信息计算,尚无成熟的基于CPU+GPU进行加速的指控功能的相关研究.

2 基于GPGPU的并行空间查询方法

通常,为了提升空间查询的效率,首先需要为空间数据建立索引机制,其中,树型索引(R树、R*树[7]、R+树等)是采用最为广泛的一类索引结构.由于树型索引在空间查询时需要大量使用指针及递归操作,这样的操作方式在显存中是不支持的,在GPU环境下,传统的树型索引并不适用.因此,本文设计了一种适用于GPU环境的并行栅格索引.

另外,为解决由于过滤和精炼过程带来的CPU与GPU频繁进行数据交换的问题,本文在栅格索引基础上,设计了一种CPU与GPU协同加速空间查询方法,将传统的粗糙过滤和精炼两个空间查询步骤进行合并,尽可能地最大限度过滤非相关战场空间要素,避免了内存和显存的频繁数据交换,提高了空间查询的整体效率.

2.1 适用于GPU环境的并行栅格索引的构造

为了最大化地利用GPU的多线程机制,本文考虑采用栅格结构建立空间索引,一方面利于计算任务均匀划分,另一方面可以极大地保证处于激活状态的线程块数量.

整个栅格索引渲染过程算法框架如图1所示.首先,输入矢量要素信息,该信息包括组成矢量要素的点坐标,并通过坐标转换通道进行坐标转换,将矢量要素坐标精确到亚像素精度.然后进入算法核心,该核心包括两部分,第1部分(图1中左侧部分)是矢量要素栅格化过程,第2部分(图1中右侧部分)是要素索引与对应栅格的解析过程.

要素渲染过程分为3步,第1步,提取矢量要素的轮廓信息;第2步,利用栅格控制器,扫描轮廓并计算填充格网区域;第3步,按照水平扫描的方式对填充区域进行跨段或跨单元填充.索引解析过程,通过输入的矢量要素获得索引信息,通过索引转换引擎将其映射到对应的栅格单元中,最后与矢量要素栅格化的结果通过渲染缓存生成栅格索引.

2.2 索引在主机端与设备端的传输与替换

由于显存空间的大小受限,因此,对于大比例尺下的空间图层要素,其栅格索引无法一次性全部加载到显存空间中.为了解决这个问题,需要将栅格索引进行瓦片分割.采用预生成及预装配技术,在初始化时,按照显存内分配的索引空间大小将部分瓦片索引预传输到设备端.如果查询请求的范围内无法获得有效的显存栅格索引,则采取瓦片替换策略对索引进行更新.

由于地理空间的连续性,用户在查询某一个区域的时候很可能再访问其相邻的区域,因此,在替换更新时需要对相邻区域进行考虑.若令索引瓦片为TileIndex(i,j),其中i代表索引所在行,j代表索引所在列.则该瓦片的相邻区域采用如图2表示,形式化定义为:

图1 并行栅格索引生成过程

图2 领域瓦片范围

每次访问显存内的索引瓦片,其对应位置的计数器都会进行累加.若当前访问的索引TileIndex(i,j)不在显存内,则按照LFU替换策略,利用CUDA并行归并排序算法统计访问频率,最先移出最少使用的索引瓦片TileIndex(s,t).其中需要多考虑的一个步骤是,如果TileIndex(s,t)∈Ωi,j,则不进行替换,按照索引访问频率选取下一个待替换的瓦片.

2.3 并行空间查询算法的设计与实现

利用上述这种适合在GPU环境下实现并行空间查询的栅格索引机制,这部分主要讨论两类最常用的空间选择查询:BBOX范围选择查询和多边形选择查询.这两种查询在GPU设备端的并行查询的区别在于,多边形选择查询需要额外再多计算一步,即将不规则的查询多边形栅格化,并且在多线程并行计算时,每个线程需要多进行一次栅格判断.其余整体算法步骤可按照BBOX范围选择查询来处理.因此,本部分先主要介绍BBOX范围选择查询的具体算法,再对该算法进行扩展,从而能处理任意多边形选择查询请求.整体流程如图3所示.

2.3.1 索引在显存内的存储模式

考虑到GPU多线程对显存空间访问的特性,索引在设备端的存储需要满足合并访问的原则.采取的具体方法如图4所示.

首先,对于给定的显存大小,确定待存储索引的范围,若假设包含m行n列瓦片,则应存储的索引集合为然后对各个瓦片之间,采用Hilbert曲线的方式进行存储,其中指定的瓦片在显存空间内的起始位置可通过(H−code−1)×size快速定位.最后针对每个瓦片内部,采取行扫描的方式存储,并且保证行宽度为(W1,W2,W3,W4)字长的16的倍数.

采用这样的方法主要基于两点考虑:由于相邻瓦片索引在实际的地理空间上具有相近性,而Hilbert曲线可以将相邻索引进行聚集,由此在各个瓦片之间极大地保证了空间数据良好的局部性,使得在空间查询时能够访问显存内连续的一段空间,为合并访问打下基础.另外对于每个瓦片的内部,上述方法能使得Half-warp中的16个线程只经过一次传输便可以完成整个Half-warp对显存的访问请求.经过上述步骤处理后,显存空间的瓦片索引既实现了线性存储,同时也满足了合并访问的原则,从而实现合并读取,提高线程对显存的访问效率.

图3 算法流程

图4 索引瓦片在显存内的存储方式

2.3.2 线程块计算过程

这个过程主要计算参与并行空间查询的线程块数量和每个块内的线程数量.其中线程数量以二维方式指定,且保证每个块内线程总数不超过512,设为TN(32,16).为了保证每个线程处理一个栅格单元,需要按照查询请求的MBR范围来计算线程块数量.

假设查询范围为Bbox(x1,y1,x2,y2),且每个索引瓦片的宽为idwidth,高为idheight.那么此范围覆盖的瓦片索引集合为SIdx(s,t),其中,

2.3.3 kernel并行查询过程

根据BN(bnx,bny)和TN(32,16)的线程启动策略,通过3层偏移关系可以计算出线程处理的空间位置.首先,在设备端并行空间查询的核函数中可利用threadIdx内建变量,确定当前线程在线程块中的偏移位置,然后再根据Grid中Block-Dim、BlockIdx内建变量计算出线程在整个Grid中的全局位置,最后根据Grid覆盖的起始瓦片在显存索引中的偏移来进一步确定当前线程应处理的单元的空间位置.空间位置的确定如图5所示.

若该空间位置在查询外包内,则需要对栅格进行提取.否则,当前线程不参与运算.对于Bbox范围选择查询,可按上述方式过滤非相关线程的计算;对于多边形选择查询,需要按照查询条件的外包进行二值栅格化,即多边形范围内栅格化为1,外包内且多边形外部栅格化为0.线程在并行提取栅格单元时,需与对应位置的查询条件进行栅格运算,如图6所示.

图5 三层偏移计算关系

图6 不规则范围查询处理情况

显存内对应位置栅格单元的具体提取过程:首先通过上述步骤得到的每个线程处理的空间位置,计算出该空间位置所属的索引瓦片,根据瓦片Hilbert曲线排序对应的H-code进一步确定该索引在显存内的偏移.由于瓦片内部采用行扫描顺序存储,因此,只需根据空间位置在瓦片内部的偏移按照行列顺序获取具体的栅格单元.

影响CUDA性能的一个重要因素是全局存储器的访问效率[8],由于在多线程提取索引单元时需要对显存进行大量的读写,因此,介绍一下本方法在对全局存储器读取进行的优化,这也是生成栅格索引采用四字节的原因.在构造栅格索引结构时,采用三字节这种存储方式对于中等数据量的空间信息来说能够保存的索引数据已经足够,但如果是这样的话无法保证同一half-warp中的线程访问显存的地址对齐到可合并访问的段内,因此,将导致数据分为16次串行传输.如图7所示.

图7 合并访问方式

由于CUDA支持32 bit字长的合并访问,第1种情况,索引的(W1,W2,W3)字长无法对齐到可合并的段内,造成线程访问显存的效率较低.而第2种情况,虽然在存储时多了一个字节,但却正好满足全局存储器内存融合访问模式,索引与合并段一一对齐,因此,half-warp中16个线程访问的数据只需进行一次传输,优化了访问效率.如图8所示.

图8 并行排除冗余索引

采取的具体做法是:利用CUDPP[9]的并行归并算法对提取的有效索引ID进行排序,对排序后的索引集合使每个线程分别处理一个元素,计算该元素与其前驱的差值并保存为标识量flag.这样,flag中不为0的元素位置,对应的索引集合中的矢量要素ID一定是不重复的.相反,flag中为0的一个或者连续一段元素对应的索引ID说明为冗余值,可忽略不计.因此,根据flag的标识性质,对索引集合利用CUDPP compact进行并行压缩,将标识量非0的要素ID前置,最后利用并行Prefix sum方法求得最终ID结果集,只将该部分集合传递回主存中,这样便可以大大减少传输量,提高并行查询效率.

本文所设计的并行空间查询方法较一般空间查询的区别如图9所示.第1种情况是使用树型索引的传统两步空间查询,即根据查询条件W的MBR,先进行空间过滤,筛选出候选集合R1和R2.再通过精炼步骤,对空间要素进行精确检测,找到查询结果为R1.第2种情况是并行栅格查询方法,只要根据查询条件W通过对栅格索引对应位置的并行提取,便能够确定查询结果集.这种方法合并了原始方法的过滤和精炼两步,既避免了候选集在显存与主存中的频繁传输,也避免了精确匹配时大量的计算过程.

图9 传统空间查询步骤与本文方法对比

3 实验结果与分析

实验部分主要分为两个部分,第1部分主要将本文所述的基于GPGPU的并行空间查询方法与传统CPU中的空间查询做对比,第2部分主要分析主机端与设备端之间数据传输对本文算法的整体性能影响.实验环境配置如下:CPU处理器采用InteCoreTMi7-9202.67GHz,内存为6GB;GPU显卡采用GeForce GTX 260(192 SP,1.3),显存为2GB.利用CUDA 4.0[10]工具集在Microsoft Visual Studio 2008进行程序开发.

3.1 与CPU方法对比

这部分实验的目标是测试本方法与采用R*空间索引的CPU中的空间查询方法,在相同benchmark测试集及相同查询请求情况下的空间检索时间及空间过滤度.测试集主要是针对点、线、面3类不同类型的地理要素,为保证查询的公平性及均匀性,空间查询请求通过随机生成,并且每个测试图层按照不同的查询范围进行8组实验(每组包含500个查询请求),查询范围则通过与整个图层的比例关系进行确定.例如查询结果横坐标120表示将原始图层范围的1/120作为查询条件.

对于稠密的点数据,虽然本文方法的计算时间较R*方法少,但若计算数据传输时间,与CPU中采用R*索引的空间查询时间相当,如图10所示.

图10 点图层查询结果

对于海量密集的线图层和面图层来说,本方法在空间查询时间方面较CPU中基于R*树型结构索引的空间查询具有绝对优势.由于数据量较大,树型索引的高度将急剧增加,又由于数据密集,导致索引中间节点的MBR重叠度较高,因此,树型索引在处理这类空间数据上效率较低.从图11中可以看出,对于每组数据采用R*索引的查询时间平均超过了1600ms,而本文提出的方法即使计算数据的传输时间也可以控制在50ms左右,此时数据传输时间与CPU中的计算时间相比可忽略不计,因此,在本文提出的方法在处理密集线数据上可以得到32倍的加速比.同样的,从图12中可以看出,由于密集面数据具有较高的MBR重叠度,采用R*索引平均需要500ms左右,而本文提出的方法需要50ms左右,其加速比也能达到10倍左右.

从以上对不同类型及不同特点空间数据的查询效率对比分析中,我们可进行如下总结:本文提出的基于栅格索引及GPGPU并行架构的空间查询方法适用于处理海量密集的线面类空间数据,与CPU中的树型索引方法相比,能达到量级的效率提升.对于传统CPU中的树型索引查询方法,数据量越大数据越密集,查询效率越低,而本文提出的方法的查询效率不受数据复杂度的影响,只受查询范围及查询到的ID集合大小影响.

图11 线图层查询结果

图12 线图层查询结果

3.2 kernel性能分析

这部分主要针对本方法的内部执行过程中,各个部分所占的时间比进行测试,意在分析主机端与设备端之间数据传输对整个算法的性能影响.我们采用的分析工具为CUDA profiler[11],它是NVIDIA公司在CUDA toolkit中提供的一种可视化工具,其中包含了一系列可以反映应用程序中CPU和GPU的活动行为的时间表,同时CUDA profiler还提供了自动分析引擎来帮助用户迅速确定性能瓶颈位置,从而找到可优化的策略.如图13所示.

通过Summary Plot图表,我们可以直观地看到kernel函数和memory copy在总GPU时间占用的百分比.不难发现,本文提出的方法窗口查询的计算部分占到了整个kernel时间的76.21%,而数据拷贝只占到了20%左右,因此,可以得出结论,本文提出的方法中数据传输部分并未造成性能瓶颈.

4 结论

图13 CUDA profiler性能分析

为了满足指挥控制系统对战场空间目标的实时查询与直观展示,本文提出了一种基于GPGPU的并行空间查询方法.首先,设计了一种基于栅格结构的并行空间索引机制,相比于树型索引,本文设计的索引结构具有更小的空间重叠度,同时能够充分利用GPU并行环境的计算特点,为细粒度的并行查询提供了有利的契机.在此并行栅格索引机制的基础上,本文提出了一种在GPGPU计算模式下的并行空间查询算法.通过与已有的研究成果进行试验对比,在执行效率方面,本方法对于处理海量密集型的空间数据具有绝对的优势,加速比可达到一个数量级.

本文提出的方法可应用于各级指挥控制系统的GIS地理信息分系统,可有效提高战场空间目标的检索效率,满足毫秒级空间查询响应时间要求,为辅助战场指挥决策起到了关键作用.该方法的研究,也为我军指挥控制系统采用CPU+GPU异构计算模式进行指控功能算法加速提供了一种探索性的尝试,为促进我军指控系统的跨越式发展打下了坚实的基础.