2. 海岛(礁)测绘技术国家测绘地理信息局重点实验室(山东科技大学), 山东 青岛 266590;
3. 公路地质灾变预警空间信息技术湖南省工程实验室(长沙理工大学), 长沙 410004
2. Key Laboratory of Surveying and Mapping Technology on Island and Reed, National Administration of Surveying, Mapping and Geoinformation(Shandong University of Science and Technology), Qingdao Shandong 266590, China;
3. Engineering Laboratory of Spatial Information Technology of Highway Geological Disaster Early Warning in Hunan Province(Changsha University of Science & Technology), Changsha Hunan 410004, China
地貌晕渲图通过模拟太阳光对地面照射所产生明暗程度的不同,采用随光度近似连续变化的灰度或彩色色调表示地面点颜色[1],从而实现以图像形式表达逼真的三维地形起伏,是地形二维表达形式中最具立体真实感的一种可视化表达方式[2]。现有文献从地貌晕渲实现机制出发,对影响地貌晕渲效果的地形因子计算、光照模型及光照模型参数设置等因素进行了深入研究,取得了许多成果[2-5],但现有算法多是基于单核单线程计算模式,计算效率相对较低。目前,高分辨遥感影像、高密度LiDAR(Light Detection And Ranging)点云数据等精细化数字高程模型(Digital Elevation Model,DEM)数据源的获取技术已经相当成熟,这使得DEM表达的精细程度越来越高。随着地貌晕渲区域的范围扩大和表达精度增加,地貌晕渲需要处理的地形网格数显著增加,计算时间消耗更多,对现有单核单线程地貌晕渲算法提出了挑战。
从硬件角度考虑,目前高性能并行计算实现的方式[6-7]主要有计算机集群、多核处理器(Multi-core processors)以及图形处理器(Graphics Processing Unit,GPU)。从软件架构来说,最常用的支持并行计算的并行计算库包括MPI(Message Passing Interface)、OpenMP(Open Multi Processing)、Intel IPP(Integrated Performance Primitives)等。计算机集群模式需要多计算机支持,成本高且开发复杂,仅限于在一些专门领域使用。随着多核处理器以及图形加速卡的普及,基于多核或GPU的并行计算模式正成为未来程序开发的发展趋势。这不仅充分发挥了新的硬件架构的计算性能,也极大提升了应用程序的运行效率。在地学领域,许多学者将单机多核并行计算引入地学计算,作了许多创造性的工作,显著提升了计算效率,如溃坝水流模拟[8]、DEM洼地识别与填充[9]、海量DEM处理[10]、基于坡度的Lidar点云简化与内插[11]、遥感图像处理算法并行化[12]等。
为了适应当前计算机多核处理的发展趋势,微软从.NET Framework 4.0开始集成TPL(Task Parallel Library)和PLINQ(Parallel LINQ)来支持并行,在统一的工作调度程序下进行硬件并行协调,在提高应用程序性能的同时降低了现存并发模型的复杂性。本文将在.NET环境下利用多核并行计算模式对现有地貌晕渲算法进行改进,提高地貌晕渲效率。
1 基于坡度坡向的地貌晕渲算法 1.1 坡度与坡向计算DEM表达主要有不规则三角网(Triangulated Irregular Network,TIN)和规则格网(Regular Square Grid,RSG)两种模型。由于RSG模型具有存储结构简单、适于多分辨率表达、计算处理方便以及易于影像等其他栅格数据融合等优势,而被广泛采用。地形的坡度和坡向的计算方法与地形的表达形式密切相关。当采用规则格网形式表达时,一般采用差分公式来计算。坡度、坡向的计算公式分别为:
$Slope = {\rm{arctan}}\sqrt {f_x^2 + f_y^2} \times 180/{\rm{ \mathsf{ π} }}$ | (1) |
$Aspent = 180 - {\text{arctan}}\left( {{f_y}/{f_x}} \right) + 90 \times {f_x}/\left| {{f_x}} \right|$ | (2) |
其中:fx为X方向上的高程变化率,fy为Y方向上的高程变化率。坡向值按照正北方向为0°,按顺时针方向依次递增,取值范围为0°~360°。
求解fx和fy方法是通过3×3的移动窗口扫描整个格网DEM,利用移动窗口中的格网点高程,按照给定的计算方法来计算中心格网点的fx和fy。常用的计算fx和fy的方法有数值分析法、快速傅里叶变换法、局部曲面拟合法和空间矢量法等。根据移动窗口中格网单元参与计算方式的不同,数值分析法又分为简单差分法、二阶差分法、三阶不带权差分法(Sharpnack算法)、三阶反距离平方权差分法(Horn算法)、三阶反距离权差分法和Frame差分法等方法。这里采用Horn算法,计算公式如下:
${f_x} = \left[ {{z_3} - {z_1} + 2 \times \left( {{z_6} - {z_4}} \right) + {z_9} - {z_7}} \right]/\left( {8g} \right)$ | (3) |
${f_y} = \left[ {{z_7} - {z_1} + 2 \times \left( {{z_8} - {z_2}} \right) + {z_9} - {z_3}} \right]/\left( {8g} \right)$ | (4) |
其中:z1~z9为移动窗口中的9个格网高程,g为DEM格网间隔。移动窗口上各个格网单元的编号如图 1所示。
当移动窗口移动到RSG边缘时,计算边缘网格的坡度坡向会缺少必需的格网单元值。如图 2(a),以z1和z12分别代表移动窗口在RSG四周和角点边缘时的两种典型情况,z1和z12为待计算单元,空白单元为移动窗口计算时的缺失单元。为了使移动窗口在边缘计算时能够获得足够的网格数,可在DEM格网周围增加一圈虚拟网格。具体步骤为:1) 对边缘格网上下左右各增加一虚拟行或一虚拟列;2) 对虚拟行、列上的在原始RSG上进行高程推算。
为了保证边缘地形的自然延展趋势,增强边缘格网晕渲处理效果,高程推算采用虚拟网格所在行、所在列或所在对角线上邻近的2个网格的高程值进行线性外推,得到虚拟网格的高程值。如图 2(b)所示,阴影部分为实际网格,空白格网部分为虚拟网格,则虚拟网格za、zb、zc的推算值分别为2×z1-z6、2×z1-z2、2×z1-z5,其他虚拟网格的高程可以类似推算。
依据推算格网单元的位置,可以将虚拟网格分为角点网格(左上、左下、右上、右下)、行网格(上、下)、列网格(左、右)等3大类8小类,然后分别处理。假设扩展后的RSG为一个n×n的矩阵,则通过虚拟网格的行列号很容易确定网格类别、具体所在行列以及参与该网格推算的两个实际网格,从而计算出虚拟网格的高程。
1.3 RSG网格光照强度计算除了地形起伏形态(主要是坡度、坡向),光源的类型、数量以及方位等对地貌晕渲效果也起着决定性作用[13]。通常,地貌晕渲时将光源设置为方向光源来模拟太阳照射效果,有时也会引入多个光源构成综合光照效果以凸显复杂地形地貌形态[4-5]。本文算法引入一个光源模拟太阳照射,并采用晕渲效果较好的相对辐射模型来计算每个网格的光照强度IR。计算公式为:
${I_R} = {G_{{\text{max}}}} \times \left( {{\text{cos}}\left( {{A_f} - {A_s}} \right){\text{sin}}\;{H_f}\;{\text{cos}}\;{H_s} + {\text{cos}}\;{H_f}\;{\text{sin}}\;{H_s}} \right)$ | (5) |
其中:Gmax为最大灰度级,一般取255;Af为网格的坡向,取值范围为[0, 360];As为太阳方位角,取值范围为[0, 360],以正北方向为0°基准,顺时针增加;Hf为网格的坡度,取值范围为[0,90);Hs为太阳高度角,取值范围为[0,90];IR的取值范围在[0, 255],0为最黑,255为最亮。
2 地貌晕渲并行算法当前,许多个人计算机和工作站都有两个或两个以上的内核。为了充分利用硬件资源,提高计算效率,需要对代码进行并行化,以将工作分摊在多个处理器上。微软从.NET Framework 4开始提供了新的运行时、类库以及诊断工具,增强对并行编程的支持。这些功能简化了并行开发,使程序员能够通过固有方法编写高效、细化且可伸缩的并行代码,而不必直接处理线程或线程池等低级操作。
规则格网DEM的数据结构类似于栅格数据,许多有关计算都是以移动窗口方式进行的局部运算,当将规则格网DEM分割后,各个子块仍然可以独立地完成相应的计算,因此,规则格网DEM具有适用并行计算的特性,并行计算模式成为提升格网DEM地形分析算法效率的有效途径,如并行计算模式下的等高线提取[14]、洼地填充[15]、三维晕渲[16]等。并行化算法在算法设计过程中需要将整个任务划分成若干个互相独立子任务,这些子任务不考虑互相的依赖和顺序。尽管规则格网DEM的晕渲过程依次包括了坡度计算、坡向计算、网格光照强度计算等3个子过程,但每个网格的坡度、坡向及光强度计算仅与其相邻网格有关,因此,地貌晕渲算法的并行化采用数据并行的策略。
规则格网DEM地貌晕渲算法并行化实现的主要步骤是:1) 获取运行计算机的CPU核数与内存,作为数据分块的依据;2) 将规则格网DEM数据根据内核数按2n进行分块,n从0~4;3) 对每个数据分块创建一个线程,并将该线程分配到不同的CPU内核上,经过坡度、坡向以及光照强度计算,获得分块数据的晕渲图像;4) 主线程接收来自各个子线程的晕渲子图像,依据数据分块时子块在整个DEM数据矩阵的行列号确定位置,对晕渲子图像进行合并,从而得到完整的地貌晕渲图。算法流程如图 3所示。该算法通过自动检测运行计算机的内核数和内存容量,从而自动适应计算机的硬件配置,灵活而充分地利用硬件资源。
常用的DEM数据分块方法有按行分割、按列分割、按格网分割3种[17]:按行分割是将DEM数据按照行数分成2n个数据块;按列分割是将DEM数据按照列数分成2n个数据块;按格网分割是将DEM数据按照行数、列数分成2n个数据块。n取整数且2n小于等于核数。为了使各个线程上的任务能够同时结束,从而保证主线程在图像合并时能得到完整的内容,需要均匀分块,保证各个子块处理时间大致相等。
地貌晕渲算法采用格网均匀分割,分割后每个子块都需要处理边缘网格(见1.2节),因此,对DEM数据进行分割时,在计算好分割线的情况下,可在分割线的基础上多增加一行或一列,即DEM相邻分块之间有2行或2列的重叠区。以16×16的原始格网DEM均匀分割为4个子块为例(图 4),每个子块应该为4×4的格网。为处理边缘,可通过增加内部重叠网格和边缘虚拟网格将每个4×4子块扩展为5×5的格网,从而既消除了因DEM数据分割而产生的边缘网格晕渲误差,又满足了数据子块内各个网格单元的移动窗口计算要求。
.NET Framework4.0通过引入并行库TPL实现了基于任务设计而不用处理重复复杂线程的并行开发框架,支持数据并行、任务并行与流水线。为了简化DEM晕渲算法的并行化过程,采用Parallel类的Invoke方法,它可以实现对给定任务实现并行开发。
算法首先判断计算机的内核数以及内存大小,依据配置确定采用的线程数及任务分解方式。在DEM数据分块后,将不同的数据子块以及按照范围分别载入相应的晕渲方法中,并传入太阳高度角、太阳方位角、画布等参数。利用Invoke方法并行化的过程就是将每个分块对应的晕渲方法作为参数传递给Invoke方法即可。当所有晕渲线程完成晕渲后,将结果图片进行拼接。需要注意的是,拼接过程要等所有子线程都处理完成后再进行,否则导致部分子块晕渲图像丢失。由于对DEM进行了均匀分割,因此每个线程执行的时间几乎相等,这样避免了单个线程运行时间过长而影响整体计算时间。以伪代码方式,并行化晕渲算法如下:
GetCPUNumber(); //获取CPU数
GetMemory(); //获取内存大小
Switch(CPUNumber){
Case 1: //单核
bitmap11=new Bitmap(DEM行数, DEM列数);
Graphics11=Graphics.FromImage(bitmap11);
parallelCalBlock11=new parallelCal(DEM, Graphics11, 太阳方向角, 太阳高度角, 子块起始行号, 子块起始列号, bitmap11高度, bitmap11宽度);
Parallel.Invoke(parallelCalBlock11.DEMHillshade);
MergeImage(); //图片拼接
Case 2: //双核,DEM被横向分割为两块
bitmap11=new Bitmap(DEM行数/2, DEM列数);
Graphics11=Graphics.FromImage(bitmap11);
parallelCalBlock11=new parallelCal(DEM, Graphics11, 太阳方向角, 太阳高度角, 子块起始行号, 子块起始列号, bitmap11高度, bitmap11宽度);
bitmap12=new Bitmap(DEM行数/2, DEM列数);
Graphics12=Graphics.FromImage(bitmap12);
parallelCalBlock12=new parallelCal(DEM, Graphics12, 太阳方向角, 太阳高度角, 子块起始行号, 子块起始列号, bitmap12高度, bitmap12宽度);
Parallel.Invoke(parallelCalBlock11.DEMHillshade, parallelCalBlock12.DEMHillshade);
MergeImage(); //图片拼接
Case 4://四核
…
}
程序运行时,主线程由系统自动生成,负责监管整个程序的运行,统一进行系统资源的调配;任务线程由主线程逐个启动,线程数量取决于数据块的数量,同时被系统分配到不同的CPU内核中,不同内核上运行的任务线程互相独立。
2.3 晕渲结果拼接单机多核计算机作并行计算时各个线程共享内存。内核中任务线程运行的结果保存在内存中,待所有任务线程结束后,需要将不同内核生成的图像按照数据块的对应位置由主线程进行拼接,形成一张图像作为最后晕渲结果。图形设备接口GDI+(Graphics Device Interface Plus)提供了3种实现图像拼接的途径:像素复制、BitBlt函数以及DrawImage函数。像素复制利用BitMap类的GetPixel和SetPixel函数在像素操作层次上实现图像拼接,效率较低,而调用Windows应用编程接口(Application Programming Interface,API)提供的BitBlt函数实现拼接又较为复杂,因此本文采用DrawImage函数。具体步骤是:首先,利用Bitmap类创建一个与DEM原始行列数相等的空白图像;其次,将空白图像作为画布;最后,按照DEM各个子块对应的位置将各个晕渲子图像绘制到该画布上就可以完成晕渲图像的拼接。
3 实验结果与分析为了验证上述地貌晕渲算法的效果和效率,在.NET多核编程环境下,使用C#语言编写了地貌晕渲并行处理算法。计算环境配置为:CPU Intel 2.6 GHz,4个物理硬件核心,内存为4 GB,Windows 7操作系统。
3.1 地貌晕渲效果分析选取地形起伏明显的某山地区域,格网DEM大小为1 001×1 001,设定正东方向为0°方位角。首先,分别测试了本文算法在不同方位角和高度角下的地貌晕渲效果;然后在相同光照条件下本文算法与ArcGIS软件的地貌晕渲效果。
3.1.1 太阳方位角对地貌晕渲效果的影响选取太阳高度角为45°,方位角分别取45°、135°、225°和315°时的地貌晕渲效果如图 5所示。
从图 5中可以看出,太阳方位角对地貌晕渲的立体感影响不大,但对地形的阴影分布和立体形态有着较为明显的影响,方位角为45°、135°时地貌晕渲的立体形态和方位角为225°、315°时的结果甚至呈现反地形现象。
3.1.2 太阳高度角对地貌晕渲效果的影响选取太阳方位角为45°,高度角分别取30°、45°、60°和90°时的地貌晕渲效果,如图 6所示。从图 6中可以看出,太阳高度角对地貌晕渲效果的影响十分明显,当高度角为45°时,地貌晕渲立体效果最好;当高度角为90°时地形立体效果几乎就消失了。
设置太阳方位角和高度角均为45°,本文算法与ArcGIS的晕渲结果对比如图 7所示。两者都对不同网格灰度值之间的过渡较平滑,立体感强。相对于ArcGIS,本文算法的晕渲效果中亮度对比适中,而ArcGIS晕渲的结果则相对偏暗,山谷河谷区域颜色太深,掩盖了部分地形细节。
为了研究CPU物理内核数目与DEM晕渲耗时的关系,使用行列数分别为500×500、1 000×1 000、5 000×5 000、10 000×10 000等四组不同大小的DEM数据,在相同实验计算环境下,对不同核数参与计算进行效率统计。算法使用多核并行化后的效果可用加速比Sp=T1/Tp以及并行效率ηp=Sp/p衡量[18]。p代表核数,T1代表最佳串行计算时间,Tp代表使用p个核并行计算的时间。实际操作时T1也往往表示使用一个核计算的时间。表 1列出了各组DEM数据在不同CPU核数参与计算情况下的耗时以及加速比,图 8为不同数据规模下不同核数参与计算时的算法耗时对比。取表 1中S2的平均值计算η2为90.31%;同理,计算η4为81.84%。
由测试统计结果可知,DEM晕渲并行化大大提高了算法效率,充分发挥了硬件性能。依据Amdahl定律,算法并行化效率的提升与算法的串行化比例有关,因此,尽管算法效率与参与CPU核数呈正相关关系,但并不会出现核数增加一倍,耗时减少一半的理想情况。另外,核数的增加需要更多的数据分块和数据拼接处理时间,CPU向内存频繁存取数据以及协调各个线程的开销增加,算法的并行效率会有所降低。
4 结语目前大多数计算机都支持多核计算,利用多核计算模式来提升传统算法的计算效率是一种可行方法。本文利用多核编程技术对现有地貌晕渲算法进行了改进和实验验证。地貌晕渲效果实验表明:1) 本文算法具有良好晕渲效果;2) 光源设置与地貌晕渲效果之间关系密切。效率实验结果表明:参与计算核数与地貌晕渲效率大致呈线性正相关关系,多核并行计算可以大大提高地貌晕渲算法的效率。另外,本文的地貌晕渲并行化算法还不能处理大规模格网DEM数据,若DEM数据规模大于计算机内存容量,则需要开辟专门的线程进行内存与硬盘数据之间的数据调度,这也是下一步要研究的内容。
[1] | 张玲, 杨晓平, 魏占玉, 等. 三维数据的二维可视化方法综述[J]. 地震地质, 2014, 36(1): 275-284. (ZHANG L, YANG X P, WEI Z Y, et al. Overview of visualization methods of three dimensional topographic data[J]. Seismology and Geology, 2014, 36(1): 275-284.) |
[2] | 王少荣, 陈星宇, 徐秋红. 基于曲率的小比例尺地貌晕渲图生成方法[J]. 计算机应用与软件, 2015, 32(7): 204-206. (WANG S R, CHEN X Y, XU Q H. Generation method of small-scale relief shading based on curvature[J]. Computer Applications and Software, 2015, 32(7): 204-206.) |
[3] | 李少梅, 阚映红, 陈艳丽, 等. 数字地貌晕渲光照模型研究[J]. 测绘科学技术学报, 2010, 27(1): 57-60. (LI S M, KAN Y H, CHEN Y L, et al. Research on illumination model of digital relief shading[J]. Journal of Geomatics Science and Technology, 2010, 27(1): 57-60.) |
[4] | 丁宇萍, 蒋球伟. 地貌晕渲图的生成原理与实现[J]. 计算机应用与软件, 2011, 28(9): 214-216. (DING Y P, JIANG Q W. Hill shading map generation principles and realization[J]. Computer Applications and Software, 2011, 28(9): 214-216.) |
[5] | VERONESI F, HURNI L. A GIS tool to increase the visual quality of relief shading by automatically changing the light direction[J]. Computers & Geosciences, 2015, 74: 121-127. |
[6] | JIN H, JESPERSEN D, MEHROTRA P, et al. High performance computing using MPI and OpenMP on multi-core parallel systems[J]. Parallel Computing, 2011, 37(9): 562-575. DOI:10.1016/j.parco.2011.02.002 |
[7] | GEPNER P, KOWALIK M F. Multi-core processors:new way to achieve high system performance[C]//PARELEC 2006:Proceedings of the 2006 International Symposium on Parallel Computing in Electrical Engineering. Washington, DC:IEEE Computer Society, 2006:9-13. |
[8] | ZHANG S H, XIA Z X, YUAN R, et al. Parallel computation of a dam-break flow model using OpenMP on a multi-core computer[J]. Journal of Hydrology, 2014, 512: 126-133. DOI:10.1016/j.jhydrol.2014.02.035 |
[9] | ZHOU G Y, LIU X L, FU S H, et al. Parallel identification and filling of depressions in raster digital elevation models[J/OL]. International Journal of Geographical Information Science, 2016:1-18[2017-01-10]. http://dx.doi.org/10.1080/13658816.2016.1262954. |
[10] | YILDIRIM A A, WATSON D, TARBOTON D, et al. A virtual tile approach to raster-based calculations of large digital elevation models in a shared-memory system[J]. Computers & Geosciences, 2015, 82: 78-88. |
[11] | SHARMA R, XU Z W, SUGUMARAN R, et al. Parallel landscape driven data reduction & spatial interpolation algorithm for big LiDAR data[J/OL]. ISPRS International Journal of Geo-Information, 2016, 5(6):97:1-14[2016-10-15]. http://dx.doi.org/10.3390/ijgi5060097. |
[12] | YANG J H, ZHANG J X. Parallel performance of typical algorithms in remote sensing-based mapping on a multi-core computer[J]. Photogrammetric Engineering & Remote Sensing, 2015, 81(5): 373-385. |
[13] | 陈艳丽, 陈欣, 魏兰. 数字地貌晕渲光源的设置和应用研究[J]. 测绘与空间地理信息, 2016, 39(1): 91-93. (CHEN Y L, CHEN X, WEI L. The research on setting and application of light sources of digital relief shading[J]. Geomatics & Spatial information Technology, 2016, 39(1): 91-93.) |
[14] | 王宗跃, 马洪超, 徐宏根, 等. 在集群多核CPU环境下的等高线并行提取方法[J]. 计算机工程与应用, 2010, 46(17): 5-7. (WANG Z Y, MA H C, XU H G, et al. Method of parallel contours extraction based on multi-core CPU[J]. Computer Engineering and Application, 2010, 46(17): 5-7. DOI:10.3778/j.issn.1002-8331.2010.17.002) |
[15] | 徐昕. 基于单机多核格网DEM填洼并行算法的研究与实现[J]. 城市地理, 2015(24): 96-97. (XU X. Research and implementation of parallel algorithm filling grid DEM based on multi-core computer[J]. Cultural Geography, 2015(24): 96-97. DOI:10.3969/j.issn.1674-2508.2015.24.073) |
[16] | 陈景广, 佘江峰, 宋晓群, 等. 基于多核CPU的大规模DEM并行三维晕渲[J]. 武汉大学学报(信息科学版), 2013, 38(5): 618-621. (CHEN J G, SHE J F, SONG X Q, et al. Parallel rendering of large-scale DEM based on multi-core CPU[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5): 618-621.) |
[17] | 黄先锋, 孙岩标, 张帆, 等. 多核计算环境下的LiDAR数据DEM内插方法研究[J]. 山东科技大学学报(自然科学版), 2011, 30(1): 1-6. (HUANG X F, SUN Y B, ZHANG F, et al. Research on DEM interpolation method with Lidar data under multi-core computing environment[J]. Journal of Shandong University of Science and Technology (Natural Science), 2011, 30(1): 1-6.) |
[18] | 王顺绪. 多核计算机上并行计算的实现与分析[J]. 淮海工学院学报(自然科学版), 2009, 18(3): 30-33. (WANG S X. Implementation and analyses of parallel computation on multi-core computers[J]. Journal of Huaihai Institute of Technology (Natural Science Edition), 2009, 18(3): 30-33.) |