CN111553980B - 基于图拉普拉斯下采样技术的dem地形综合方法 - Google Patents

Abstract

本发明公开了一种基于图拉普拉斯下采样技术的DEM地形综合方法，首先提取DEM地形中高程点的三维坐标；其次是具体的DEM地形关键点提取方法，包括使用D8算法提取DEM地形骨架和图拉普拉斯下采样方法提取DEM地形的局部细节特征；图拉普拉斯下采样方法是基于图模型进行的，在这个过程中构造包含DEM高程点间空间关联信息的DT图模型，进而得到对应的初始的图拉普拉斯模型，最终使用拉普拉斯极特征向量方法完成图拉普拉斯下采样操作；为了获取不同粗化层次的DEM地形综合结果，图拉普拉斯下采样操作过程中还需要使用Kron‑Reduction方法获取不同粗化层次的图拉普拉斯模型；最后是使用提取到的高程点重构代表DEM地形特征的三维不规则格网，即可实现多尺度DEM地形综合。

Description

基于图拉普拉斯下采样技术的DEM地形综合方法

技术领域

本发明涉及地理信息科学领域，更具体地说，涉及一种基于图拉普拉斯下采样技术的DEM地形综合方法。

背景技术

数字高程模型(DEM)可以提供多尺度的地理景观信息，因此对于这一基础地理数据的研究是十分有必要的。然而，现阶段对DEM尺度的研究还不够充分，实现不同尺度的DEM地形综合还存在许多问题。例如生成不同尺度的DEM需要考虑地形骨架和地形结构细节，以满足不同粗化层次的综合结果对地形特征需求。现有技术(例如，最大Z-容差法)在生成粗化的DEM时，主要是通过选择具有最小误差的高程点来重建TIN(三角不规则网)。然而，在这一过程中，许多结构细节(如坡度和局部地形特征)却被忽略，最终导致DEM地形综合的精度随着尺度减少而剧烈降低，这在一定程度上削弱了DEM作为基础地理数据的应用价值。

发明内容

本发明要解决的技术问题在于，针对现有技术在生成粗化的DEM时，通过选择具有最小误差的高程点来重建TIN的这一过程中，许多结构细节却被忽略，最终导致DEM地形综合的精度随着尺度减少而剧烈降低，削弱了DEM作为基础地理数据的应用价值的技术缺陷，提供了一种基于图拉普拉斯下采样技术的DEM地形综合方法。

本发明提供的基于图拉普拉斯下采样技术的DEM地形综合方法，首先对目标DEM地形数据的遥感影像进行预处理，提取DEM地形中高程点的三维坐标，以便进行后续的DEM处理工作；其次是具体的DEM地形关键点提取方法，包括使用D8算法提取DEM地形骨架和图拉普拉斯下采样方法提取DEM地形的局部细节特征；图拉普拉斯下采样方法是基于图模型进行的，在这个过程中需要构造包含DEM高程点间空间关联信息的DT(Delaunay Triangle)图模型，进而得到对应的初始的图拉普拉斯模型，最终使用拉普拉斯极特征向量方法完成图拉普拉斯下采样操作；为了获取不同粗化层次的DEM地形综合结果，图拉普拉斯下采样操作过程中还需要使用Kron-Reduction方法获取不同粗化层次的图拉普拉斯模型；最后是使用提取到的高程点重构代表DEM地形特征的三维不规则格网(TIN)，即可实现多尺度DEM地形综合。

本发明提供的基于图拉普拉斯下采样技术的DEM地形综合方法，该方法可以对通过图拉普拉斯下采样技术提取DEM地形结构细节(如坡度和局部地形特征)，再结合DEM的骨架重新构造不同粗化层次的TIN(三角不规则网)进而实现不同尺度的DEM地形综合。本发明适于研究具有更为复杂的地形特征的DEM数据。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明的基于图拉普拉斯下采样技术的DEM地形综合方法一实施例的流程图；

图2是研究区域图；

图3是提取的DEM地形骨架图；

图4是图拉普拉斯下采样图；

图5是Kron-Reduction图；

图6是提取DEM地形结构特征示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

参考图1，本实施例的基于图拉普拉斯下采样技术的DEM地形综合方法，其包括以下步骤：

S1、使用ArcGIS 10.0软件提取研究区域的DEM中地形高程点的三维坐标(X,Y,Z)，具体通过ArcGIS中Toolbox工具来完成，研究区域请参考图2。

S2、提取研究区域的DEM中的地形骨架，所述骨架地形包括山谷和山脊，具体请参考图3中的a部分和b部分。对于山谷，基于步骤S1中提取的三维坐标(X,Y,Z)，使用D8算法进行直接提取；对于山脊，将步骤S1中提取的三维坐标(X,Y,Z)中的Z坐标乘以-1生成初始DEM的反地形，然后使用D8算法进行提取，其中Z坐标表示高程，提取出来的为骨架地形的高程点，具体可参考图3的c部分。提取研究区域的DEM中的地形骨架具体方法包括如下步骤：

S21：对于山谷V_t，计算相邻8个高程点单元之间坡降Slope的公式为：

Slope＝D_Z/D_i(1)；

其中，D_Z是相邻两个高程点之间的高程差值，D_i是相邻两个高程点之间的欧氏距离，l是两个高程点欧式距离在地表平面上的投影长度；

如果任意目标高程点P_t1周围相邻8个高程点(即以目标高程点为正方形的中心，正方形的四个顶点以及四条边的中点共8个点)单元之间最大坡降Slope小于0，则赋以负值以表明此格网方向未定；如果最大坡降Slope大于或等于0，则将对应最大值的方向值作为中心格网处的方向值；最后可以通过D8算法计算出研究区域的山谷走向，即完成DEM地形中山谷的提取。

S22：对于山脊B_S，先将步骤S1中提取的三维坐标(X,Y,Z)中的Z坐标乘以-1生成初始DEM的反地形，然后采用步骤S21中相同的方法进行计算。

S3、提取行研究区域的DEM中的地形细节特征，具体包括如下子步骤：

S31、构造包含DEM中地形高程点之间空间关联信息的DT图模型G_t。步骤S31具体包括：

构造包含地形高程点之间空间关联信息的DT图模型G_t来存储DEM地形的结构特征，构造的DT图模型中的每个顶点代表对应的DEM地形高程点，每一条边都存储相连的两个节点所代表的DEM地形高程点的空间关联关系，得到无向加权图G＝<V,E,W>，G由有限的顶点集V，一组边E和加权邻接矩阵W组成；

其中，高程点P_ta和P_tb之间的加权邻接矩阵W_ab的定义如下：

其中，D_ab是高程点P_ta和P_tb之间的欧式距离，N_Ha和N_Hb分别为高程点P_ta和P_tb的Z坐标值，W_ab为加权邻接矩阵W中第a行第b列的元素。

S32、对于DT图模型，通过拉普拉斯变换获取对应的图拉普拉斯模型

步骤S32具体包括如下步骤：

根据子步骤S31得到的无向加权图G＝<V,E,W>，结合下述公式计算图G的度矩阵D：

其中D_ii表示第i个顶点的度数，n表示无向加权图G中顶点的个数，w_ij为加权矩阵W中第i行第j列的值；然后，根据下述公式构建拉普拉斯矩阵：

其中，

为拉普拉斯矩阵，度矩阵D是列矩阵，D_ii为D中的第i个元素，W是图的加权矩阵。

S33、依据图拉普拉斯下采样方法获取当前粗化尺度下的DT图模型G_t中对应的极大特征值对应着DEM地形细节特征的高程点，从而提取特定尺度下的DEM地形细节特征，图拉普拉斯下采样图如图4所示，DEM地形结构特征示意图如图6所示。其中，DT图模型中对应着DEM地形细节特征的高程点对应图拉普拉斯模型对应的极大特征值；步骤S33中，获取拉普拉斯图模型对应的极大特征值对应在的DT图模型G_t中的高程点的方法如下：

V₁＝V₊:＝{i∈V:u_max(i)≥0}

其中，V₁是保留下来的代表着DEM地形细节特征的高程点，V₊是构造图拉普拉斯模型

中极特征值u_max(i)≥0对应的DT图模型G_t中的高程点。

S34、接下来使用Kron-Reduction图缩减方法获取更加粗化层次的图拉普拉斯模型

具体如图5所示。步骤S34中，Kron-Reduction图缩减方法如下公式所示：

其中，

保留图拉普拉斯模型的全部信息，/>

是某个特定粗化层次保留的节点，/>

表示图缩减后的图拉普拉斯矩阵，/>

表示图拉普拉斯下采样保留的节点/>

对应的图拉普拉斯矩阵，/>

表示图拉普拉斯下采样保留的节点/>

对应的逆图拉普拉斯矩阵，/>

表示图缩减前节点V₁对应的图拉普拉斯矩阵；

进一步，将V1与一个子图

构造如下关联：

G^Kron-reduced＝{V₁,E^Kron-reduced,W^Kron-reduced}

其中，E^Kron-reduced表示图缩减后边对应的矩阵，W^Kron-reduced中任意第i行第j列的元素

的定义如下：

/>

其中，

表示/>

中的第i行第j列的元素，通过Kron-Reduction缩减的图G^Kron-reduced保留了原始图模型的连接关系。

S35、将所述更加粗化层次的图拉普拉斯模型作为新的所述当前粗化尺度下的图拉普拉斯模型，重复S33和S34，每次均获取比前一次更加粗化层次的图拉普拉斯模型过程，直至获取到目标粗化层次的DEM地形综合尺度，最终完成提取行研究区域的DEM中的地形细节特征。

S4、通过使用步骤S2中提取的DEM中的地形骨架的高程点和步骤S3中提取的能够表达特定尺度下DEM地形细节特征的高程点，重新生成新的TIN来表达对应尺度下的三维地形地貌特征。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (7)

1.一种基于图拉普拉斯下采样技术的DEM地形综合方法，其特征在于，包括以下步骤：

S1、获取研究区域的DEM中地形高程点的三维坐标(X,Y,Z)；

S2、提取研究区域的DEM中的地形骨架，所述地形骨架包括山谷和山脊；对于山谷，基于步骤S1中提取的三维坐标(X,Y,Z)，使用D8算法进行直接提取；对于山脊，将步骤S1中提取的三维坐标(X,Y,Z)中的Z坐标乘以-1生成初始DEM的反地形，然后使用D8算法进行提取，其中Z坐标表示高程；

S3、提取行研究区域的DEM中的地形细节特征，具体包括如下子步骤：

S31、构造包含DEM中地形高程点之间空间关联信息的DT图模型；

S32、对于DT图模型，通过拉普拉斯变换获取对应的图拉普拉斯模型；

S33、依据图拉普拉斯下采样方法获取当前粗化尺度下的DT图模型中对应的极大特征值对应着DEM地形细节特征的高程点，从而提取特定尺度下的DEM地形细节特征；其中，DT图模型中对应着DEM地形细节特征的高程点对应图拉普拉斯模型对应的极大特征值；

S34、接下来使用Kron-Reduction图缩减方法获取更加粗化层次的图拉普拉斯模型；

S35、将所述更加粗化层次的图拉普拉斯模型作为新的所述当前粗化尺度下的图拉普拉斯模型，重复S33和S34，每次均获取比前一次更加粗化层次的图拉普拉斯模型过程，直至获取到目标粗化层次的DEM地形综合尺度，最终完成提取行研究区域的DEM中的地形细节特征；

S4、通过使用步骤S2中提取的DEM中的地形骨架的高程点和步骤S3中提取的能够表达特定尺度下DEM地形细节特征的高程点，重新生成新的TIN来表达对应尺度下的三维地形地貌特征。

2.根据权利要求1所述的基于图拉普拉斯下采样技术的DEM地形综合方法，其特征在于，步骤S1中，所述获取研究区域的DEM中各个地形高程点的三维坐标具体是指使用ArcGIS提取研究区域的DEM中各个地形高程点的三维坐标。

3.根据权利要求1所述的基于图拉普拉斯下采样技术的DEM地形综合方法，其特征在于，步骤S2中，所述提取研究区域的DEM中的地形骨架具体方法包括如下步骤：

S21：对于山谷，计算相邻8个高程点单元之间坡降Slope的公式为：

Slope＝D_Z/D_i(1)；

其中，D_Z是相邻两个高程点之间的高程差值，D_i是相邻两个高程点之间的欧氏距离，l是两个高程点欧式距离在地表平面上的投影长度；

如果任意目标高程点P_t1周围相邻8个高程点单元之间最大坡降Slope小于0，则赋以负值以表明此格网方向未定；如果最大坡降Slope大于或等于0，则将对应最大值的方向值作为中心格网处的方向值；最后可以通过D8算法计算出研究区域的山谷走向，即完成DEM地形中山谷的提取；

S22：对于山脊B_S，先将步骤S1中提取的三维坐标(X,Y,Z)中的Z坐标乘以-1生成初始DEM的反地形，然后采用步骤S21中相同的方法进行计算。

4.根据权利要求1所述的基于图拉普拉斯下采样技术的DEM地形综合方法，其特征在于，步骤S31具体包括：

构造包含地形高程点之间空间关联信息的DT图模型来存储DEM地形的结构特征，构造的DT图模型中的每个顶点代表对应的DEM地形高程点，每一条边都存储相连的两个节点所代表的DEM地形高程点的空间关联关系，得到无向加权图G＝<V,E,W>，G由有限的顶点集V，一组边E和加权邻接矩阵W组成；

其中，高程点P_ta和P_tb之间的加权邻接矩阵W_ab的定义如下：

其中，D_ab是高程点P_ta和P_tb之间的欧式距离，N_Ha和N_Hb分别为高程点P_ta和P_tb的Z坐标值，W_ab为加权邻接矩阵W中第a行第b列的元素。

5.根据权利要求2所述的基于图拉普拉斯下采样技术的DEM地形综合方法，其特征在于，步骤S32具体包括如下步骤：

根据子步骤S31得到的无向加权图G＝<V,E,W>，结合下述公式计算图G的度矩阵D：

其中D_ii表示第i个顶点的度数，n表示无向加权图G中顶点的个数，

为加权矩阵W中第i行第j列的值；然后，根据下述公式构建拉普拉斯矩阵：

其中，

为拉普拉斯矩阵，度矩阵D是列矩阵，D_ii为D中的第i个元素，W是图的加权矩阵。

6.根据权利要求4所述的基于图拉普拉斯下采样技术的DEM地形综合方法，其特征在于，步骤S33中，获取拉普拉斯图模型对应的极大特征值对应在的DT图模型中的高程点的方法如下：

V₁＝V₊:＝{i∈V:u_max(i)≥0}

其中，V₁是保留下来的代表着DEM地形细节特征的高程点，V₊是构造图拉普拉斯模型

中极特征值u_max(i)≥0对应的DT图模型中的高程点。

7.根据权利要求6所述的基于图拉普拉斯下采样技术的DEM地形综合方法，其特征在于，步骤S34中，Kron-Reduction图缩减方法如下公式所示：

其中，

保留图拉普拉斯模型的全部信息，/>

是某个特定粗化层次保留的节点，

表示图缩减后的图拉普拉斯矩阵，/>

表示图拉普拉斯下采样保留的节点/>

对应的图拉普拉斯矩阵，/>

表示图拉普拉斯下采样保留的节点/>

对应的逆图拉普拉斯矩阵，/>

表示图缩减前节点V₁对应的图拉普拉斯矩阵；

进一步，将V1与一个子图

构造如下关联：

G^Kron-reduced＝{V₁,E^Kron-reduced,W^Kron-reduced}

其中，E^Kron-reduced表示图缩减后边对应的矩阵，W^Kron-reduced中任意第i行第j列的元素

的定义如下：

/>

其中，

表示/>

中的第i行第j列的元素，通过Kron-Reduction缩减的图G^Kron-reduced保留了原始图模型的连接关系。/>

Images