丘陵地形特征线提取方法、丘陵地DEM精细化生产方法与流程

丘陵地形特征线提取方法、丘陵地dem精细化生产方法
技术领域
1.本发明涉及特征线提取领域，尤其涉及一种丘陵地形特征线提取方法及丘陵地dem精细化生产方法。


背景技术：

2.全球地理信息资源获取主要采用光学立体影像，通过立体像对密集匹配自动提取地形特征。通过数字表面模型(digitial surface modal，dsm)进行dem生产能够实现大面积区域三维地理信息的快速获取，同时能够弥补境外或无人区等可达性较差区域的数据缺失或分辨率不高等问题。具备立体观测能力的高分辨率国产卫星资源三号(zy-3)，以高空间分辨率、高几何定位精度，谱段丰富的优势，现已经作为全球地理信息资源建设工程的主要数据源，广泛用于全球重点区域的dem生产。
3.高精度dem生产中，丘陵、山地、高山地等地貌的滤波及编辑必须通过地形特征线辅助来实现，因此对特征线的准确快速提取是保证dem成果质量和效率的重要前提。地形特征线包括山脊线和山谷线，作为地形起伏变化的分界线，描述了陆地表面基本地貌形态，在地貌综合、水文模型、流域分析等研究中具有重要作用，在高精度dem生产中更起着极其重要的作用。目前，在生产实践过程中，特征线的提取主要有利用地形几何特征自动提取算法和人工立体采集两种方法。在立体环境下，依靠人工采集的方法虽然提取精度高，但需要耗费较大的人力物力，采集时间较长。基于地形几何特征的自动提取算法虽然自动化程度较高，但山谷植被覆盖遮挡等因素造成提取特征线提取不全或提取错误，容易引起dem编辑错误。因此，在利用dsm数据进行高精度dem生产时，现有研究未能提出针对丘陵地高分辨率dsm数据的特征线快速准确提取的方法。
4.因此，如何提供一种高精度的丘陵地形特征线提取方法及丘陵地dem精细化生产方法是本领域技术人员亟需解决的问题。


技术实现要素：

5.有鉴于此，本发明提供了一种丘陵地形特征线提取方法及丘陵地dem精细化生产方法，对比人工采集山谷线和自动提取结果的空间一致性和重叠度；统计对比不同图幅自动提取方法和传统的人工采集方法所需时间；在平面和立体环境下，检查利用提取特征线优化编辑后山谷区域的滤波效果；采用检查点进行质量评定，检测山谷线上dem成果数据的高程值是否满足技术规程的数学精度。
6.为了实现上述目的，本发明提供如下技术方案：
7.一种丘陵地形特征线提取方法，具体步骤包括：
8.构建特征集和样本库，提取山谷处的植被覆盖范围；
9.运用距离变换和曲率分析方法得到山谷特征线的初步提取结果；
10.利用拓扑分析方法进行特征线的连续性优化，同时利用局部空间关系和经验知识移除异质线段，得到优化后的山谷特征线。
11.可选的，在上述的一种丘陵地形特征线提取方法中，所述特征集中包含纹理特征、高程特征和光谱特征，其中根据dsm和dom成果数据进行预处理分别得到高程特征和光谱特征；所述样本库存储待提取植被区的所有植被样本。
12.可选的，在上述的一种丘陵地形特征线提取方法中，提取山谷处的植被覆盖范围具体步骤包括如下：
13.获取遥感图像分类所需的特征集，包括光谱特征、纹理特征、高程特征，波段叠加得到叠加图像；
14.对所述叠加图像进行基于非监督的分割，生成分割对象；
15.获取遥感图像分类所需的样本数据；
16.将特征集和样本做为监督分类的输入进行面向对象的图像分类，进而得到目标区域的植被面分布图，进而提取得到山谷处的植被覆盖范围。
17.可选的，在上述的一种丘陵地形特征线提取方法中，山谷特征线提取步骤如下：
18.将面向对象的图像分类结果作为输入，利用欧式距离，计算山谷植被区域内所有像素到范围边界的最小距离值，得到距离场数据；
19.利用曲率分析方法对所述距离场数据求导，计算得到梯度数据，然后对梯度数据进行重分类和掩膜运算，提取梯度极大值区域，获得山谷特征线的初步提取结果。
20.可选的，在上述的一种丘陵地形特征线提取方法中，所述距离场数据具体计算步骤如下：
21.利用欧式距离计算二分类图中山谷处的植被覆盖范围内所有像素到范围边界的最小距离值，得到距离场数据；
22.de(x，y)＝[(x
2-x1)2+(y
2-y1)2]
1/2
；
[0023]
式中，xi为经度坐标值，yi为纬度坐标值。
[0024]
可选的，在上述的一种丘陵地形特征线提取方法中，所述曲率分析方法具体步骤如下：
[0025]
计算距离场数据的梯度大小s，取梯度较大值点连接得到骨架线；
[0026]
s(x，y，de)＝[(d
de
/d
x
)2+(d
de
/dy)2]
1/2
；
[0027]
式中，d
de
为距离场数据的导数，dx为经度方向上的导数，砂为纬度方向上的导数。
[0028]
可选的，在上述的一种丘陵地形特征线提取方法中，对所述山谷特征线的连续性优化，引入栅格扩张法，填充缺失值，缓解多个小型闭环状骨架线的问题；进而，引入拓扑细化法，迭代剥离最外层像素，提取出单个像素宽度的山谷线；对所述单个像素宽度的山谷线进行缓冲区分析；然后引入悬挂节点的策略，移除非特征线，进而得到优化后的山谷特征线。
[0029]
一种利用丘陵地形特征线提取方法的丘陵地dem精细化生产方法，具体步骤包括：
[0030]
对dsm分幅成果数据进行滤波编辑，将地表高程降至地面高程，从而得到dem初步滤波结果；
[0031]
对初步滤波结果进行优化处理，对所述山谷特征线，利用林窗量测植被高度并赋予特征线高程属性，并通过三维特征线对dem数据进行二次滤波编辑处理，得到地貌形态特征表达合理的dem成果数据。
[0032]
经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种丘陵地形特
征线提取方法及丘陵地dem精细化生产方法，利用地表覆盖特征与地貌形态特征的关联性，提出一种适用于丘陵地形的特征线自动提取方法，替代传统的人工立体采集作业方式。通过对不同区域不同类型的八幅1:5万标准分幅数据开展生产试验，从定性和定量两个角度证明了应用本发明所提方法获取的山谷线几何位置精度高，连通性好，利用特征线优化编辑后的dem成果数据与立体模型套合良好，地貌形态表达合理，数据质量良好，完全满足数学精度要求。本发明从生产的实际需求出发，借助遥感和gis的技术手段，对dem关键生产环节提出一种新的作业方法，在保证产品精度的前提下，提高生产效率，减少人工作业量。该方法可为全球其他重点区域的高精度dem生产提供参考借鉴。
附图说明
[0033]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
[0034]
图1为本发明的实施例试验数据地理分布和地形地貌特征示意图；
[0035]
图2为本发明的实施例试验数据描述示意图；
[0036]
图3为本发明的丘陵地dem生产流程图；
[0037]
图4为本发明的arcgis python工作流流程图；
[0038]
图5为本发明的精度评价示意图；
[0039]
图6为本发明的人工采集山谷线与提取结果对比示意图；
[0040]
图7为本发明的人工采集与自动提取时效对比示意图；
[0041]
图8为本发明的丘陵地山谷植被滤波效果图；
[0042]
图9为本发明的dem成果质量评定结果图。
具体实施方式
[0043]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0044]
本发明实施例公开了一种丘陵地形特征线提取方法及丘陵地dem精细化生产方法，利用地表覆盖特征与地貌形态特征的关联性，提出一种适用于丘陵地形的特征线自动提取方法，替代传统的人工立体采集作业方式。通过对不同区域不同类型的八幅1:5万标准分幅数据开展生产试验，从定性和定量两个角度证明了应用本发明所提方法获取的山谷线几何位置精度高，连通性好，利用特征线优化编辑后的dem成果数据与立体模型套合良好，地貌形态表达合理，数据质量良好，完全满足数学精度要求。本发明从生产的实际需求出发，借助遥感和gis的技术手段，对dem关键生产环节提出一种新的作业方法，在保证产品精度的前提下，提高生产效率，减少人工作业量。该方法可为全球其他重点区域的高精度dem生产提供参考借鉴。
[0045]
本实施例选择中非共和国(以下简称“中非”)丘陵地区作为试验区，中非是非洲大
陆中部的内陆国家，国土面积约为62.3万平方千米，范围介于14
°
25
′
w～27
°
28
′
w和2
°
13
′
n～11
°
01
′
n之间(图1)。测区内海拔高度约为500-900米，约80％的地形为低山丘陵，地形起伏较小，山脊特征不明显，山谷纵深较小并向两侧山坡面平缓过渡。测区属于热带草原气候，年平均气温约26℃，5月到10月为雨季，11月到次年4月为旱季。丘陵地表覆盖以稀疏丛林及季节性河沼为主，山谷处分布茂密灌木和乔木。
[0046]
需要了解的是：zy-3立体测图卫星影像的平面和高程精度与国际同类型卫星相当甚至更优，完全满足1:5万测图精度需要，现已作为全球地理信息资源建设工程的主要数据源。
[0047]
本实施例所用数据包括zy-3立体像对和全球地理信息资源建设工程中的dsm、dom成果数据，坐标系统为wgs-84，存储单元为1:5万标准分幅。
[0048]
本实施例采用自动匹配和人工交互式编辑方法，获得与立体模型套合、地貌晕渲状态合理、满足数字精度要求的dsm成果数据，并将其作为dem生产的初始数据和面向对象分类输入特征。
[0049]
dom成果数据由原始全色和多光谱影像经过正射纠正、融合、匀光匀色及裁切等预处理生成，用于目视解译参考和面向对象分类输入特征。
[0050]
考虑到测区内丘陵地形的复杂性和多样性，为了充分验证本研究所提方法的可行性，选择八幅1:5万标准分幅数据，位置分布均匀，涵盖不同的地形地貌特征类型(图1)。根据山谷区域和非山谷区域的地物要素特征，大致将中非丘陵地貌划分为四类，山谷区植被覆盖类型包括茂密乔木和较稀疏灌木，非山谷区地物要素包括稀疏草地和较稀疏灌木，具体见图2。
[0051]
需要了解的是：考虑到中非丘陵地形的地表覆盖特征具有相似性和特殊性，在遥感影像上，山谷区域的植被生长茂密且以乔木为主，呈现带状分布特征，山坡上植被稀疏且以草本植物为主，山谷区域和非山谷区域存在明显光谱、纹理和地表高程等特征差异，利用遥感解译方法，能够快速准确地提取出山谷范围。再者，利用gis空间关系分析可以实现山谷线的自动化提取。基于上述思路，本研究设计了中非测区丘陵地形的特征线自动提取流程，如图3所示，主要包括面向对象分类、特征线提取以及基于特征线的dem二次滤波编辑三个步骤，并将其应用于dem生产流程中。
[0052]
具体流程如下：
[0053]
1)影像分割和面向对象分类。根据dsm和dom成果数据，构建特征集(包括光谱特征、纹理特征和高程特征)和样本库，提取山谷处的植被覆盖范围；
[0054]
进一步，构建特证集具体步骤如下：
[0055]
波段提取：从多光谱四波段数据(包含红波段、蓝波段、绿波段和近红外波段)中抽取出三个波段，重组成新的多光谱数据(包含红波段、蓝波段和绿波段)。
[0056]
重采样：采用均值聚类方法，将2米格网间距的数据重采样至10米格网间距。
[0057]
纹理特征计算：基于多光谱数据，利用灰度共生矩阵计算生成纹理特征。
[0058]
镶嵌：将相邻几景的原始卫星影像，采用无缝镶嵌技术生成镶嵌线，利用镶嵌线完成的数字表面模型数据镶嵌。
[0059]
裁切：按照gbt 13989-2012规定，根据全球地理信息资源建设工程中要求，裁切范围按标准1:5万分幅内图廓四个角点坐标的最小外接矩形向外扩展50个格网间距。
[0060]
像素位转换：将图像的像素值范围由32位转换成16位深度。
[0061]
更进一步，利用遥感图像分类方法提取山谷处的植被覆盖范围(植被面)，主要步骤如下：
[0062]
步骤s1:获取遥感图像分类所需的特征集，包括光谱特征(来源于数字正射影像成果数据，简称dom成果数据,含3个波段)、纹理特征(基于数字正射影像，采用灰度共生矩阵计算得到，含1个波段)、地表高程信息(来源于数字表面模型成果数据，简称dsm成果数据,含1个波段)，波段叠加成一个含5个波段的图像。
[0063]
步骤s2:对上述图像进行基于非监督的分割，生成分割对象。
[0064]
步骤s3:获取遥感图像分类所需的样本数据。
[0065]
步骤s4:将特征集和样本做为监督分类的输入进行面向对象的图像分类，进而得到目标区域的植被面分布图，进而提取得到山谷处的植被面范围。
[0066]
2)山谷特征线提取。运用距离变换和曲率分析方法得到山谷线的初步提取结果；
[0067]
3)优化特征线初步提取结果。利用拓扑分析等方法保证特征线的连续性，利用局部空间关系和经验知识移除异质线段；
[0068]
4)对dsm分幅成果数据进行滤波编辑，将地表高程降至地面高程，从而得到dem初步滤波结果；基于特征线对dem进行二次滤波编辑。赋予二维特征线高程信息，并利用三维特征线优化植被初步降高结果，合理表达山谷地貌形态；
[0069]
5)特征线提取结果的几何位置和dem滤波效果的精度评价。
[0070]
为了进一步优化上述技术方案，面向对象分类的基本处理单元是通过分割算法生成的对象，相较于传统的以像素为最小处理单元的分类，面向对象分类更适用于高空间分辨率影像地物识别。通过影像分割，整景影像被分割成独立的区域或对象，具有相似光谱和空间属性的像素被合并至一个对象内，相邻对象之间差异性增大，从而能够有效去除噪声干扰，提高分类精度。面向对象分类主要包括分割和分类两个步骤。采用现有技术中软件，软件中面向对象分类模块提供了针对高空间分辨率遥感影像的分割和分类算法。
[0071]
通过对分割特征的权重值设置，将分割过程应用于zy-3 dom和dsm成果数据。影像的分割尺度大小根据植被带状分布宽度调整，一般来说，宽度越大，其分割尺度参数值越大。特征集包括光谱、纹理以及地表高程特征。本实施例通过测试关键参数的不同组合方案并分析每个图幅的山谷区提取效果以确定各参数值。此外，本发明实施例通过交互式界面构建样本库，可根据预测结果对样本库进行迭代修改以及再预测。
[0072]
为了进一步优化上述技术方案，将距离场数据看做是地形高程数据，高程变化斜率最大值像素连接而成的特征线即为骨架线。利用欧式距离计算二分类图中山谷范围内所有像素到范围边界的最小距离值，得到距离场数据de，如式(1)所示。
[0073]
de(x，y)＝[(x
2-x1)2+(y
2-y1)2]
1/2
；
ꢀꢀ
(1)
[0074]
式中，xi为经度坐标值，yi为纬度坐标值。
[0075]
采用地表曲率分析方法，计算距离场数据的梯度大小s，取梯度较大值点连接得到骨架线，如式(2)所示。
[0076]
s(x，y，de)＝[(d
de
/d
x
)2+(d
de
/dy)2]
1/2
；
ꢀꢀ
(2)
[0077]
式中，d
de
为距离场数据的导数，dx为经度方向上的导数，dy为纬度方向上的导数。
[0078]
通过以上方法实现山谷线的初步提取。然而，初步提取的山谷线并未考虑实际的
地表几何形态，提取结果主要存在空间连通性不高、特征线误判的问题，因此对其进行优化编辑必不可少。空间连通性不高主要表现为相邻特征线之间不连通、多个闭环状骨架线相连两种情况。前者是由于地表曲率分析结果中，骨架线上存在梯度相对较小的区域，导致阈值判断难以将其划分正确。因此在相邻线段端点坐标距离较短的前提下，本实施例引入格栅扩张分析策略，使相邻线段之间连通。针对后者的多闭环线段情况，引入拓扑细化法，迭代剥离最外层像素，从而提取出由单个像素宽度组成的山谷线。由于山区范围轮廓复杂，基于距离变换的方法提出的骨架线存在一些较短的分支，但根据地貌形态的数学定义，这些分支并非山谷线。本实施例通过遍历所有线段端点，判断线段是否存在悬挂节点的策略，移除非特征线段，得到矢量山谷线。
[0079]
由于提取的山谷范围会影响山谷线提取精度，为了改善山谷线提取精度，可根据经验知识对提取结果进行交互式后处理。本研究使用arcgis python创建可重复工作流如图4所示，满足特征线自动化提取的需求。通过对工作流中各模块数据处理速度和效果进行多次实验对比，进一步缩短山谷线提取时间。
[0080]
为了进一步优化上述技术方案，对dsm分幅成果数据进行滤波编辑，将建筑、桥梁和植被等地表高程降至地面高程，从而得到dem数据。在进行丘陵地区的dem生产时，还需在保证数学精度符合技术规定的前提下，正确表达山谷的地貌形态特征。由于中非丘陵地区的山谷及邻近坡地分布茂密灌木和乔木，常规滤波编辑后的dem数据会出现明显的毛刺、跳点等噪声，无法正确表达山谷形态特征。因此，需要对初步滤波后的结果进行优化处理，即采集体现重要山谷形态的特征线，利用林窗量测植被高度并赋予特征线高程属性，并通过三维特征线对dem数据进行二次滤波编辑处理，修正噪声影响，得到地貌形态特征表达合理的dem成果数据。
[0081]
需要了解的是：人工立体采集方法是目前生产中使用较多的一种特征线提取方法，是一种被普遍接受并且精度相对较高的方法。本实施例将人工采集的山谷线定义为参考特征线，采用基于混淆矩阵的精度评价方法，分析自动提取山谷线和参考特征线的几何位置关系。该方法的主要思想是计算提取山谷线和参考特征线分别落于两者缓冲区内和缓冲区外的像素个数，并根据各统计结果构建参考指标(图5)。具体方法如下。
[0082]
首先以参考特征线(人工采集结果)为匹配标准，建立其半径为1个像素或2个像素的缓冲区，判断自动提取山谷线和缓冲区的空间关系，混淆矩阵如表1所示。若自动提取结果落入缓冲区中，则认为其与参考特征线相匹配，并将落入缓冲区内的像素记为tp1(true positive)，反之则认为其与参考特征线不匹配，并将缓冲区外的像素记为fp(false positive)；其次，以自动提取山谷线为匹配标准，建立缓冲区，并将落入缓冲区内的参考特征线像素记为tp2(true positive)，反之则认为其与提取山谷线不匹配，记为fn(false negative)。
[0083]
表1混淆矩阵
[0084][0085]
基于混淆矩阵分析，可定义查准率(correct)、查全率(complete)两个参数，其中查全率反应了自动提取山谷线的完整度，查准率反应了自动提取山谷线的准确性，通过结合查准率和查全率两个参数可以综合评价山谷线的提取质量(quality)，分别如式(3)—(5)所示。
[0086][0087][0088][0089]
采用多种检验方法对山谷线提取精度及其在dem精细化生产上的应用效果进行分析与评价，具体包括：1)对比人工采集山谷线和自动提取结果的空间一致性和重叠度；2)统计对比不同图幅自动提取方法和传统的人工采集方法所需时间；3)在平面和立体环境下，检查利用提取特征线优化编辑后山谷区域的滤波效果；4)采用检查点进行质量评定，检测山谷线上dem成果数据的高程值是否满足技术规程的数学精度。
[0090]
根据山谷线提取精度评价方法，定量评价本实施例提出的山谷线自动提取方法的位置精度，结果如表2所示(将缓冲区内的像素记为匹配，缓冲区外的像素记为非匹配,计量单位为像素px，像素分辨率为10米)。本研究对每幅测试数据分别生成半径为1像素和2像素的缓冲区，并分析其精度结果。
[0091]
表2山谷线提取结果精度分析
[0092][0093]
分析对比表2可知：1)当缓冲区半径为1像素，即自动提取的山谷线的几何定位精度为1像元时，nb33e020014、nb34e017017两景影像的提取精度最高，其质量分别为70.42％和71.60％，两景影像的山谷区域与非山谷区域的地表覆盖差异显著，其中nb34e017017影像中山谷线上植被分布对称性更好。值得注意的是，本发明实施例使用的dem格网间距为10米，该缓冲区约束条件较为严格；2)当缓冲区半径为2个像素时，所有影像的提取结果的查全率和查准率均高于87％，质量均高于77.8％，表明自动提取山谷线的精度(几何精度、准确性)较高。
[0094]
图6是本文方法提取的山谷线与人工采集山谷线的对比。从图6可以看出，自动提取的结果与人工采集结果的几何位置基本一致，体现重要山谷形态的特征线均完整提取，并且具有良好的连通性。自动提取方法存在少量的错提和漏提现象，主要表现为山谷主线附近小分支的遗漏提取，或是由于植被面提取错误造成的山谷线提取错误，因此需要通过少量的人工干预来修正错误。
[0095]
本实施例统计了所有试生产数据的特征线人工采集和自动提取时间，结果如图7所示。从图中可以看出，单景影像的人工采集时间是自动提取方法用时的2～4倍。对比影像累积处理时间，人工采集方法和自动提取方法的所需总时间的差值随着处理影像数量的增加而不断增大。因此，使用自动提取方法批量化提取山谷线较传统生产方式而言，生产效率大大提高，优势明显。
[0096]
进一步，应用自动提取特征线的丘陵地植被滤波效果如图8所示。在dsm成果中，山谷区域的植被形态明显，山谷本身地貌形态未体现。利用特征线编辑后的dem成果，灌木、乔木的冠层高程被降至地面高程，植被基本被完全滤除，改善了初步滤波结果的噪声现象。同时，dem数据能够合理体现中非测区内不同地形地貌的山谷特征。此外，在立体环境中检查发现，山谷处的等高线与立体模型的套合情况较好。
[0097]
利用人工采集的检查点进一步评价dem产品质量。检查点是由经验丰富的作业员在立体环境下采集，并利用林窗或邻近裸地量测获取其高程值，每幅图采集15～20个，均匀分布于图幅内山谷处。八幅1:5万标准分幅数据的质量评定结果如图9所示。从图中可以看
出，1类地貌形态(即山谷处覆盖茂密乔木，山坡覆盖稀疏草地或裸地)的检查点高程中误差最小，因为其山谷和山坡的植被覆盖类型差异明显，山谷线自动提取位置较准确。2类和4类的高程中误差相对较小，为2.4～3.1米，少量高程较差大的检查点常分布于植被较稀疏的山谷区域，由降高过多造成，但误差值显著低于最大误差限度(12米)且不会造成高程中误差(6米)超限。第3类地貌中高程较差大的点基本分布于山谷和山坡植被覆盖类型较相似，斑块破碎化较大的区域，如nb33e012012的质量评定结果所示，其高程中误差为3.5米，最大误差为5.9米，该类地貌的山谷自动提取结果需要后处理修正，也可在立体环境下对滤波结果进一步优化编辑。
[0098]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0099]
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。