资源环境承载能力评价作为区域发展规划的限制性概念，其研究主题首先是自然资源^[1]。目前，大部分学者的研究集中在地下水资源^[2-3]、矿产资源^[4-5]等地质要素。但是，自然资源也包括覆盖地表的森林、河流、山川等，且这些地表自然资源与地下水、矿产等地质要素联系密切，但这部分自然资源的研究在地质研究领域并不多见。其原因主要有3个方面：①未对地表自然资源予以研究重视；②数据源获取困难，人工数据采集的局限性大，数据不完整，且人为因素干扰较大；③地表自然资源的监测与调查涉及的学科面较广，存在一定的技术难题。

数据获取手段的不断进步，遥感对地观测技术的不断成熟，以及遥感交叉学科的蓬勃发展，为以上问题提供了解决方案。遥感技术可以快速、大范围、非人工干预地直接获取地表覆盖卫星影像数据，这种卫星影像包含丰富的光谱和空间信息。遥感交叉学科领域的学者们，提出了诸多卫星影像数据自然资源识别方法，如植被指数NDVI^[6]、水体指数NDWI^[7]等。这些技术进步使全自动监测地表自然资源成为可能。

山体边界及山体信息提取是山地区域进行诸多研究的前提条件。例如，为研究山区生态效应，韩芳等^[8]针对青藏高原山体范围展开生态环境分析，该研究假设已获得山体几何真值；山体基面高度是影响山体效应最重要和关键的地形因子，在秦巴山地山体效应定量化研究中^[9]，基于DEM提取了山体特征线，未提取山脚线；在山区高速公路两侧山体的形变监测及山体滑坡提取，需首先提取山体外轮廓信息^[10-12]。在山体边界提取技术方面，激光雷达技术、InSAR和倾斜摄影测量技术已应用于三维建模、山体创面建模等，其中涉及山体边界、轮廓等几何信息的提取，但这些方法自动化程度低，需要借助野外控制测量，不适用于大范围山体确界^[10-14]。此外，有研究针对高山地区通过提取检测山体阴影^[15-16]，进而智能搜索山体相关信息，如冰川的识别等^[17]。这种间接识别山体的方法适合特定区域的研究。

其中，韩芳等^[8]围绕青藏高原山体范围展开生态效应分析，但是该研究并没有对山体信息的识别进行探讨，仅在假设已经获得山体几何真值的情况下，进行了生态环境方面的论述。对于变形分析和山体滑坡监测，李勇等^[9]采用InSAR(Interferometric Synthetic Aperture Radar)技术对高速公路两旁的山体进行了形变监测。虽然形变监测属于山体几何信息提取范畴，但是高速公路位于山体内部，其方法并未对大范围山体的外轮廓信息进行提取和监测。李超^[13]借助激光雷达技术在三维建模方面的优势，将激光雷达获取的点云数据应用于山体形变监测任务。陈兴芳等^[10]和杨燕等^[11]分别采用InSAR和倾斜摄影测量技术，对滑坡区域进行建模。上述方法虽然在其研究主题上取得了预期效果，但均存在相同的问题，即不适用于大范围山体确界。倾斜摄影测量技术亦被应用于山体创面建模^[12]，其通过传统野外测量的方法获得崩塌体控制点，实现了精确的几何信息提取，但是该方法自动化程度低，需要借助野外控制测量。杨奇勇等^[14]针对山体阴影开展研究，并对阴影区域的像素灰度值采用克里格方法进行了修复，但是其研究侧重于图像辐射信息的改正，未对山体几何信息提取进行探讨。岳照溪等^[15]围绕山体进行了智能化方法探索，并引入DEM(Digital Elevation Model，数字高程模型)作为辅助，但是该方法主要用于检测山体阴影。都伟冰等^[16]针对新疆陡峭的高山地区展开研究，通过Landsat TM影像提取了山体阴影区域，进而识别冰川。

山体边界的划定已用于生态控制线的研究和规划中。珠海市分类划定生态用地用以构建生态安全格局^[17]，其中生态空间划分为山体、湿地滩涂、滨河地带和农田耕地，而山体范围的确定直接影响了生态控制线的划定。冯泉霖等^[18]等研究了济南市山体保护边界问题，济南城市化进程的加快使得山体被侵占，山体的保护控制成为“生态文明”建设的重要内容，其保护边界的划定分为自然山体边界与人工开采山脚线，然而主要依赖于野外实际调查，未形成理论方法。

综上所述，山体边界的划定在自然资源监测与生态文明建设中有重要的意义。遥感和统计学习技术的进步为研究山体提供了基础，但是目前基于该思路的相关研究仍然较少，尤其是全自动识别山体及山体边界的研究十分匮乏，亟待探索。本文利用30 m分辨率的安庆地区Landsat 8遥感卫星数据，围绕山体的全自动识别和边界提取，设计了一种DEM辅助的、基于统计学原理的智能算法。

1.   研究区与数据

安庆市坐落于安徽省西南部，位于北纬29°47′~ 31°17′、东经115°46′~117°44′之间(图 1)，地处安徽、湖北和江西三省交界, 同时也是长江经济带和“一带一路”双节点城市。安庆北倚大别山，南临长江，地势由西北向东南逐步降低呈三级台地，属于典型山水城市。区内自然资源丰富，森林资源主要位于山地丘陵区，有巨大的生态系统服务价值。自动提取山体边界，对于分析山地自然资源、城市规划建设有重要的意义。

图  1  研究区位置

Figure  1.  Location of the study area

下载: 全尺寸图片 幻灯片

本文使用的实验数据是安庆地区30 m分辨率的Landsat 8遥感卫星图像数据和该区域对应的DEM数据。需要注意的是，本文只采用遥感图像数据所包含的蓝、绿、红3个波段光谱信息。这3个波段图像数据和DEM数据如图 2所示。安庆地区西北部存在一片大范围连续分布的山体区域，东南部存在多个小范围山体，错落分布于东南部各地，无论是大范围山体和小范围山体区域，其边界均极复杂。

图  2  安庆市遥感图像和DEM数据

a—蓝波波段图像；b—绿波波段图像；c—红波波段图像；d—DEM数据

Figure  2.  Remote sensing image of Anqing City and DEM

下载: 全尺寸图片 幻灯片

2.   DEM辅助的山体识别与边界提取

本文采用统计学习领域的聚类分析方法设计山体识别算法。该算法首先利用Landsat 8卫星遥感图像的光谱信息和DEM高程信息作为特征，采用高斯混合模型^[19]在特征空间进行聚类，发掘潜在的位于山区的像素。然后，利用DBSCAN(density-based spatial clustering of applications with noise)算法^[20]，实现图像空间聚类，使离散的山区像素聚合成山体，再利用后处理手段，剔除噪声和小图斑。最后，利用边缘检测算法，提取各聚类的边缘，进而检测山体边界。具体流程如图 3所示。

图  3  山体及山体边界提取流程图

Figure  3.  Flowchart of the mountain bodies and their borders extraction

下载: 全尺寸图片 幻灯片

2.1   基于高斯混合模型的山体识别

该步骤主要包括：图像特征空间构建、高斯混合模型表示、高斯混合模型求解。

(1) 图像特征空间构建

利用Landsat遥感图像蓝波、绿波和红波波段的光谱信息及相应DEM高程信息，构建具有4个维度的特征空间。即在后续的建模过程中，每一个像素均对应一个四维向量，用于刻画该像素的光谱和几何信息。

(2) 高斯混合模型表示

假设特征空间存在2个高斯分布，即山体从属于第一个高斯分布(参数为μ₁，Σ₁)，非山体从属于另一个高斯分布(参数为μ₂，Σ₂)。2个高斯分布的先验服从伯努利分布，因此仅存在一个参数ϕ。则该模型的对数似然函数L(ϕ, μ, Σ)为：

其中，z(i)是隐含变量，表示每一个特征可能的类别；p(?)代表离散的概率值。

(3) 高斯混合模型求解

从高斯混合模型的对数似然函数可以发现，对数函数内部有求和表达式，因此，该似然函数无法利用常规方法实现最大化。因此，算法采用一种迭代求解的方法，逐渐逼近极大值解。

首先，给定隐含变量z(i)的初始值，则式(1)中的似然函数可以简化成如下形式：

然后，得到关于参数ϕ, μ, Σ的极大似然解：

以上过程需要迭代执行，不断估计隐含变量z(i)的值，进而不断更新参数ϕ, μ, Σ的极大似然解。其中隐含变量z(i)的更新公式为：

以上过程收敛后，最大的p (z(i)=j|x(i); ϕ, μ, Σ)值对应的类别标号，即该模型所预测的当前像素类别标号。

2.2   DBSCAN算法的山体区域提取

通过高斯混合模型结算，并识别出来的山体是独立的像素，若提取山体区域，独立的像素需要进行追踪，进而形成区域整体以表示山体结构。DBCAN方法是一种基于密度的聚类方法，可以在含有噪声的集合中追踪出任意形状的簇(区域)，其原理详见参考文献[10]。鉴于DBSCAN方法极小的计算复杂度和内存消耗，本文采用该方法实现独立像素到区域整体的追踪。在实现过程中，由于山体与非山体像素已经在前述步骤中进行了标记，所有山体的像素值是相同的。因此，不需要考虑密度大小问题，也不需要筛选核心点，即任何一个像素均可作为聚类中心。

基于DBSCAN算法的山体区域提取具体实现步骤如下。

(1) 在图像中选择一个从属于山体类别的像素，并创建一个队列，将该像素加入其中。

(2) 根据图像空间相邻关系，检索该像素的相邻像素，若相邻像素中含有同样从属于山体类别的像素，则将这些像素加入到当前队列；若不存在相邻像素从属于山体类别，则返回第(1)步，并将当前队列视为一个山体区域。

(3) 循环执行第(1)、(2)步，直到所有像素均被遍历为止，位于同一个队列的像素具有相同的标号。不同队列的标号按照追踪完成的顺序确定。

2.3   山体边界提取

整个图像覆盖范围的所有山体区域提取完成以后，整幅图像是由几个独立区域组成的，各独立区域由具有一致性标号的像素构成。基于以上分析，本文并不需要借助于目前较全面的边缘检测算子^[21]。这里利用一个具有差分作用的模板来起到区域边界检测的效力，本文使用如图 4所示的差分边缘检测算子。具体边缘检测的实现过程即首先利用图 4所示的算子，以窗口滑动的方式遍历整幅图像；将图像像素与该算子的卷积结果，作为新的像素值，进而生成一幅新的特征图像；最后设定一个经验阈值，特征图像中像素值大于该阈值的像素被认为位于山体边界。

图  4  边缘检测差分算子

Figure  4.  Difference operators of edge detection

下载: 全尺寸图片 幻灯片

3.   实验及结果分析

实验采用定性和定量化分析的方式，验证本文设计方法的有效性。其中，定量化分析的精度评价指标包括正确率(P_c)、误判率(P_e)和漏检率(P_m)，具体计算公式如下：

其中，R_a是本文方法自动提取的山体区域；R_m是人工勾选的山体区域；R_a∩R_m是本文方法自动提取的山体区域与人工勾选的山体区域的交集；R_a∪R_m是相应的并集。

3.1   可视化分析结果

本文方法分为3个步骤，即基于高斯混合模型的山体识别、基于DBSCAN算法的山体区域提取和山体边界提取。图 5是本文方法各步骤识别或提取结果。

图  5  山体识别与边界提取各步骤结果

a—安庆地区DEM信息可视化效果图；b—本文方法对安庆地区山体像素与非山体像素的识别结果；c—本文方法对安庆地区山体区域提取结果；d—安庆地区山体区域边界与遥感图像叠加结果

Figure  5.  Results of the mountain bodies and their borders extraction

下载: 全尺寸图片 幻灯片

图 5-b是基于高斯混合模型的山体识别生成的山体与非山体聚类图。可以直观地发现，该聚类图已经识别出山体像素，其中，绿色赋值的像素位于山体。但是该结果存在较多噪声(聚类图中包含很多离散的绿色赋值的孤立像素或小范围区域)，因此，后续步骤需要剔除这些孤立像素或小范围区域。

图 5-c是基于DBSCAN算法的山体区域提取结果。从该结果可以看到，不同的山体区域被识别并提取出来。并且，该步骤有效剔除了图 5-b中存在的孤立像素或小范围区域。

图 5-d是山体边界提取结果，其中红色曲线标识出了边界。从该结果可以发现，本文方法检测出的边界细致，有效刻画了山体的自然形态。图 6是基于本文方法制作的安庆地区山体地图，局部放大图表明，该方法提取的边界符合人类视觉规律，并与人类视觉规律识别出的山体边界重合。

图  6  安庆地区山体制图

Figure  6.  Map of mountains in Anqing City

下载: 全尺寸图片 幻灯片

3.2   定量化分析结果

表 1给出了本文方法的定量化分析结果。其中本文方法的准确率达到88.61%，对于非监督识别方法而言，该精度较高。此外，本文方法的误判率较低，仅为5.53%，进一步说明了该方法的有效性。虽然漏检率较高(10.76%)，但是对于一种全自动算法而言，取得该表相应高的精度及相应低的误差，已经达到了预期目的。

表  1  定量化结果分析

Table  1.  Analysis result of the quantitative results

评价指标	准确率(pc)	误判率(pe)	漏检率(pm)
精度/误差	88.61%	5.53%	10.76%

下载: 导出CSV 

| 显示表格

4.   讨论

本文设计的山体提取方法，以遥感对地观测技术为支撑，采用非监督识别策略，实现了全自动的山体识别及边界提取。该方法包括3个子部分，其中第一部分是方法主体——基于高斯混合模型的山体识别，其采用统计学习思想，是一种非监督的学习方法，不需要借助人工标记样本对模型进行训练。这一点有别于目前的深度学习^[22]技术，基于深度学习的方法虽然精度高，但其结果严重依赖于大量人工标记样本，即Data hungry问题^[23]。本文使用的数据仅覆盖安庆市区域，数据量有限，难以通过人工标记的方式获取大量训练样本，因此并不适用基于深度学习的监督学习方法。本文方法在未给定标记样本的情况下，实现了山体识别。基于前文的可视化和定量化分析，说明该方法在数据量支撑不充分的情况下，达到了预期目标。

此外，本文方法的第二和第三部分，从离散的山体像素中提取独立的山体结构，并实现山体边缘的精确检测。其中山体区域提取采用DBSCAN思想，算法复杂度为O(N)，其中，N是图像像素个数，即计算复杂度与像素总数N呈线性关系。边界提取采用差分算子，其算法复杂度与DBSCAN相同。山体区域和山体边界提取算法均作用于第一部分的山体像素识别结果，该结果以二值图像的形式存储，像素同质性强，因此不需要额外设计相似性测度，不增加算法复杂度。因此，本文方法不仅达到了预期的识别精度，且计算复杂度低，适用于缺乏人工标记样本的识别任务。

5.   结论

实验采用安庆地区30 m分辨率的Landsat 8遥感卫星图像和对应的DEM作为实验验证数据。根据定性和定量化分析，可以得出以下结论。

(1) 经典的统计学习和图像处理理论，可用于辅助生态地质调查任务。本文提出的山体边界自动化提取技术方法对这些理论的应用研究，有效识别了山体区域，并提取了符合人类视觉特征的山体边界。

(2) 本文方法能够为安庆市的山体保护提供依据，有助于指导城市山体保护和规划建设工作，为安庆市山地区域开发和保护工作奠定了基础。

(3) 本文方法实现了山体确界任务的自动化，有效减轻了山体识别过程中的人为因素干扰，为后续的自然资源承载能力分析提供了可靠的前期支撑。

在未来的研究工作中，更多有效的统计学习、机器学习和模式识别方法将被引入到基于遥感对地观测的地质调查任务中，实现宏观地质调查的全自动化。