1 引    言

沟沿线是黄土地貌最重要的一类结构线，研究其空间展布特征不仅有助于揭示黄土地貌的沟谷发育状况和演化趋势，而且也是黄土高原地区土地覆盖状况调查、制定水土保持规划的基础工作与必要条件^[1,2,3]。近年来，广泛使用的高分辨率DEM和不断发展的数字地形分析技术为定量研究黄土地貌结构提供了丰富的数据基础和可靠的技术保障，也使得自动获取连续且具精确定位特征沟沿线的技术成为研究热点^[4]。目前，基于DEM提取沟沿线主要有两种实现思路：一是形态学方法，即由DEM衍生出的坡度、坡向、曲率、水流路径等地形参数进行规则判断进而获得沟沿线的空间分布特征^[5,6,7]；二是边缘检测方法，即将DEM视为灰度数字图像，通过对比、检测其“亮度”突变点以达到提取沟沿线的目的^[8,9,10,11]。然而，受黄土地貌景观复杂性和提取算法适用性等因素影响，虽然基于上述两种研究思路的各类技术手段均能初步实现数字条件下黄土地貌沟沿线的自动提取，但从现有试验结果来看，以DEM为主要数据源的沟沿线自动提取技术在研究方法上仍有待改进和优化。 就形态学方法而言，其主要依据是研究区坡度、坡向、曲率等地形参数的指标差异，因此，这些参数本身的适用性和准确性对沟沿线最终的提取结果具有重要影响。而在明确DEM精度、结构和计算模型的条件下，分析窗口大小成为影响地形参数计算精度的首要因素^[12,13]。如文献^[12]在分析特定DEM分辨率下不同计算窗口对坡度计算的影响规律后发现，坡度计算精度随分析窗口的增大而提高，且应在拟合曲面与地形曲面相匹配的前提下通过适当扩大局部窗口获得(图 1(a))；文献^[13]则在样地实测的验证下，通过比较DEM不同分析窗口与地貌特征识别准确度间的数量关系揭示出DEM精度确定条件下，其衍生品(地形参数)的分析窗口大小与地貌特征识别准确度间的对应关系(图 1(b))。而目前基于形态学沟沿线自动提取方法中有关地形参数计算的局部窗口绝大多数还限制在3×3或5×5水平，这种固定的、无差别的局部窗口计算模式并未顾及地貌特征空间分布的多尺度性和地形参数自身的有效性差异，这也是目前以此为基础的沟沿线自动提取技术还有待改进的一个重要原因。

图 1 DEM分析窗口与地貌特征准确度间的关系和边缘检测自动提取沟沿线结果比较 Fig. 1 The landform accuracy with different neighborhood of DEM and edge detection methods

图选项

而就边缘检测方法而言，尽管该方法能保持一定的识别精度，但由于它们从数学计算角度利用不含地学意义的梯度算子或其他数学方法提取沟沿线，提取过程和算法设计尚未充分考虑沟沿线的空间形态特性，因此，借助这类方法获得的沟沿线一般缺少比较明确的地学含义，提取沟沿线结果包含大量随机噪声，容易出现封闭性、连续性差等问题，而现有源于图像处理的自动连接算法并不能有效解决这一问题(图 1(c)和(d))^[8,9,10,11]。文献^[8]在不考虑地形特征因素的前提下，借助膨胀腐蚀算法对一定阈值范围内的断裂点进行连接，但实际效果并不理想。文献^[9,11]虽然通过引入改进的Snake模型抑制DEM随机噪声，一定程度上改善了以往沟沿线局部断裂点的自动连接效果，但由于该模型及其改进均涉及多次迭代，实际应用中需耗费大量的计算资源和计算时间，从而难以推广至高分辨率DEM或较大的空间范围。

为此，本文在高分辨率DEM数据的辅助下，基于地形特征综合分析思想^[14]，提出了一种借助地形开度(topographic openness)和差值图像阈值分割原理提取黄土地貌沟沿线的技术方法。主要研究目的是，在现有研究基础上探寻新的黄土地貌沟沿线自动提取方法，以进一步改善现有技术方法的适用性与识别精度，优化提取结果与效率。研究结果表明，这一方法不仅强调了地貌基本形态特征在沟沿线提取过程中的作用与影响，同时在一定程度上也实现了沟沿线提取的精确化和自动化。

2 基于地形开度和差值图像阈值分割的沟沿线提取方法 2.1 地形正负开度的概念、特性及其与黄土地貌形态特征的对应关系

地形正负开度(positive and negative openness)由文献^[15]首次提出。其最初设计目的是为继传统地貌晕渲技术之外，寻找新的、能够增强DEM可视性效果的技术手段^[16,17]。根据定义，地形开度通过地面上一点8方向天顶角(zenith angles)、天底角(nadir angles)的平均值衡量一定范围内地形表面相对水平面间凹凸起伏程度的大小，因此，其实质又是一个描述地表形态的定量指标(见图 2)^[18,19]。

开度的计算原理为：若图 2(a)中A为待计算地形正负开度的中心像元，在某一观察半径L范围内，该点正开度(Φ_A,L)、负开度(Ψ_A,L)可由式(1)计算求得

点击浏览公式

点击浏览公式

图 2 地形开度计算示意图 Fig. 2 The topographic openness

图选项

利用上述公式对规则格网的DEM逐像元进行运算便可得到该DEM的正、负开度栅格图像，其中，正开度描述地表凸起程度大小，负开度代表地表凹陷程度大小(见图 2(c)和图 2(d))。

就黄土地貌坡面形态及其分界特征而言，具有明显地理分割界限的正负地形区是目前黄土高原最为常见且普遍被人们所认可的黄土地貌形态，其典型之处在于上下坡面之间存在明显的坡度变化，这些坡面转折点的集合则是构成黄土地貌沟沿线的基本要素。对上下坡面景观形态进行比较，通过判断地表要素集是否具有高程、坡度等差异是目前定义沟沿线的主要依据^[20]，因此，笔者认为以DEM为主要数据源的沟沿线自动提取技术应该以上述形态学特性作为其分析、提取的判别准则和理论支撑。而相比传统方法，地形正负开度能更为精准地量化表达研究区地貌形态及坡面变化的本质特征，主要原因如下：

首先，相比边缘检测中的各类梯度算子，地形正负开度是具有明确地貌形态学含义的特征指标，以此为基础的沟沿线自动提取技术能够最大限度地顾及地表要素集在高程、坡度上的形态差异，避免由于使用无地学含义的边缘检测算子引起的沟沿线随机噪声多、结果连续性差、提取效率低等问题。根据定义，地表上一点正开度值的大小从数量上代表该点上凸(convex-upward)程度，负开度则代表该点下凹(concave-upward)程度大小，其结果不仅使得利用正负开度生成的图像具有良好的可视化效果，而且使其具有对地形坡位精确度量的能力。正开度越大，说明求取开度的位置倾向于正地形；反之，负开度越大，越倾向于负地形(对于山顶点，其正开度将大于90°，负开度小于90°；对于洼地点，其正开度小于90°，负开度则大于90°)。

其次，相比形态学方法中的各类地形参数因子，地表一点的正负开度在一个选定的观察半径L内而非一个固定的分析窗口(3×3或5×5)定量标识地表高程、坡面形态特征，可以有效地顾及地貌特征空间分布的多尺度特性。事实上，地貌形态具有非常强烈的多尺度特性早已是不争的事实^[21,22,23]，使用传统固定分析窗口求取地形参数时，由于比较、分析的范围仅限于中心像元及其临近的有限个高程点(8或24)，这也使计算结果带有明显的局部特性，难以准确反映地貌真实特征(图 3(a))。而地形正负开度的计算则在给定观察半径L内比较并确定中心像元相对于其中最高点、最低点的高程、坡面形态特征差异(图 3(b))，特别在L值较大的情况下，地形正负开度能够定量表达中心像元在“整体”地貌形态中的起伏特征，准确量化表征地貌“多尺度”特性，因此，地形正负开度在一定程度上能够避免传统地形参数计算中使用非最佳特征分析窗口造成的地貌特征识别效果不佳的问题。

图 3 中心像元固定窗口分析 Fig. 3 The landform accuracy with different neighborhood

图选项

以上述分析为主要依据，本文借助高分辨率DEM数据，通过引入地形正负开度因子这一对描述地表形态的定量指标，结合差值图像阈值分割原理，实现黄土地貌沟沿线的自动提取。总体技术路线如图 4所示。

图 4 基于地形开度和差值图像阈值分割的沟沿线提取总体流程 Fig. 4 The overall flow chart of loess shoulder line extraction

图选项

(1) 原始数据的预处理。为减少由高分辨率DEM地形细节反映太过明显而造成的细节误差和数字化造成的误差，需要对原始DEM进行均值滤波等图像处理。

(2) 计算地形正负开度及其差值(DΔ_L)。确定正负开度计算的最大观察半径L后，计算研究区的正负开度值进而获得研究区域上正开度和负开度的栅格图像，其中，正开度反映该区域正地形空间分布状况，负开度则反映该区域负地形的空间分布状况。而二者的差值——开度差值图像由于负值倾向于负地形，正值倾向于正地形，因此，平地和正负地形交界地区的开度差值则集中在0°值附近，从黄土地貌坡面形态及其分界特征来看，若能假定设置一个合适的正阈值t，则当开度差值小于-t时，地貌特征表现为负地形；开度差值大于t时，地貌特征表现为正地形；介于-t与t之间即为平地和正负地形交界区域(见式(3)和图 5)。

图 5 区域地貌高程分布 Fig. 5 The elevation of an area

图选项

点击浏览公式

(3) 最佳阈值选取和差值图像阈值分割。对开度差值图像进行阈值处理可以得到正负地形分布图像，其关键是选择合适的分割阈值t。为此，本文在综合考虑黄土地貌形态明显性和特定尺度两个限定条件下，参照文献^[20]对黄土地貌正、负地形的基本定义结果并根据试验研究区地貌景观实际情况，确定初始分割阈值t为，高原沟壑区正负地形坡度约以15°为界，丘陵沟壑区正负地形坡度则约以25°为界^[20]。考虑到该初始分割阈值往往未必是最佳分割阈值，本文在初始分割阈值基础上采用迭代方法，借助不同分割阈值条件下获得的沟沿线的平均坡度、平均表面曲率和平均上下坡度差等地形因子的统计分析结果^[24]以确定开度差值图的最佳分割阈值。基本原则为：分析每次迭代得到的沟沿线的平均坡度、平均曲率和平均上下坡度差3个量化指标，将平均曲率和平均上下坡度差最大且平均坡度符合研究区域实际情况的阈值确定为最佳阈值。

(4) 提取沟沿线及其后处理。笔者采用数学形态学原理^[25]对二值化开度差值图像进行正负地形的边界提取以得到研究区沟沿线的空间分布特征(见式(4))。然后对得到的沟沿线进行检查，如果提取结果存在自相交或破碎等情况，需进行后处理工作

点击浏览公式

3 试验及结果分析 3.1 试验数据

已有研究表明，对于1∶1万比例尺DEM，5m是保证该地区地形描述精度的理想分辨率尺度^[26]，为此，本次试验使用了购自陕西测绘地理信息局基础地理信息中心的洛川县堡子头5m分辨率DEM数据。数据来源为数字线划图(DLG)，高程中误差为1.24m。该试验区的地貌类型为黄土塬，位于洛川塬地区，具有明显的正负地形特征，利于进行沟沿线提取(见图 6)。

图 6 研究区位置及其地貌特征 Fig. 6 The study area and its geomorphological characteristics

图选项

要计算DEM的地形正负开度，首先需要确定最大观察半径L。图 7为研究区开度不同最大观察半径L得到的正负开度栅格均值曲线。由曲线变化趋势可以看出，开度随着最大观察半径的增大而快速收敛。据此结合DEM的分辨率及计算效率，确定本研究区开度最大观察半径L为50个像元，即在栅格主方向长度为250m，在对角线方向长度为250m。

图 7 最大观察半径L与正负开度的关系曲线图 Fig. 7 The relationship curve graph of lookup distance L and the average of positive openness or negative openness

图选项

参照文献^[20]对黄土正、负地形的基本定义结果并根据试验研究区地貌景观实际情况，本文确定初始分割阈值t为15°(高原沟壑区)，即选择-15°~15°为初始分割阈值的迭代范围。通过迭代得到研究区沟沿线各量化指标的变化趋势如图 8所示，从图上可以看出利用10°~12°的阈值提取得到的沟沿线上下坡度差均值和表面曲率均值取得最大值。而坡度均值在该范围内处在30°左右，符合研究区黄土塬地区沟沿线的坡度特性，据此选定10°为本研究区二值化差值图像的最佳分割阈值。

图 8 各沟沿线量化因子与开度差值阈值的关系曲线图 Fig. 8 The relationship curve graph of the shoulder line quantifying factor and the threshold of openness difference

图选项

利用t=10°作为最佳阈值对开度差值图像进行阈值处理，得到开度差值的二值图像。在二值图像上因为平地的误分，而导致正地形区域产生的孔洞见图 9(a)。利用图 9(b)所示的平地栅格与二值图像进行并操作去除孔洞。另外还需要对正负开度二值化图像进行形态学碎斑消除处理，以去除碎小斑块。图 9(c)为经过处理后的研究区正负地形二值栅格图像。

图 9 研究区正负开度差值图像二值化图像 Fig. 9 The binary image obtained by threshold processing openness difference image

图选项

利用正负地形二值图像所表示的正负地形，基于数学形态学中的边界提取算子(见式(4))对研究区正负地形二值栅格图像进行边界提取得到沟沿线。提取后的沟沿线还需要进行后处理，以去除沟沿线自相交和破碎沟沿线。图 10(a)、(b)分别为上述研究区提取得到的沟沿线及其局部区域放大效果图。

图 10 利用10°阈值提取的研究区沟沿线 Fig. 10 The resultant shoulder line extracted by the threshold 10°

图选项

为比较并检验本研究提取结果的准确度，笔者分别采用LOG算子结合目视解译、地形位置指数(topographic position index，TPI)两种方法对相同试验区提取的沟沿线进行结果对照，综合评判3种方法提取结果的准确度及其相互关系。分析结果如图 10(c)、(d)所示，主要分析结论归纳如下：

(1) 从直观的比较结果不难看出，3种方法获取的沟沿线在空间形态特征上能够保持基本一致，也具有较好的吻合程度，由此可见3种方法均可用于研究区沟沿线的自动获取。

(2) 就LOG算子结合目视提取方法而言，由于提取过程中借助了目视判读解译，在获取沟沿线时人为地考虑了地貌形态特征的细微差别，因此，提取的沟沿线形态最为平滑(图 10(c))，但整个过程较为费时。而就TPI方法而言，其本质基于固定分析窗口的局部地形因子处理技术，因此，提取结果具有较强的局部特性，很难在整个研究区尺度上顾及沟沿线整体形态特征，这也是造成沟沿线提取结果在形态上表现较为复杂曲折的重要原因(图 10(d))。

(3) 基于地形开度方法，由于正负开度本身源于可视线规则，加之其蕴含的尺度灵活性特点，可以在提取过程中适当顾及研究区整体地貌形态分布特征，因此，本文研究方法的沟沿线提取结果在保持其形态特征方面更接近目视解译的提取结果。同时，从提取结果的线性特征来看，本文研究结果较TPI方法获得的线形具有更好的综合性。

(4) 相对于上述两种对比方法的结果，本文的研究结果获得沟沿线较为向正地形方向收缩，也即使用本研究方法得到的地貌蚕食度更大。笔者认为产生这种结果的原因主要受选定开度差值阈值大小的影响。随着开度差值阈值的增大，阈值分割得到的负地形相对于正地形的面积会增大。相反，随着阈值的减小，阈值分割得到的负地形相对于正地形的面积会减小。由此可见，开度差值阈值是影响本论文研究结果准确性的重要因素之一，如何寻找最佳分割阈值、增强最佳分割阈值针对不同地域沟沿线形态的适用性是本文后续研究工作需要着重解决的关键问题。 4 结 论 4.1 地形开度的特征分析及其对地貌特征的保留作用

传统的沟沿线提取方法主要包含两类。第1类是基于数字图像处理的方法，其劣势在于提取的沟沿线缺少明确地学含义并且完整性和连续性较差。第2类是基于局部地形因子表示形态特征的沟沿线提取，其缺点是分析尺度微观，提取的沟沿线事实上并未顾及地貌特征空间分布的多尺度性和地形参数自身的有效性差异。考虑到在明确原始DEM数据精度和地形因子提取模型已经确定的条件下，局部分析窗口大小是决定地形因子有效性的关键因素，而这些因子的准确性又将直接影响沟沿线本身提取结果的精度，因此，针对上述两种研究方法存在的主要问题，本研究首先着重探讨了地形正负开度的概念、特性，并重点阐述了地形正负开度与黄土地貌形态特征的对应关系，以便为本研究后续开展提供必要的理论依据和判别准则。

与局部地形因子分析不同，地形开度的计算在于其设定了一个最大观察半径L而不是一个固定的分析窗口尺寸。因此，与传统的微观地形因子相比，开度的分析窗口会根据中心像元的地理位置和地形特征作出相应调整，其分析尺度也更为宏观，能够自适应地貌形态变化，从而避免了采用任意分析窗口产生的分析误差和对地貌特征的消隐作用。主要因为给定观察半径L内比较并确定中心像元相对于其中最高点、最低点的高程、坡面形态特征差异，特别是在L值较大的情况下，地形正负开度能够定量表达中心像元在整体地貌形态中的起伏特征，准确量化表征地貌多尺度特性，因此，地形正负开度在一定程度上能够避免传统地形参数计算中使用非最佳特征分析窗口造成的地貌特征识别效果不佳的问题

4.2 基于地形开度和差值图像阈值分割原理的沟沿线提取性能表现

结合地形正负开度差值图像阈值分割方法，本文实现了黄土地貌沟沿线自动提取。由于选择合适的差值图像分割阈值是决定沟沿线提取结果准确度优劣的关键，因此，本文借助不同分割阈值条件下获得的沟沿线平均坡度、平均表面曲率和平均上下坡度差等地形因子的统计分析结果，最终确定开度差值图最佳分割阈值。试验结果表明，据此确定的分割阈值能够较好地保存沟沿线形态特征的连续性和完整性，所提取的沟沿线几乎不用手动修改，在一定程度上提高了沟沿线自动提取的工作效率。

通过3种提取方法结果的检验比较，可知本文方法在提取沟沿线过程中能够综合考虑地貌本身形态特征，试验结果在沟沿线线形上优于传统基于局部地形因子提取沟沿线形态特征技术(如TPI方法)，比较接近目视解译提取结果，在提取效率上较目视解译有一定优势。最后，由于开度差值阈值是影响本文研究结果准确性的重要因素之一，因此，如何寻找最佳分割阈值、增强最佳分割阈值针对不同地域沟沿线形态的适用性将是本文后续研究需要着重解决的一个关键问题。