CN113450456B - 高可识别度的dem制作方法 - Google Patents

Abstract

本发明公开了一种高可识别度的DEM制作方法，其包括镶嵌拼接得到研究区影像；对研究区影像进行边缘检测得到待研究区域内的地物轮廓线；获取由地物轮廓线围成的面矢量多边形和地物轮廓线构成的面矢量多边形，并结合分辨率为30m的DEM数据统计每个多边形的特征值；将特征值相似度大于预设阈值时对应的多边形合并为同一分区块；结合地形高程点数据采用地统计学插值方法生成每个分区块对应的高精度DEM数据；对每个分区块的栅格数据进行镶嵌处理，将所有分区块的高精度DEM数据合并为待研究区域的高精度DEM数据。

Description

高可识别度的DEM制作方法

技术领域

本发明涉及DEM制作领域，具体涉及一种高可识别度的DEM制作方法。

背景技术

DEM全称为数字高程模型，是用一组有序数值阵列形式表示地面高程的一种实体地表模型，现阶段在测绘、水文、气象、工程建设、通讯、军事等方面都具有广泛的应用。而针对水利工程，现存的一些水文径流、降雨、防洪减灾等方面的模拟均需要DEM的支撑，地形数据是水文模型模拟中最基本也最核心的部分。

现有市场上存在的DEM产品，其存在一定的误差，若是直接采用现有的DEM产品进行水文模型模拟难以准确模拟实际水文径流过程，使得防洪减灾存在一定的模拟精度误差。对此目前在进行水文模拟的过程中常需要进行实地勘测地形数据，耗费大量的经济成本，以满足精度要求。

而通过地形测量高程点进行DEM的制作过程中，针对高密度地形高程点数据，通常采用TIN三角网进行地形插值，插值效果较好，但在实际转化过程中仍会出现地形坦化的现象，导致地物转折信息缺失。而针对低密度高程点数据，在进行地形插值的过程中，由于数据量过小的原因，地形高程数据在空间上的偏差现象则更为严重。两种方法均难以反映出地形地貌特征，生成的DEM精度难以满足要求，且会对研究过程中的模拟结果造成一定的影响。

发明内容

针对现有技术中的上述不足，本发明提供的高可识别度的DEM制作方法解决了现场勘测地形后直接插值在转折点及轮廓附近形成的DEM误差大的问题。

为了达到上述发明目的，本发明采用的技术方案为：

提供一种高可识别度的DEM制作方法，其包括步骤：

S1、获取待研究区域的多光谱影像和全色影像，并对多光谱影像和全色影像进行预处理，之后将预处理后的影像进行镶嵌拼接得到研究区影像；

S2、采用Sobel算子对研究区影像进行边缘检测，得到待研究区域内的地物轮廓线；

S3、获取由地物轮廓线围成的面矢量多边形和地物轮廓线构成的面矢量多边形，并结合分辨率为30m的DEM数据统计每个多边形的特征值；

S4、计算多边形之间特征值的相似度，并将相似度大于预设阈值时对应的多边形合并为同一分区块；

S5、根据每个分区块中的所有地形勘测点的高程数据，采用地统计学插值方法生成每个分区块对应的高精度DEM数据；

S6、对每个分区块的栅格数据进行镶嵌处理，将所有分区块的高精度DEM数据合并为待研究区域的高精度DEM数据。

本发明的有益效果为：本DEM制作方法能够很好的适用于高密度地形高程点和低密度地形高程点进行DEM的制作过程中，且能够很好的反映出地貌特征。

本方案通过遥感技术进行地物解译，以实现合理的分区，再结合公开的30m空间分辨率的DEM(30m空间分辨率的DEM的获取来源可以参考网址http://www.gscloud.cn/)进行特征值比较，实现对分区数据进行合理的二次验证，以使相似地貌划分至同一区。

根据分区特点将地形插值方法与分区地貌特点进行合理的匹配，可以进一步减少地形插值的波动，最终进行DEM分区的镶嵌合并，得到高可识别度的DEM。

从多个待研究区生成的DEM看，其内部的轮廓线较为清晰，能够很好的反映出地形地貌特征，且能够清晰的反映出不同区域边缘处的转折特征，其在高程参数上展现的细节更能接近实际情况。

本方案得到的DEM能够大大提升了地形的描绘能力，可以较好的应用到工程设计方面，整体效果较佳，精度较高，基于最终生成的DEM能够很好的反映出地形地貌特征，尤其在地物转折区域的高程特征更贴合实际，能够减少现有插值中地物转折区域的误差。

附图说明

图1为高可识别度的DEM制作方法的流程图。

图2为实例一中DEM制作区域。

图3为实例一中DEM效果对比图，其中(a)为直接插值DEM，(b)为高可识别度DEM。

图4为实例二DEM插值区域的遥感影像。

图5为实例二中DEM效果对比图，其中(a)为直接插值DEM，(b)为高可识别度DEM。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

参考图1，图1示出了高可识别度的DEM制作方法；如1所示，该方法包括步骤S1至步骤S6。

在步骤S1中，获取待研究区域的多光谱影像和全色影像，并对多光谱影像和全色影像进行预处理，之后将预处理后的影像进行镶嵌拼接得到研究区影像。

实施时，本方案优选步骤S1进一步包括：

S11、获取待研究区域的多光谱影像和全色影像，并对多光谱影像依次进行辐射定标、大气校正和正射校正及对全色影像进行辐射定标和正射校正；

S12、对正射校正后的多光谱和全色影像进行图像融合，以得到更高精度的待研究区域影像，之后采用分辨率高于多光谱影像和全色影像分辨率的遥感影像对融合后的影像进行几何校正；采用多光谱影像和全色波段影像进行波段融合，进而提升遥感影像的分辨率。

S13、对几何校正后的影像进行镶嵌拼接，得到最终的研究区影像。

在步骤S2中，采用Sobel算子对研究区影像进行边缘检测，得到待研究区域内的地物轮廓线。

在本发明的一个实施例中，所述步骤S2进一步包括：

S21、将研究区影像转换为由R、G、B三波段组成的jp2格式的图片，该种格式的图片可以使转换后图片分辨率不变。

S22、将图片转化为灰度图后，进行归一化处理；

S23、调用sobel算子分别对归一化后的图像进行0°、45°、90°、135°四个方向上的卷积，分别得到图像中每个像素点在四个方向的梯度值：

其中，A为归一化后的像素点的像素值；

S23、根据像素点在四个方向的梯度值，计算每个像素点卷积后的梯度值：

其中，G_0°、G_45°、G_90°、G_135°分别为像素点在0°、45°、90°、135°方向上的梯度值；G_x和G_y分别为像素点在x和y方向上的梯度值；G为像素点卷积后的梯度值；

S24、判断像素点的梯度值是否大于预设梯度阈值；若是，将像素点标记为边缘点，否则将像素点标记为非边缘点；

本方案的预设梯度阈值优选为待研究区域所有像素点卷积后的梯度值的平均值。

为了进一步提高标记的准确性，本方案的预设梯度阈值还可以做如下优化：

S241、选取研究区域的遥感影像中地物类型大于预设数量的区域作为试验区域，并采用人工标记的方式描绘地物轮廓线作为验证样本；

S242、获取待研究区域中与试验区域对应区块提取的地物轮廓线，并采用提取的地物轮廓线构成的试验样本；

S243、分别计算验证样本和试验样本构成的区域的覆盖面积F1和F2，并判断覆盖面积F2与覆盖面积F1的比值是否大于等于面积阈值；

S244、当比值大于等于面积阈值时，输出预设梯度阈值；当比值小于面积阈值时，令初始预设梯度阈值＝2*初始预设梯度阈值，之后返回步骤S243。

S25、采用标记的所有边缘点作为待研究区域内的地物轮廓线，为了便于后续查看，可以把边缘点置为白点，非边缘点置为黑点。

在步骤S3中，获取由地物轮廓线围成的面矢量多边形和地物轮廓线构成的面矢量多边形，并结合分辨率为30m的DEM数据统计每个多边形的特征值。

在本发明的一个实施例中，所述步骤S3进一步包括：

S31、采用S13获得的遥感影像进行地理配准，并定义地理坐标系为GCS_WGS_1984，得到具有空间参考的边缘影像；

S32、将边缘影像中地物轮廓线的像素值赋值为1，将非边缘点的像素赋值为0，并将赋值的边缘影像由栅格转换为面矢量，得到由地物轮廓线围成的面矢量多边形和地物轮廓线构成的面矢量多边形；

S33、根据30m分辨率的DEM数据，统计每个多边形内部高程的特征值均值、方差、标准差、中位数、众数和最值：

X_min＝min(X₁,X₂…X_i)，X_max＝max(X₁,X₂…X_i)

其中，X_min为多边形内高程值的最小值；X_max为多边形内高程值的最大值；M_O为多边形内高程值的众数；为多边形内高程值的平均值；X_i为多边形内的第i个最小高程值，M_d为多边形内高程值的中位数，σ为多边形内高程的标准差，n为多边形内高程的数量；

当n为奇数时，当N为偶数时，/>

在步骤S4中，计算多边形之间特征值的相似度，并将相似度大于预设阈值时对应的多边形合并为同一分区块；

在本发明的一个实施例中，所述步骤S4进一步包括：

S41、将所有多边形构成一个集合，在集合中随机选取一个多边形；

S42、计算选取的多边形与集合中余下多边形间相应特征值的比值：

Q_Ak2＝σ_A/σ_k，Q_Ak3＝M_{O_A}/M_{O_k}

Q_Ak4＝X_{min_A}/X_{min_k}、Q_Ak5＝X_{max_A}/X_{max_k}、Q_Ak6＝M_{d_A}/M_{d_k}

其中，σ_A、M_{O_A}、X_{min_A}、X_{max_A}和M_{d_A}分别为选取的多边形A高程值的平均值、标准差、众数、最大值、最小值和中位数；/>σ_k、M_{O_k}、X_{min_k}、X_{max_k}和M_{d_k}分别为余下多边形中第k个多边形高程值的平均值、标准差、众数、最大值、最小值和中位数；Q_Ak1、Q_Ak2、Q_Ak3、Q_Ak4、Q_Ak5和Q_Ak6分别为多边形A与余下多边形中第k个多边形的平均值比值、标准差比值、众数比值、最大值比值、最小值比值；

之后将两个多边形的6个比值(Q_Ak1、Q_Ak2、Q_Ak3、Q_Ak4、Q_Ak5和Q_Ak6)作为一组数据。

S43、判断每组数据中的6个比值是否全部大于预设阈值，预设阈值可以设定为0.9，若是，则标记余下多边形中的当前多边形，否则，不标记；

S44、当每组数据均与预设阈值比较后，判断是否存在被标记的多边形，若是，进入步骤S45，否则进入步骤S46；

S45、将被标记的多边形与选取的多边形合并成一个多边形，采用未进行合并的多边形更新集合，并随机选取一个多边形进入步骤S42；

S46、采用集合中的余下多边形更新集合，并随机选取一个多边形进入步骤S42。

在步骤S5中，根据每个分区块中的所有地形勘测点的高程数据，采用地统计学插值方法生成每个分区块对应的高精度DEM数据；

在本方案中，当分区块的方差小于预设最小方差时，地统计学插值方法为多项式法；当分区块的方差大于预设最大方差时，所述地统计学插值方法为反距离权重插值法；当分区块的方差位于预设最小方差和预设最大方差之间时，所述DEM数值插值方法为样条函数法。

在步骤S6中，对每个分区块的栅格数据进行镶嵌处理，将所有分区块的高精度DEM数据合并为待研究区域的高精度DEM数据。

下面结合具体的实例对本方案方法制作的DEM的效果进行说明：

实例一：高密度地形数据下高可识别度DEM生成

下面以白洋淀内部村落圈头乡及其内部复杂水域为例，同时结合图2对采用本方案制作的DEM进行说明。

所选区域位于白洋淀湖泊内部，内部水系构造复杂，通过已有勘测的地形高程点进行地形插值难以得到高精度的DEM，且插值后的误差较大，难以在DEM上反映整体地貌轮廓特点。

对此采用本方案提供的高可识别度的DEM制作方法进行白洋淀湖泊DEM的制作：

首先，进行地物轮廓线检测阶段；

接着，统计特征值阶段制作：空间分辨率30mDEM下载、DEM高程点提取、统计特征值比较，从而将该区域分为开阔水体、苇田、田埂和浮岛四大类。

在统计特征比较阶段完成后，需要进行地形勘测点准备工作，本实例应用基于勘测的1:2000地形数据，根据分区结果，分别筛选出开阔水体、苇田、田埂和浮岛四大类内部的高程数据。

接着，进入地形插值阶段，根据分区特点，分别选择不同的插值方法。在此开阔水面区域选用多项式法，苇田区域选用反距离权重法，浮岛和田埂区域选用样条函数法。

最后，进行DEM拼接阶段的制作，将各分区插值后的DEM进行镶嵌为一个整体，空间分辨率为1m，结果见图3中的(b)。

采用现有比较常规的DEM制作方法进行白洋淀湖泊DEM的制作：利用地形高程点进行TIN地形三角网，然后进行三角网格的转换，转换为空间分辨率相同的DEM，结果参见如图2中(a)。

通过对比图2中的(a)和(b)，从直观效果来看，本方案的方法制作的DEM，其内部的轮廓线较为清晰，能够很好的反映出地形地貌特征，且能够清晰的反映出不同区域边缘处的转折特征，其在高程参数上展现的细节更能接近实际情况。

直接插值得到的DEM在地形的表现上则对地形高程坦化程度比较严重，同时丢失了因为地物不同在高程上表现的关键节点。而通过对比直接插值的DEM和本方案DEM两种的数据特征，其中直接插值的最大值为7.46m、最小值为3.84m、平均值为5.97m，其标准差为1.073。而本方案制作的DEM最大值为11.93m、最小值为-0.3181m、平均值为5.78m、标准差为1.294。

由此可见，本方案的高可识别度DEM制作方法可以在一定区域范围内对地形进行很好的精确还原，其最大值和最小值可以在数值上得到更好的还原。

实例二：稀疏典型高程控制点下高可识别度DEM生成

本实例选定在黄河三角洲附近区域，DEM插值区域的遥感影像如图4所示，该区域以农田、堤防和河道为主，在同区域内部地势相对较为平坦。

下面基于测量的高程点对该区域进行地形插值，所测高程点位29个，均为典型控制点。现有技术的插值方法与实例一相同，本方案和现有技术插值后得到的DEM的效果图如图5所示，空间分辨率为10m。

对比图5中的(a)和(b)，高可识别度制作的DEM在地形高程点较少区域仍能很好的反映出地貌形态，其分区特征较为明显，能够大大弥补直接插值所需的高密度高程点成本。

Claims (5)

1.高可识别度的DEM制作方法，其特征在于，包括步骤：

S1、获取待研究区域的多光谱影像和全色影像，并对多光谱影像和全色影像进行预处理，之后将预处理后的影像进行镶嵌拼接得到研究区影像；

S2、采用Sobel算子对研究区影像进行边缘检测，得到待研究区域内的地物轮廓线；

S3、获取由地物轮廓线围成的面矢量多边形和地物轮廓线构成的面矢量多边形，并结合分辨率为30m的DEM数据统计每个多边形的特征值；

S4、计算多边形之间特征值的相似度，并将相似度大于预设阈值时对应的多边形合并为同一分区块；

S5、根据每个分区块中的所有地形勘测点的高程数据，采用地统计学插值方法生成每个分区块对应的高精度DEM数据；

S6、对每个分区块的栅格数据进行镶嵌处理，将所有分区块的高精度DEM数据合并为待研究区域的高精度DEM数据；

步骤S3进一步包括：

S31、采用S13得到的研究区影像进行地理配准，并定义地理坐标系为GCS_WGS_1984，得到具有空间参考的边缘影像；

S32、将边缘影像中地物轮廓线的像素值赋值为1，将非边缘点的像素赋值为0，并将赋值的边缘影像由栅格转换为面矢量，得到由地物轮廓线围成的面矢量多边形和地物轮廓线构成的面矢量多边形；

S33、根据30m分辨率的DEM数据，统计每个多边形内部高程的特征值均值、方差、标准差、中位数、众数和最值：

X_min＝min(X₁,2…X_i)，X_max＝max(X₁,2…X_i)

其中，X_min为多边形内高程值的最小值；X_max为多边形内高程值的最大值；M_O为多边形内高程值的众数；为多边形内高程值的平均值；X_i为多边形内的第i个最小高程值，M_d为多边形内高程值的中位数，σ为多边形内高程的标准差，n为多边形内高程的数量；

当n为奇数时，当N为偶数时，/>步骤S4进一步包括：

S41、将所有多边形构成一个集合，在集合中随机选取一个多边形；

S42、计算选取的多边形与集合中余下多边形间相应特征值的比值，将两个多边形的6个比值作为一组数据；

S43、判断每组数据中的6个比值是否全部大于预设阈值，若是，则标记余下多边形中的当前多边形，否则，不标记；

S44、当每组数据均与预设阈值比较后，判断是否存在被标记的多边形，若是，进入步骤S45，否则进入步骤S46；

S45、将被标记的多边形与选取的多边形合并成一个多边形，采用未进行合并的多边形更新集合，并随机选取一个多边形进入步骤S42；

S46、采用集合中的余下多边形更新集合，并随机选取一个多边形进入步骤S42。

2.根据权利要求1所述的高可识别度的DEM制作方法，其特征在于，所述步骤S1进一步包括：

S11、获取待研究区域的多光谱影像和全色影像，并对多光谱影像依次进行辐射定标、大气校正和正射校正及对全色影像进行辐射定标和正射校正；

S12、对正射校正后的多光谱和全色影像进行图像融合，并采用分辨率高于多光谱影像和全色影像分辨率的遥感影像对融合后的影像进行几何校正；

S13、对几何校正后的影像进行镶嵌拼接，得到最终的研究区影像。

3.根据权利要求1所述的高可识别度的DEM制作方法，其特征在于，所述步骤S2进一步包括：

S21、将研究区影像转换为由R、G、B三波段组成的jp2格式的图片；

S22、将图片转化为灰度图后，进行归一化处理；

S23、调用sobel算子分别对归一化后的图像进行0°、45°、90°、135°四个方向上的卷积，分别得到图像中每个像素点在四个方向的梯度值：

其中，A为归一化后的像素点的像素值；G_0°、G_45°、G_90°、G_135°分别为像素点在0°、45°、90°、135°方向上的梯度值；

S23、根据像素点在四个方向的梯度值，计算每个像素点卷积后的梯度值：

其中，G_x和G_y分别为像素点在x和y方向上的梯度值；G为像素点卷积后的梯度值；

S24、判断像素点的梯度值是否大于预设梯度阈值；若是，将像素点标记为边缘点，否则将像素点标记为非边缘点；

S25、采用标记的所有边缘点作为待研究区域内的地物轮廓线。

4.根据权利要求3所述的高可识别度的DEM制作方法，其特征在于，所述预设梯度阈值为待研究区域所有像素点卷积后的梯度值的平均值。

5.根据权利要求1所述的高可识别度的DEM制作方法，其特征在于，当分区块的方差小于预设最小方差时，所述地统计学插值方法为多项式法；当分区块的方差大于预设最大方差时，所述地统计学插值方法为反距离权重插值法；当分区块的方差位于预设最小方差和预设最大方差之间时，所述地统计学插值方法为样条函数插值法。