一种基于CUDA的DEM特征点提取并行化方法

摘要

本发明公开了一种基于CUDA的DEM特征点提取并行化方法，首先对目前比较常用于DEM地形特征点提取的maximum z‑tolerance方法进行优化，得到改进maximum z‑tolerance方法，之后利用CUDA对改进的maximum z‑tolerance方法进行并行化处理，从而得到高精度高效率的DEM特征点提取并行化方法。本发明不仅完成了maximum z‑tolerance方法的优化，而且将并行技术CUDA应用于数字地形分析领域，在保证DEM特征点提取算法的精度情况下，提高了DEM地形特征点的提取效率，保证了多尺度DEM的生成效率，满足了实际应用中的尺度需求。

权利要求

1.一种基于CUDA的DEM特征点提取并行化方法，其特征在于，包括以下步骤：步骤1：对maximum z-tolerance方法进行优化，得到改进的maximum z-tolerance方法；

所述对maximum z-tolerance方法进行优化，是通过在每次迭代过程中，将寻找TIN的所有三角面中具有最大高程差的DEM像元点优化为寻找TIN中每个三角面中具有最大高程差的所有DEM像元点，以减少迭代次数来优化maximum z-tolerance方法，从而得到改进的maximum z-tolerance方法；

步骤2：对改进的maximum z-tolerance方法进行并行化处理，得到DEM特征点提取并行化方法；

步骤2的具体实现包括以下步骤：

步骤2.1：对步骤1中改进的maximum z-tolerance方法的串行程序进行程序性能分析，得到其中计算量最大耗时最多且具有可并行性的函数；

步骤2.2：利用CUDA编程模型和CUDA存储模型分析并行过程中数据的传输过程；

步骤2.3：根据CUDA并行模式选择适宜的线程划分方式；

具体实现过程是，当计算每个三角面内具有最大高程差的特征点时，由于TIN中待计算的三角面均存储在线性结构中，并且每个线程负责处理一个三角面，所以线程划分方式是：threads<<>>,blocks<<<(triangleNumber+TC*TC-1)/(TC*TC),1>>>；

其中TC是每个线程块中的x方向和y方向上的线程数量，triangleNumber是TIN中三角面的数量；

步骤2.4：将步骤2.1中分析得到的函数改进成CUDA中的核函数，从而得到DEM特征点提取并行化方法；

步骤3：根据步骤2中的DEM特征点提取并行化方法，构建不规则三角网，从而生成数字地形模型。

2.根据权利要求1所述的基于CUDA的DEM特征点提取并行化方法，其特征在于：DEM特征点提取并行化方法的效率和精度评估，包括以下子步骤：步骤4.1：评估DEM特征点提取并行化方法的效率，包括maximum z-tolerance方法与DEM特征点提取并行化方法运算时间的对比，计算出加速比；

步骤4.2：评估生成的数字地形模型的精度，包括将步骤3中生成的数字地形模型与国家基础地理信息中心现有的数字地形模型进行比较和评估地形参数和地形特征。

3.根据权利要求2所述的基于CUDA的DEM特征点提取并行化方法，其特征在于：所述将步骤3中生成的数字地形模型与国家基础地理信息中心现有的数字地形模型进行比较，计算出均方根误差。

4.根据权利要求2所述的基于CUDA的DEM特征点提取并行化方法，其特征在于：所述地形参数包括平均坡度和地表粗糙度，所述主要地形特征包括流域匹配度和流域线匹配误差。

说明书

一种基于CUDA的DEM特征点提取并行化方法

技术领域

[0001] 本发明属于数字地形分析技术领域，特别涉及一种基于CUDA的DEM特征点提取并行化方法。

背景技术

[0002] 目前数字高程模型(DEM)综合方法主要包括基于栅格DEM的综合方法(如重采样，地形特征保持方法)和基于矢量TIN的综合方法(如提取特征点、线方法)。其中较为常用的重采样方法具有“削峰填谷”的平滑作用，地形参数随着尺度的变化会发生较大变化，无法保持地形特征一致性；地形特征保持方法在制图领域基于等高线进行了精度比较，在地形分析领域有待验证；基于不规则三角网(TIN)的DEM综合方法是较为常用且精度较高的DEM综合方法，该综合方法通常采用DEM特征点和特征线来构建TIN，生成多尺度DEM，而DEM特征点提取是其最为重要的部分。

[0003] 现有比较经典的DEM特征点提取方法有道格拉斯-普克算法、VIP算法、Heuristic算法和maximum z-tolerance算法。道格拉斯-普克算法(Douglas and Peucker,1973,参考背景文献1)是将曲线近似为一系列的点，通过设定阈值，减少曲线上点的数量从而来获取重要地形特征点。VIP(Very Important Point)算法是在利用3x3局部窗口进行平移的过程中，利用中心点周围的8个像元的高程值内插出中心点的高程值，最后通过内插值与实际值偏差来衡量点的重要性，重要性高的点即为地形特征点(Chen and Guevara,1987,参考背景文献2)。Heuristic算法(Lee,1991,参考背景文献3)是通过设定固定的特征点数量，逐渐删除真实值与估算值差值最小的特征点，直到保留设定的特征点数量，该算法又叫启发式丢弃算法，通过与阈值的比较逐步删除点来构建TIN拟合地形表面。Hierarchy算法(DeHaemer and Zyda,1991,参考背景文献4)首先构建覆盖区域的少数三角形，然后采用四叉树、k-d树等方式拆分三角形构建树型层次结构，具有最大误差值的点作为拆分点，即地形特征点。Maximum z-tolerance算法(Chang,2007,参考背景文献5)采用一种迭代添加点的方法，通过比较点与所在三角面的高程差与阈值的大小来提取地形特征点，如果大于阈值，便将该点标记为特征点，重新构建TIN，再进行下一次的特征点提取，直到所有的格网点到对应三角面的高程差均在指定阈值内。3D Douglas-Peucker算法(Fei等,2009,参考背景文献6)将道格拉斯-普克算法从二维扩展到三维领域用于特征点提取。

[0004] 研究表明maximum z-tolerance方法应用最为广泛且精度高，但其计算速度慢。因此，本发明提出对maximum z-tolerance方法的优化，得到改进的maximum z-tolerance方法，并将并行计算嵌入改进的maximum z-tolerance方法进行并行化处理，得到DEM特征点提取并行方法，从而在保证计算精度的前提下，提高DEM地形特征点的提取效率。相对于CPU的多核并行，GPU并行计算的性能更好，适合于并行处理大数据量、高并行性、低数据耦合、高计算密度、固定化流程的问题。

[0005] CUDA(统一计算设备架构)是2006年NVIDIA公司研发的基于GPU进行并行加速的开发环境和软件体系，通过创建大量的线程来计算大量具有重复计算的算法以提高可并行性算法的计算效率。自推出之后，CUDA就已被广泛应用于矩阵运算、水文分析、遥感图像增强等应用领域。万单领(万单领，2007,参考背景文献7)利用CUDA对Delaunay三角剖分中三角网如何生成和分割进行并行处理，研究表明该种并行算法不仅能得到较准确的实验结果且计算速度较快。Rueda(Rueda and Ortega，2008，参考背景文献8)等采用CUDA进行区域内点的包含和自交检测，研究表明加速比能够提高到36.0。Micikevicius(Micikevicius，2009，参考背景文献9)基于CUDA完成了三维有限差分计算的并行化改进，表明并行模式的不同会影响差分计算的结果以及计算效率。Ortega(Ortega and Rueda，2010，参考背景文献10)等利用CUDA进行河流网提取算法(D8算法)进行并行处理，研究表明河流网提取的加速比能够高达8.0。徐赛华(徐赛华，2011，参考背景文献11)等基于CUDA进行大规模三维数据的并行可视化，并行将三维数据转换成图形或图像，显示在屏幕上，在图像显示质量不变的情况下，提高了计算效率。Qin(Qin and Zhan，2012，参考背景文献12)等利用基于CUDA对水流量的计算进行并行处理，分别利用D8算法和基于图形理论的MFD算法来计算水流量，表明前者的加速比能够达到22.3，后者的加速比能够达到10.9。Mielikainen(Mielikainen等，2013，参考背景文献13)等利用CUDA对WRF模型进行加速改进得到高效率的WRF模型，基于CUDA的WRF模型加速比能够达到70.0，当利用4个GPU时，加速比能够达到132.0。Wu(Wu等，2015，参考背景文献14)等基于CUDA实现了三维地形实时渲染，研究表明并行算法能够很大程度地提高计算效率，加速比能够达到212.3。Rueda(Rueda等，2016，参考背景文献15)等利用CUDA对D8算法提取河流网进行并行化改进，并将其与OpenACC并行进行对比分析，发现基于CUDA的D8并行算法，提取效果很好且代码简单，计算效率有很多提升。

[0006] 通过对上述内容的了解，发现还没有任何将CUDA用于格网DEM生成TIN的DEM特征点提取方法。因此，本发明旨在改进的maximum z-tolerance方法的基础上，利用CUDA对其进行并行化，最后达到在保证精度的前提下快速提取DEM特征点的目的。

[0007] 文献1.Douglas D H,Peucker T K.1973.Algorithms for the reduction of the number of points required to represent a digitized line or its caricature[J].Canadian Cartographer,10(2):112–122；

[0008] 文献2.Chen Z,Guevara J.1987.Systematic selection of very important points(VIP)from digital terrain model for constructing triangular irregular networks[C].Proc.8th International Symposium on Computer-Assisted Cartography(AutoCarto 8),Baltimore,Maryland,29:50–56；

[0009] 文献3.Lee,J.1991.Comparison of existing methods for  building triangular irregular network models of terrain from grid digital elevation models[J].International Journal of Geographical Information Systems,5(3):267–285；

[0010] 文献4.DeHaemer Jr M,Zyda M J.1991.Simplification of objects rendered by polygonal approximations[J].Computers and Graphics,15(2):175–184；

[0011] 文献5.Chang K.2007.Introduction to Geographic Information Systems[M].4th edition.New York:McGraw-Hill；

[0012] 文献6.Fei L,He J.2009.A three-dimensional Douglas-Peucker algorithm and its application to automated generalization of DEMs[J].International Journal of Geographical Information Science,23(6):703–718；

[0013] 文献7.万单领.2007.基于GPU加速的并行粒子群算法及其应用[D].大连:大连理工大学；

[0014] 文献8.Rueda A J,Ortega L.2008.Geometric algorithms on CUDA[J].In:Proceedings of the International Conference on Computer Graphics,Theory and Applications,Funchal,Madeira,Portugal,107–112；

[0015] 文献9.Micikevicius P.2009.3D Finite Difference Computation on GPUs using CUDA[J].Proc.of 2nd GPGPU2,Mar:79-84；

[0016] 文献10.Ortega L,Rueda A J.2010.Parallel drainage network computation on CUDA[J].Computers&Geosciences,36:171-178；

[0017] 文献11.徐赛华,张二华.2011.基于CUDA的三维数据并行可视化[J].CCT理论与应用研究,20(1):47-54；

[0018] 文献12.Qin  Chengzhi,Zhan  Lijun.2012.Parallelizingflow-accumulationcalculationsongraphicsprocessing  units—from 
iterativeDEMpreprocessingalgorithmtorecursive multiple-flow-
directionalgorithm[J].Computers&Geosciences,43:7-16；

[0019] 文献13.Mielikainen Jarno,HuangBormin,WangJun,HuangHLA,Goldberg Mitchell D.2013.Compute unified device architecture(CUDA)-based parallelization of WRF Kessler cloud microphysics scheme[J].Computers&Geosciences,52:292-299；

[0020] 文献14.Wu Jiaji,Deng Long,Anand Paul.2015.3D Terrain Real-time Rendering Method Based on CUDA-OpenGL Interoperability[J].IETE Technical Review,32(6):471-478；

[0021] 文献15.RuedaAntonioJ,Jose  MNoguera,AdrianLuque.2016.A comparisonofnativeGPUcomputingversusOpenACCfor implementing flow-
routingalgorithmsinhydrologicalapplications[J].Computers&Geosciences,87:91-
100.

发明内容

[0022] 为了解决现有的DEM特征点提取方法无法满足日益增长的数据处理需求问题，本发明提供了一种能够高精度高效率地提取DEM特征点的并行化方法。

[0023] 本发明所采用的技术方案是：一种基于CUDA的DEM特征点提取并行化方法，其特征在于，包括以下步骤：

[0024] 步骤1：优化目前比较常用于DEM特征点提取的maximum z-tolerance方法，得到改进的maximum z-tolerance方法；

[0025] 步骤2：对改进的maximum z-tolerance方法进行并行化处理，得到DEM特征点提取并行化方法；

[0026] 步骤3：根据步骤2中的DEM特征点提取并行化方法，构建不规则三角网，从而生成数字地形模型。

[0027] 作为优选，步骤1中所述对maximum z-tolerance方法的优化，是通过在每次迭代过程中，将寻找TIN的所有三角面中具有最大高程差的DEM像元点优化为寻找TIN中每个三角面中具有最大高程差的所有DEM像元点，以减少迭代次数来优化maximum z-tolerance方法，从而得到改进的maximum z-tolerance方法。

[0028] 作为优选，步骤2的具体实现包括以下步骤：

[0029] 步骤2.1：使用程序性能分析工具对步骤1中改进的maximum z-tolerance方法的串行程序进行程序性能分析，得到其中计算量最大耗时最多且具有可并行性的函数；

[0030] 步骤2.2：利用CUDA编程模型和CUDA存储模型分析并行过程中数据的传输过程；

[0031] 步骤2.3：根据CUDA并行模式选择适宜的线程划分方式；

[0032] 步骤2.4：将步骤2.1中分析得到的函数改进成CUDA中的核函数，从而得到高精度高效率的DEM特征点提取并行化方法。

[0033] 作为优选，步骤2.3的具体实现过程是，当计算每个三角面内具有最大高程差的特征点时，由于TIN中待计算的三角面均存储在线性结构中，并且每个线程负责处理一个三角面，所以其线程划分方式是threads<<>>,blocks<<<(triangleNumber+TC*TC-1)/(TC*TC),1>>>,其中TC是每个线程块中的x方向和y方向上的线程数量，triangleNumber是TIN中三角面的数量。

[0034] 作为优选，DEM特征点提取并行化方法的效率和精度评估，包括以下子步骤：

[0035] 步骤4.1：评估DEM特征点提取并行化方法的效率，包括maximum z-tolerance方法与DEM特征点提取并行化方法运算时间的对比，计算出加速比；

[0036] 步骤4.2：评估生成的数字地形模型的精度，包括将步骤3中生成的数字地形模型与国家基础地理信息中心现有的数字地形模型进行比较和评估地形参数和主要地形特征。

[0037] 作为优选，所述将步骤3中生成的数字地形模型与国家基础地理信息中心现有的数字地形模型进行比较，计算出均方根误差。

[0038] 作为优选，所述地形参数包括平均坡度和地表粗糙度，所述主要地形特征包括流域匹配度和流域线匹配误差。

[0039] 依照本发明所提供基于CUDA的DEM特征点提取并行化方法，可以高精度高效率地进行DEM地形特征点的提取,从而解决现有算法无法满足日益增长的数据处理需求问题，以快速生成DEM。因此本发明特别应用于数字地形分析领域，快速实现DEM的地形建模。

附图说明

[0040] 图1为本发明实施例的流程图；

[0041] 图2为本发明实施例的改进的maximum z-tolerance方法进行并行化处理流程图；

[0042] 图3为本发明实施例的生成数字地形模型流程图；

[0043] 图4为本发明实施例的DEM特征点提取并行化方法的效率和精度评估流程图。

[0044] 具体实施方法

[0045] 为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

[0046] 本发明要解决的核心问题是：得到一个高效率高精度的DEM特征点提取并行化方法，根据多尺度DEM的应用需求从精细尺度的DEM数据中快速提取DEM特征点，从而快速地完成DEM的地形建模。

[0047] 请见图1，本发明提供的一种基于CUDA的DEM特征点提取并行化方法，包括以下步骤：

[0048] 步骤1：优化目前比较常用于DEM特征点提取的maximum z-tolerance方法，得到改进的maximum z-tolerance方法；

[0049] 为了便于实施，本发明中现有的DEM特征点提取方法采用的是应用最为广泛且精度高但计算速度慢的maximum z-tolerance方法。首先，在此介绍maximum z-tolerance方法的思想：首先将DEM的四个角点设为有效值，连接任意一条对角线，得到初始的不规则三角网TIN并使其覆盖整个DEM，然后遍历TIN中每个三角面，计算其所包含的所有像元点在三角面上的高程值与像元点原始高程值的差值，确定高程差值最大的点，并将最大高程差值与阈值进行比较。如果最大高程差值小于选定的阈值，结束DEM特征点的提取，用已经提取的DEM特征点构建TIN即可。反之，如果最大高程差值大于选定的阈值，则将该点标识为DEM特征点并将其添加到原来的DEM特征点集。当现存TIN中的每个三角面都被遍历后，重新构建TIN，持续执行以上迭代过程直到最大高程差值均在阈值范围内，得到最终的DEM特征点。由此可知，maximum z-tolerance方法在每次迭代过程中，只在TIN的所有三角面中找到一个具有最大高程差的点。

[0050] 具体实施时，通过在每次迭代过程中，寻找TIN的所有三角面中具有最大高程差的点，以减少迭代次数来优化maximum z-tolerance方法，从而得到改进的maximum z-tolerance方法。本发明对此选择了陕西韭园沟地区的DEM数据为实验数据，其空间分辨率为5m，像元大小为2112x1792，数据范围为东经110°15′-110°22′，北纬37°32′-37°37′，并选取了阈值分别为0.1，0.5，1.0，1.5，2.0，2.5和3.0。通过实验验证，表明改进的maximum z-tolerance算法的加速比分别是22.29，17.52，13.02，12.04，12.21，10.73和8.89，并且该算法的精度还略有提高。

[0051] 步骤2：对改进的maximum z-tolerance方法进行并行化处理；

[0052] 请见图2，步骤2的具体实现包括以下步骤：

[0053] 步骤2.1：使用程序性能分析工具对步骤1中改进的maximum z-tolerance方法的串行程序进行程序性能分析，得到其中计算量最大耗时最多且具有可并行性的函数；

[0054] 具体实施时，可以使用Intel@V TuneTM性能分析器对步骤1中的改进的maximum z-tolerance方法的串行程序进行了程序性能分析，发现其中计算每个三角面内的最大高程差是其计算量最大耗时最多且具有可并行性的函数。

[0055] 步骤2.2：利用CUDA编程模型和CUDA存储模型分析并行过程中数据的传输过程；

[0056] 研究CUDA编程模型和CUDA存储模型，可知CUDA编程模型中最主要的运行设备是主机(host)和设备(device)；其通常用CPU作为host(主要是因为GPU分支控制能力低下)，其主要任务是设备和数据的初始化；将GPU这种能执行多个线程的计算设备看作device，让GPU来并行处理处理大量的计算密集任务。另外，一个系统中可以存在一个host，但可以存在多个device。CUDA可以将函数定义为内核函数(kernel)，kernel的计算单元包括线程、线程块和网格，其中若干个线程(thread)构成一个线程块(block)，若干个线程块组成一个网格(grid)。核函数运行时，若干个线程并行执行核函数，以提高运行效率。kernel需要使用_global_或_device_来定义，以说明此函数在GPU上运行，同时使用<<>>来定义表明每次最多有x*y线程(受GPU的计算能力限制)同时进行核函数计算。执行内核的每个线程都有一个独一无二属于自己的线程号，该线程号可以通过CUDA的内置变量在kernel中对其进行访问和查询。CUDA中的存储器种类有：全局存储器(global memory)、寄存器(Registers)、共享存储器(shared memory)、本地存储器(local memory)、常数存储器(constant memory)和纹理存储器(texture memory)。CUDA中每个grid中的所有线程都能访问同块全局存储器，每个block都拥有自己的共享存储器，每个thread都有自己独立的寄存器和本地存储器，所有线程都能够以只读形式访问常数存储器和纹理存储器。另外，不同的存储器具有不同的访问速度。

[0057] 具体实施时，根据CUDA编程模型和CUDA存储模型分析并行过程中数据的传输过程。该过程是首先在CPU中利用DEM数据的四个角点构建TIN，然后将DEM数据、TIN中的三角面坐标集合和指定的阈值通过CUDA中的传递函数cudaMemcpy()拷贝到GPU的global memory，最后TIN中所有的三角面并行进行核函数计算，将计算结果通过CUDA的传递函数cudaMemcpy()传回CPU。其中，利用寄存器存储访问速度的优势，将计算中无需传出的过程数据临时存储于寄存器。如果传到CPU的结果不为空，则将传出的DEM特征点插入原始的DEM特征点集合。利用新DEM特征点集合重新构建TIN，并重新向GPU中传入最新TIN中的三角面坐标集合，再持续执行以上过程直到核函数运行结果为空，即可得到最终的DEM特征点。

[0058] 步骤2.3：根据CUDA并行模式选择适宜的线程划分方式；

[0059] 研究CUDA的两层并行模型，包括网格(grid)中线程块(block)间的并行和线程块(block)中线程(thread)间的并行。CUDA使得GPU能够实现线程间的消息通信，在同一个线程块中的线程能够在并行执行的同时还能够通过共享存储器实现线程块内线程间的通信。详细地划分方式可参考Shane Cook.2014.CUDA programming A Developer’s Guide to Parallel Computing with GPUs[M].China Machine Press。

[0060] 具体实施时，CUDA中的核函数调用方式为kernel<<>>,其中blocks和threads分别为满足计算需要启动的线程块和每个线程块中的线程数量。当计算每个三角面内具有最大高程差的特征点时，由于TIN中待计算的三角面均存储在线性结构中，并且每个线程负责处理一个三角面，所以其线程划分方式是threads<<>>,blocks<<<(triangleNumber+TC*TC-1)/(TC*TC),1>>>,其中TC是每个线程块中的x方向和y方向上的线程数量(受GPU计算能力的限制)，triangleNumber指TIN中三角面的数量(通过对TIN进行遍历得到)。

[0061] 步骤2.4：将步骤2.1中分析得到的函数改进成CUDA中的核函数，从而得到高精度高效率的DEM特征点提取并行化方法。

[0062] 具体实施时，将计算三角面内的最大高程差函数改造成CUDA中的kernel。计算最大高程差kernel的思想：首先根据步骤2.3的线程划分方式，利用CUDA内置变量定义两个线程索引来访问TIN中三角面。例如：idx＝(blockIdx.x*blockDim.x)+threadIdx.x；idy＝(blockIdx.y*blockDim.y)+threadIdx.y；tid＝gridDim.x*blockDim.x*idx+idy。其中threadIdx,blockIdx和blockDim分别是用来记录当前线程的位置索引信息，当前线程块的位置索引信息和定义的线程块位置索引信息。然后，根据定义的线程索引号tid获取待计算的三角面三点的坐标，并根据三点坐标进行平面方程的求解，记为z＝ax+by+c。为了减少传输并加快计算速度，三角面的平面方程计算也在GPU段进行并临时存储于寄存器中。a,b和c可以通过三角面三点的坐标计算得到。接下来，判断位于被计算的三角面内的DEM像元点，判断方法是面积判断法。根据三角面的外接矩形范围，可缩小寻找范围。利用海伦公式计算原始三角面和DEM像元点到三角面三个顶点重新构成的三个三角形的面积，分别记为S0，S1,S2,S3。如果S0和(S1+S2+S3)相等，则认为DEM像元点位于三角面内，否则DEM像元点不在三角面内。再次，计算在三角面内的DEM像元点到三角面平面的最大高程差值，而最大高程差值根据像元点的行列号跟平面方程进行求得，并将最大高程差值与指定的阈值进行比较。如果大于阈值，则将该像元点标记为DEM特征点，并将该DEM特征点的行列号分别存储在数组中。否则不进行标记。最后，将存储行列号的两个数组传到CPU，即可得到所有DEM特征点。到此，完成了改进的maximum z-tolerance方法的并行化，得到了DEM特征点提取并行化方法。

[0063] 步骤3：根据步骤2中的DEM特征点提取并行化方法，构建不规则三角网，从而生成数字地形模型；

[0064] 参考附图3，具体实施时，基于河流网约束TIN的方法可参见Zhou Q,Chen Y,2011.Generalization of DEM for terrain analysis using a compound method.ISPRS Journal of Photogrammetry and Remote Sensing.Vol 66,1.38-45。

[0065] 根据这种CPE方法，实施例利用DEM特征点提取并行化方法提取DEM特征点。为了突出关键流域特征，再将由D8算法提取的河流网嵌入DEM特征点中，最后构建基于河流网约束的TIN，从而生成DEM。

[0066] 本实施例的DEM特征点提取并行化方法源代码是预定的开发环境下的源代码,在该开发环境下实现编译生成高精度高效率的DEM特征点提取方法的动态链接库。其中预定的开发环境为Visual Studio和CUDA相结合的开发环境。

[0067] 请见图4，本实施例针对DEM特征点提取并行化方法的效率和精度进行评估，具体步骤如下：

[0068] 步骤4.1：评估DEM特征点提取并行化方法的效率，包括maximum z-tolerance方法与DEM特征点提取并行化方法运算时间的对比，计算出加速比；

[0069] 在实例中，加速比S的公式如下：

[0070]

[0071] 其中，Tp是maximum z-tolerance方法的运行时间，Ts是DEM特征点提取并行化方法的运行时间，该值越大表示效率越高。

[0072] 步骤4.2：评估生成的数字地形模型的精度，包括将步骤3生成的数字地形模型与国家基础地理信息中心现有的数字地形模型进行比较(计算出均方根误差)和评估地形参数和主要地形特征，所述地形参数包括平均坡度和地表粗糙度，所述主要地形特征包括流域匹配度和流域线匹配误差。

[0073] 在具体实施时，精度评估包括均方根误差的计算和地形参数与主要地形特征的评估，具体实施过程如下：

[0074] 在实例中，均方根误差的计算公式如下：

[0075]

[0076] 其中，Z′是原始高程值，Z是预测高程值，n表示像元总数，下标i表示第i个单元，该值越小表示精度越高。

[0077] 地形参数的评估主要计算平均坡度 和平均地表粗糙度K的公式如下：

[0078]

[0079]

[0080] 其中，S表示坡度，A表示投影面积，n表示单元总数，下标i表示第i个单元，sec表示正割三角函数。

[0081] 主要地形特征的评估主要包括流域匹配度SMR和流域线匹配误差SME。为了计算SMR，首先对从原始DEM提取出的河网生成缓冲区，然后将生成的DEM提取的河流网和该缓冲区进行叠置分析，最后得到SMR，该值越大表示精度越高。计算公式为： 其中，L'表示生成后的DEM提取出的河流网落在缓冲区内的长度，L表示河流网的总长度。

[0082] 为了计算SME，需要计算原始DEM和生成后的DEM提取出的河流网所交叉的面积，该值越小表示精度越高。计算公式为： 其中，ΔA表示原始DEM和生成后的DEM提取出的河流网所交叉的面积，L表示河流网的总长度。

[0083] 应当理解的是，本说明书未详细阐述的部分均属于现有技术。

[0084] 应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。