本发明涉及地理空间数据并行处理技术领域，具体涉及一种基于分布式系统的海量dem金字塔切片并行构建方法。
背景技术：
目前，随着dem（数字高程模型digitalelevationmodel）的需求日益增大，高精度、高分辨率dem数据量呈几何级的增长，dem金字塔切片的理论与方法也趋于成熟，但如何对海量的dem数据高效切片处理，是dem广泛应用亟待解决的关键问题。随着并行计算、分布式计算平台的快速发展和趋于成熟，为海量dem金字塔切片的并行计算提供高效的算力。结合并行计算和成熟的dem采样算法、基于四叉树的分割算法实现dem金字塔切片的并行构建方法，可提高dem数据计算的效率。
近年来，国内外学者开始重点研究海量空间数据的分布式计算和分析，在构建dem金字塔方面，一般还是集中在线程并行策略，主要是结合高性能计算机下利用多线程来构建海量dem金字塔切片，这种方式能够充分发挥机器的性能，但是如果dem数据量较大，一般的计算机内存不能够承受，达不到计算的效果。关于dem分析计算方面也有一些专利。中国专利文献089098709949.5公开了一种面向并行数字地形分析海量dem部署与调度方法，解决了海量dem的地形分析的效率，此专利未涉及dem金字塔切片的内容。中国专利文献cn201710107308.x公开了一种基于dem的内流湖集水区划定方法，其包括获取设定区域内高程数字模型的dem数据，并对其进行切割，得到基础dem数据，根据基础dem数据确立设定区域内的山脊线，并对空间单元进行分割计算空间栅格数据。本专利没有涉及并行和dem金字塔的算法。中国专利文献201410687961.4公布了一种适于规则网格dem的结构化多尺度重构方法，同样没有涉及金字塔建立和并行方法。中国专利文献201410788404.1公开了一种面向关系数据库扩展的自适应影像金字塔切片方法，提出了自适应影像切片模型，根据自适应影像切片模型自适应的选择影像切片的大小。影像切片涉及影像元的空值，不涉及dem切片的拼接方面的方法。在软件方面，主要是商业gis软件具有dem金字塔切片的功能，这些方法的主要采用不同的采样算法来切分dem，未实现分布式的dem金字塔切片等内容。
因此，有必要开发一种基于分布式系统的海量dem金字塔切片并行构建方法，能够对大规模或海量dem的金字塔切片进行并行化切分，能够在高性能计算机集群环境下进行分布式计算，减少金字塔构建的时间，提高计算效率，同时可以结合临近平均值方法计算dem切片边缘的值，以解决相邻切片间的平滑过渡。
技术实现要素：
本发明要解决的技术问题是提供一种基于分布式系统的海量dem金字塔切片并行构建方法，能够对大规模或海量dem的金字塔切片进行并行化切分，能够在高性能计算机集群环境下进行分布式计算，减少金字塔构建的时间，提高计算效率，同时可以结合临近平均值方法计算dem切片边缘的值，以解决相邻切片间的平滑过渡。
为了解决上述技术问题，本发明采用的技术方案是：该基于分布式系统的海量dem金字塔切片并行构建方法，具体包括以下步骤：
s1获取源dem的元数据metadata；
s2建立所述源dem的局部坐标系，对所述源dem进行切分形成若干个可并行计算的切片tile；并根据所述源dem的参数计算每一个所述切片tile的索引值记为key，以<key，tile>键值对为计算单元组成一个保存源dem的元数据的链表list<key，tile>，将<key，tile>按照索引值key分布到不同的计算节点进行计算；
s3根据是否定义了最大层次maxzoom，进行设置或计算最大层次maxzoom，再设置需要的输出切片tile的参数；再根据设定的坐标系和切分规则计算输出切片tile的参数，所述输出切片tile的参数包括切片的长、宽、分辨率等；
s4计算节点并行构建所分配的键值对<key，tile>，根据索引值key合并同一计算节点上的高程值value；分别根据索引值key按照所述输出切片tile的长、宽、分辨率参数计算所对应的切片tile的外包矩形extend，计算输入切片tile中每个网格单元cell的坐标，再根据设定的输出坐标系，计算每个所述网格单元cell在输出坐标系的坐标；根据每个所述网格单元cell的坐标计算相应的所述网格单元cell所在的输出切片tile的索引值key，并将其加入一个映射map中，组成一个映射map<key，list<cell>>，其中key表示切片的索引值，list表示链表，cell表示网格单元；
s5将计算节点上相同索引值key的链表list<cell>合并；根据索引值key合并所有计算节点上的链表list<cell>，形成一个键值对<key，list<cell>>；
s6根据输出切片tile的参数和坐标系，计算输出切片tile的网格单元cell的坐标，接着根据采样方式从合并的键值对<key，list<cell>>中计算网格单元cell的高程值value；
s7利用8邻域算法，修正输出切片tile边缘的高程值；根据生成的所述链表list<key，tile>，循环获取每个切片tile的四边网格单元cell和与之相邻的8个切片tile，取相邻切片tile的网格单元高程值的平均值作为相应网格单元ell的高程值value，达到切片tile拼接的平滑；
s8判断当前层次zoom是否满足终止条件，若满足终止条件，则结束；若不满足终止条件则根据上一级层次计算下一级层次和下一级层次的切片tile的外包矩形extend的宽度和高度，计算方法：nz=cz–1，wh=2(nz 1)，其中nz表示下一级层次，cz表示当前层次，wh表示切片tile的宽和高度；
s9根据上一级层次切片tile和采样算法，并行计算生成下一级层次切片tile和索引值key；根据索引值key、下一级的层次和切片tile的参数计算当前层次切片tile所对应的下一级层次的索引值key，利用所述采样算法从当前层次切片tile采样下一级层次切片的高程值，构成下一级层次切片tile的子切片和索引值key，在每个节点上合并相同的索引值key，构成键值对<key，list<tile>>；接着根据索引值key，合并所有节点上相同的索引值key的切片tile，组成下一级层次的切片tile，再执行步骤s7；
s10循环执行步骤s8~s9，直到判断层次zoom满足终止条件，所述终止条件为：若设定了最小层次minzoom，则判断当前层次zoom是否达到设定最小层次minzoom，若达到，则终止；若未设定最小层次minzoom，若当前层次zoom小于0，则终止。
其中tile表示切片，key表示切片的索引值，value表示高程值，list表示链表，cell表示网格单元，zoom表示层次，map表示映射，extend表示源dem的外包矩形；采用上述技术方案，获取待切分源dem的元数据metadata，作为数据切分和投影转换的基础；建立源dem左上角为原点，向下向右为正方向的坐标系，根据默认的切片tile的宽度和高度对源dem图幅进行切分，根据初步切分的tile的行列号组成索引值key，构建一个以索引值key为键，切片tile为值的键值对；将键值对分布在不同的计算节点进行以下的计算；根据设定的maxzoom或计算的maxzoom，以及设定的输出切片tile的宽高、坐标系和采样算法，每个计算节点对所分配的输入切片tile进行切分计算，输出索引值key和一个包含网格单元cell的链表；接着在每个计算节点上合并相同索引值的网格单元；最后合并不同计算节点上相同的索引值key的cell，输出多个key/value对；根据索引值计算相邻的8个索引值key，获取相邻切片tile与当前切片tile的边缘高程值，将相应的边缘值的平均值作为当前切片tile的边缘值，达到切片tile拼接时的平滑；接着循环计算下一级的层次zoom、切片tile的宽高等数据，以上一级的输出链表list<key，tile>作为输入，计算下一级的链表list<key，tile>；循环终止条件是当前层次zoom小于设定的minzoom或0；其中采样方式为利用采样算法进行采样，采样算法可以是最邻近像元法，双线性内插法，三次立方卷积法等，均为现有技术。
作为本发明的优选技术方案，所述步骤s1中所述元数据metadata包括原始分辨率、坐标系，以及在源坐标系下源dem的外包矩形extent数据。获取源dem的metadata的目的是为了后续计算的需要，这些计算包括将源数据进行初步的切分，以便于后续计算的并行化。
作为本发明的优选技术方案，所述步骤s2中建立源dem的局部坐标系具体为：建立以源dem左上角为原点，向下向右为正方向的坐标系；根据设定切片tile的宽和高，将源dem数据切分为多个切片tile，为了将源dem进行初步的切分，以构成链表list<key，tile>，其中索引值key的计算方法：
其中，c、r表示切片tile的行列号，x、y表示切片tile任意一点的瓦片坐标，xmin表示源dem的外包矩形extent的x最小值，ymin表示源dem的外包矩形extent的y最小值，th是初始切片tile的高，tw是初始切片tile的宽。
作为本发明的优选技术方案，所述步骤s3若未设定层次maxzoom，则需要计算层次maxzoom，计算过程是根据源dem的外包矩形extend的左下点坐标和右上点的坐标，并将左下点坐标和右上点的坐标进行投影转换，转换为经纬度坐标，其中左下点坐标是（slon，slat），右上点的坐标是（elon，elat）；接着根据公式计算层次maxzoom，计算公式如下：
其中，pi表示圆周率，maxzoom为最大层次，dlat表示纬度差，dlon表示经度差，lat1表示slat的弧度，lat2表示elat的弧度，a表示计算的中间值，c表示计算的中间值。计算输出切片tile的extend，需要适当的扩大，目的是为了避免输出切片tile拼接不会出现裂缝等问题。
作为本发明的优选技术方案，所述步骤s4中网格单元cell=<x,y,value>，其中x，y是在切片tile的坐标，value是高程值。
作为本发明的优选技术方案，所述步骤s4并行计算输入切片tile，根据输出切片tile的参数将输入切片tile的dem值划分到不同的输出切片tile中；其中首先根据输入切片tile的索引值key计算tile的外包矩形extend，计算公式：
其中，minx、miny、maxx、maxy分别表示外包矩形extend的最小最大xy坐标值，col、row是输入切片tile的行列号，exminx、exminy表示源dem的外包矩形extend的最小的xy坐标值，h表示切片tile的高，w表示切片tile的宽。
作为本发明的优选技术方案，所述步骤s7利用8邻域算法，修正切片tile边缘的高程值value时，对输出的链表list<key，tile>并行计算相邻的切片tile，取四周除交点外的相邻网格单元cell值，并对上下或左右的cell取值，将取到的cell值取平均值作为指定网格单元cell的高程值value；四个角点的高程值的计算是取切片tile的相邻的另外的三个切片tile与相应角点的值，并取平均值作为相应网格单元cell的高程值value。
与现有技术相比，本发明具有的有益效果为：该基于分布式系统的海量dem金字塔切片并行构建方法具有以下优点：利用高性能计算集群，实现了对大规模dem数据构建dem金字塔切片的并行化，可减少dem金字塔构建的时间，提高计算效率；对切片的边缘高程值进行8邻域平滑处理，以避免dem可视化和分析计算的平滑。
附图说明
图1为本发明的基于分布式系统的海量dem金字塔切片并行构建方法的工作流程图；
图2为本发明的基于分布式系统的海量dem金字塔切片并行构建方法的局部坐标系切分dem示意图；
图3为本发明的基于分布式系统的海量dem金字塔切片并行构建方法的单个计算节点计算示意图；
图4为本发明的基于分布式系统的海量dem金字塔切片并行构建方法的多个计算节点合并示意图；
图5为本发明的基于分布式系统的海量dem金字塔切片并行构建方法边缘平滑示意图；
图6为本发明的基于分布式系统的海量dem金字塔切片并行构建方法可视化图；
其中：图（a）为17级图；图（b）为17、16级镶嵌显示图；图（c）为17级的wireframe图；图（d）为17、16级镶嵌wireframe图。
具体实施方式
下面将结合本发明的实施例图中的附图，对本发明实施例中的技术方案进行清楚、完整的描述。
实施例：如图1~6所示，该基于分布式系统的海量dem金字塔切片并行构建方法，具体包括以下步骤：
其中tile表示切片，key表示切片的索引值，value表示高程值，list表示链表，cell表示网格单元，zoom表示层次，map表示映射，extend表示源dem的外包矩形，<key，tile>表示键值对；
s1使用分辨率30米的源dem数据和一部分高精度源dem数据，数据格式是tif，设定的坐标系为wgs84；分别研发基于flink分布式平台的金字塔切片程序；
获取全国30米源dem的元数据metadata，包括原始的分辨率、坐标系wgs84以及extend；
s2对所述输入源dem数据建立局部坐标系，设定初始划分切片的高和宽均为126（网格），利用滑动窗口方法对源数据进行切分，接着根据tile在局部坐标系的行列号计算key，构建一系列的<key，tile>键值对，按照flink的算法将键值对数据分布到不同的计算节点上进行计算；建立以源dem左上角为原点，向下向右为正方向的坐标系；其中索引值key的计算方法：
其中，c、r表示切片tile的行列号，x、y表示切片tile任意一点的瓦片坐标，xmin表示源dem的外包矩形extent的x最小值，ymin表示源dem的外包矩形extent的y最小值，th是初始切片tile的高，tw是初始切片tile的宽；
s3根据是否定义了最大层次maxzoom，进行设置或计算最大层次maxzoom，再根据设定的坐标系和切分规则计算输出切片tile的参数；所述输出切片tile的参数包括切片的长、宽、分辨率等；若未设定最大层次maxzoom，则需要计算最大层次maxzoom，计算过程是根据源dem的外包矩形extend的左下点坐标和右上点的坐标，并将左下点坐标和右上点的坐标进行投影转换，转换为经纬度坐标，其中左下点坐标是（slon，slat），右上点的坐标是（elon，elat）；接着根据公式计算最大层次maxzoom，计算公式如下：
其中，pi表示圆周率，maxzoom为最大层次，dlat表示纬度差，dlon表示经度差，lat1表示slat的弧度，lat2表示elat的弧度，a表示计算的中间值，c表示计算的中间值；
计算输出切片tile的extend，需要适当的扩大，目的是为了避免输出切片tile拼接不会出现裂缝等问题；所述步骤s3中设置输出切片tile的参数为：宽高均为65（网格）、坐标系为wgs84、索引值key的命名遵循tms标准；
s4计算节点并行构建所分配的键值对<key，tile>，根据索引值key合并同一计算节点上的高程值value；分别根据索引值key按照所述输出切片tile的长、宽、分辨率参数计算所对应的切片tile的外包矩形extend，计算输入切片tile中每个网格单元cell的坐标，对初步切分的tile分别计算；首先根据索引值key计算tile的extend；计算输入切片tile中每个网格单元cell的坐标，其中cell=<x,y,value>，其中x，y是在tile的坐标，value是高程值的坐标，根据设定的输出坐标系wgs84，计算cell在输出坐标系的坐标，根据计算结果计算此cell所在的输出切片tile的索引值key，并将其加入map中，组成一个map<key，list<cell>>；将计算节点上相同索引值key的list<cell>合并，原理见图3所示；所述步骤s4并行计算输入切片tile，根据输出切片tile的参数将输入切片tile的dem值划分到不同的输出切片tile中；其中首先根据输入切片tile的索引值key计算tile的外包矩形extend，计算公式：
其中，minx、miny、maxx、maxy分别表示外包矩形extend的最小最大xy坐标值，col、row是输入切片tile的行列号，exminx、exminy表示源dem的外包矩形extend的最小的xy坐标值，h表示切片tile的高，w表示切片tile的宽；
s5根据索引值key合并所有计算节点上的list<cell>；由于输入切片tile分布在不同的计算节点上，而输出切片tile可能包含不同的输入切片tile，也就是说，输出切片tile的数据可能在不同的计算节点，需要将这些数据拉到同一个节点上进行合并，也就是flink的reduce操作，原理见图4所示。
s6根据输出切片tile的参数和坐标系，计算输出切片tile的网格单元cell的坐标，使用设置的采样方法，从合并的<key，list<cell>>中计算cell的高程值，本实施例中使用的是最邻近采样算法；
s7利用8邻域算法，修正长切片tile边缘的高程值value；根据生成的所述链表list<key，tile>，循环获取每个切片tile的四边网格单元cell和与之相邻的8个切片tile，取相邻切片tile的网格单元高程值的平均值作为相应网格单元ell的高程值value，达到切片tile拼接的平滑；根据输出的键值对<key，tile>，根据索引值key的tms命名规则计算每个切片tile相邻的8个切片tile，取相邻tile的cell值，并作平均值作为相应cell的值，达到tile拼接的平滑；原理示意见图5所示；
s8判断当前层次zoom是否满足终止条件，若满足终止条件，则结束；若不满足终止条件则根据上一级层次计算下一级层次和下一级层次的切片tile的外包矩形extend的宽度和高度，计算方法：nz=cz–1，wh=2(nz 1)，其中nz表示下一级层次，cz表示当前层次，wh表示切片tile的宽和高度；
s9根据上一级层次切片tile和采样算法，并行计算生成下一级层次切片tile和索引值key；根据索引值key、下一级的层次和切片tile的参数计算当前层次切片tile所对应的下一级层次的索引值key，利用所述采样算法从当前层次切片tile采样下一级层次切片的高程值，构成下一级层次切片tile的子切片和索引值key，在每个节点上合并相同的索引值key，构成键值对<key，list<tile>>；接着根据索引值key，合并所有节点上相同的索引值key的切片tile，组成下一级层次的切片tile，再执行步骤s7；本实施例中使用的是最邻近采样算法；
s10循环执行步骤s8~s9，计算下一级的list<key，tile>，直到判断当前层次zoom满足终止条件，所述终止条件为：若设定了最小层次minzoom，则判断当前层次zoom是否达到设定最小层次minzoom，若达到，则终止；若未设定最小层次minzoom，则当前层次zoom小于0，则终止。
以上所述仅为本发明的较佳实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
1.一种基于分布式系统的海量dem金字塔切片并行构建方法，其特征在于，具体包括以下步骤：
s1获取源dem的元数据metadata；
s2建立所述源dem的局部坐标系，对所述源dem进行切分形成若干个可并行计算的切片tile；并根据所述源dem的参数计算每一个所述切片tile的索引值记为key，以<key，tile>键值对为计算单元组成一个保存源dem的元数据的链表list<key，tile>，将<key，tile>按照索引值key分布到不同的计算节点进行计算；
s3根据是否定义了最大层次maxzoom，进行设置或计算最大层次maxzoom，设置需要的输出切片tile的参数；根据设定的坐标系和切分规则计算输出切片tile的参数；
s4计算节点并行构建所分配的键值对<key，tile>，根据索引值key合并同一计算节点上的高程值value；分别根据索引值key按照所述输出切片tile的长、宽、分辨率参数计算所对应的切片tile的外包矩形extend，计算输入切片tile中每个网格单元cell的坐标，再根据设定的输出坐标系，计算每个所述网格单元cell在输出坐标系的坐标；根据每个所述网格单元cell的坐标计算相应的所述网格单元cell所在的输出切片tile的索引值key，并将其加入一个映射map中，组成一个映射map<key，list<cell>>；
s5将计算节点上相同索引值key的链表list<cell>合并；根据索引值key合并所有计算节点上的链表list<cell>，形成一个键值对<key，list<cell>>；
s6根据输出切片tile的参数和坐标系，计算输出切片tile的网格单元cell的坐标，接着根据采样方式从合并的键值对<key，list<cell>>中计算网格单元cell的高程值value；
s7利用8邻域算法，修正输出切片tile边缘的高程值；根据生成的所述链表list<key，tile>，循环获取每个切片tile的四边网格单元cell和与之相邻的8个切片tile，取相邻切片tile的网格单元高程值的平均值作为相应网格单元ell的高程值value，达到切片tile拼接的平滑；
s8判断当前层次zoom是否满足终止条件，若满足终止条件，则结束；若不满足终止条件则根据上一级层次计算下一级层次和下一级层次的切片tile的外包矩形extend的宽度和高度，计算方法：nz=cz–1，wh=2(nz 1)，其中nz表示下一级层次，cz表示当前层次，wh表示切片tile的宽和高度；
s9根据上一级层次切片tile和采样算法，并行计算生成下一级层次切片tile和索引值key；根据索引值key、下一级的层次和切片tile的参数计算当前层次切片tile所对应的下一级层次的索引值key，利用所述采样算法从当前层次切片tile采样下一级层次切片的高程值，构成下一级层次切片tile的子切片和索引值key，在每个节点上合并相同的索引值key，构成键值对<key，list<tile>>；接着根据索引值key，合并所有节点上相同的索引值key的切片tile，组成下一级层次的切片tile，再执行步骤s7；
s10循环执行步骤s8~s9，直到判断层次zoom满足终止条件，所述终止条件为：若设定了层次minzoom，则判断当前层次zoom是否达到设定层次minzoom，若达到，则终止；若未设定层次minzoom，若当前层次zoom小于0，则终止。
2.根据权利要求1所述的基于分布式系统的海量dem金字塔切片并行构建方法，其特征在于，所述步骤s1中所述元数据metadata包括原始分辨率、坐标系，以及在原坐标系下dem的外包矩形extent数据。
3.根据权利要求1所述的基于分布式系统的海量dem金字塔切片并行构建方法，其特征在于，所述步骤s2中建立源dem的局部坐标系具体为：建立以源dem左上角为原点，向下向右为正方向的坐标系；根据设定切片tile的宽和高，将源dem数据切分为多个切片tile，为了将源dem进行初步的切分，以构成链表list<key，tile>，其中索引值key的计算方法：
其中，c、r表示切片tile的行列号，x、y表示切片tile任意一点的瓦片坐标，xmin表示源dem的外包矩形extent的x最小值，ymin表示源dem的外包矩形extent的y最小值，th是初始切片tile的高，tw是初始切片tile的宽。
4.根据权利要求3所述的基于分布式系统的海量dem金字塔切片并行构建方法，其特征在于，所述步骤s3若未设定层次maxzoom，则需要计算层次maxzoom，计算过程是根据源dem的外包矩形extend的左下点坐标和右上点的坐标，并将左下点坐标和右上点的坐标进行投影转换，转换为经纬度坐标，其中左下点坐标是（slon，slat），右上点的坐标是（elon，elat）；接着根据公式计算层次maxzoom，计算公式如下：
其中，pi表示圆周率，maxzoom为最大层次，dlat表示纬度差，dlon表示经度差，lat1表示slat的弧度，lat2表示elat的弧度，a表示计算的中间值，c表示计算的中间值。
5.根据权利要求3所述的基于分布式系统的海量dem金字塔切片并行构建方法，其特征在于，所述步骤s4中网格单元cell=<x,y,value>，其中x，y是在切片tile的坐标，value是高程值。
6.根据权利要求5所述的基于分布式系统的海量dem金字塔切片并行构建方法，其特征在于，所述步骤s4并行计算输入切片tile，根据输出切片tile的参数将输入切片tile的dem值划分到不同的输出切片tile中；其中首先根据输入切片tile的索引值key计算tile的外包矩形extend，计算公式：
其中，minx、miny、maxx、maxy分别表示外包矩形extend的最小最大xy坐标值，col、row是输入切片tile的行列号，exminx、exminy表示源dem的外包矩形extend的最小的xy坐标值，h表示切片tile的高，w表示切片tile的宽。
7.根据权利要求5所述的基于分布式系统的海量dem金字塔切片并行构建方法，其特征在于，所述步骤s7利用8邻域算法，修正切片tile边缘的高程值value时，对输出的链表list<key，tile>并行计算相邻的切片tile，取四周除交点外的相邻网格单元cell值，并对上下或左右的cell取值，将取到的cell值取平均值作为指定网格单元cell的高程值value；四个角点的高程值的计算是取切片tile的相邻的另外的三个切片tile与相应角点的值，并取平均值作为相应网格单元cell的高程值value。
技术总结