基于均值变点分析法的地形线提取研究*

李廷豪，仲济强

(南京海河测绘科技有限公司，江苏 南京 210008)

0 引言

数字高程模型(DEM)的出现使得地形的表达更加的多元化，人们可以更加快速、简便、精确地了解地形信息。数字高程模型是地面高程信息的数字表示形式，对于地学领域的研究具有广泛的应用价值[1]。与传统的地形图相比，DEM可以用于工程规划建设、景观分析、日照分析等。山脊线和山谷线代表了地形最为重要的结构信息，在水文分析、地形重建等方面具有重要的意义[2]。随着地理信息系统的快速发展，提取地形线的方法也越来越多，提取的准确性和便捷性等方面也得到了一定的提升。主要方法有：基于图像处理技术的方法[3]、基于地形表面几何分析的方法[4-5]、基于地表流水物理模拟的方法[6]以及基于地形表面几何分析与地表流水物理模拟相结合的方法[7]。不同的提取方法具有不同的优势以及缺陷，本文基于地表流水物理模拟方法结合数理统计方法来提取山脊线与山谷线，以期得到更加精准的提取效果。

1 基于地表流水物理模拟方法

水从高处向低处流，在山谷处流量大，形成汇水，在山脊处水量小，形成分水等，因此山谷线与山脊线也被称为合水线和分水线。通过分析地表水流的这种自然规律，结合合水线和分水线的汇水量的累计特征，即合水线(山谷线)的汇水量单调递增，而分水线(山脊线)上的汇水量是不累积的，实现对地形线进行提取，称为基于地表流水物理模拟的地形特征提取方法。

其基本原理是首先从每个点8个方向找出坡度下降最大的方向，然后根据高程从大到小依次计算每个点的汇水量，其中区域汇水量最大的点为汇水线的终结点，从终结点出发，找出与该汇水线的终结点相邻近的点，并且该点的汇水量远远小于终结点的汇水量，同时大于其他点的汇水量，确定这个点为新的汇水线上的终结点，直到判定点的所有邻近点的汇水量差异较小，即可得到矢量化的汇水线数据。通过计算流向汇水线的点构成的区域的边界线得到汇水区域的边界线。大量研究表明，这种方法提取的地形线能保持较好的连贯性，但提取结果的质量与邻域窗口的大小有较大的关系。基于此，本文以该方法为基础，提出采用统计学中的均值变点分析法选取最佳邻域窗口，以提高提取的准确性。

2 基于均值变点分析法的邻域窗口确定

正负地形是两种相反的地形，分别高于或者低于邻区。基于地表流水物理模拟的地形特征提取方法需要对原始DEM进行正负地形的提取。在提取正负地形时首先需要进行邻域分析，计算邻域均值DEM，邻域窗口的大小直接影响均值DEM的计算，从而也就影响地形线的提取结果[8]。邻域分析的基本原理是以待计算栅格为中心，向其周围扩展一定范围，基于这些扩展栅格数据进行函数运算，从而得到此栅格的值。这里采用邻域分析的焦点统计法，通过访问栅格中的每个像元，根据识别出的邻域范围来计算高程平均值，得到均值DEM。

由于邻域窗口的大小直接关系到地形线的提取效果，同时顾及山脊线和山谷线的特点，本文提出采用均值变点分析法来选取合适的窗口。均值变点分析法是一种处理非线性数据的数理统计方法，特别对只有一个变点的情况更为有效。对于一非线性的系统输出序列，若在某未知时刻或空间位置突然发生变化，这一点就称为变点[9]。主要计算过程为：

设有样本序列M，令i=2，3，…，N，第i个样本将序列分成两段：X1，X2，…，Xi-1与Xi，Xi+1，...，XN，分别统计两段样本的算术平均值：

(1)

(2)

计算两段样本的离差平方和之和为：

(3)

计算样本的算术平均值与样本的离差平方和：

(4)

(5)

计算期望值：

(6)

由于变点可以使S-Sn的差距变大，因此可以根据计算S-Sn值变化来判断变点存在的位置。

3 实验及结果分析

3.1 正负地形提取

结合地形起伏程度等相关因素，本文选取安徽省黄山市黄山区进行研究，黄山区位于东经117°50′与118°21′，北纬30°00′与30°32′之间。黄山区整体地势起伏较大，其中最高峰海拔达到1 864 m，最低海拔为85 m。本文将地理空间数据云上下载的SRTM 90 m分辨率的黄山区数字高程数据经过裁剪作为实验DEM数据。首先对原始DEM数据进行正负地形的提取，以获得两种相反的地形来对两种相反的地形线分别进行提取分析。在进行正负地形提取时首先要进行邻域分析来获取均值DEM数据，其中邻域窗口的选择至关重要，邻域窗口选取过大或者过小都影响着山脊线和山谷线的提取，而均值变点分析法的使用则较好地解决了窗口大小的选取问题。不同邻域窗口下均值变点分析法统计结果见表1。

表1 不同邻域窗口下均值变点分析法统计结果Tab.1 Statistical results of mean change point analysis under different neighborhood windows

如表1所示，为确定最佳邻域窗口，将邻域窗口设为2×2一直到30×30依次进行运算，得到不同窗口下的均值DEM。再通过前文公式分别计算样本离差平方和与两段样本的离差平方和。通过计算可以得到S=20.345 6，进一步计算在不同邻域窗口大小下两者的差值。

如图1所示，将格网大小与均值DEM进行拟合发现随着格网面积的变化，均值DEM也单调变化，刚开始均值DEM上升速度较慢，而到了一定点，增长速率急剧上升，这个点称为变点。通过分析S-Sn的变化趋势来确定变点的具体位置，发现当邻域窗口为22×22时，S-Sn的值达到最大，因此在进行邻域分析的时候，将22×22窗口作为最佳邻域窗口来得到最合适的均值DEM，以达到更好的提取效果，提高结果的准确性与精密性。

图1 格网大小与均值DEM及S-Sn的关系Fig.1 Relationship between grid size， mean DEM and S-Sn

通过以上方法计算均值DEM，将原始DEM与均值DEM作差可以得到研究区域的概略DEM数据，对概略DEM数据进行两次重分类可以得到研究区域的正负地形结果。为了验证采用均值变点分析法所得到的正负地形具有更好的效果，本文选取最佳邻域窗口大小的同时选取其他两组不同的窗口大小进行实验，获得不同的正负地形结果进行对比分析，图2(a)、(b)、(c)分别是邻域窗口为15×15、22×22、28×28时得到的正地形结果图，图2(d)、(e)、(f)为对应的负地形结果图。从图中可以发现，不同邻域窗口下的正负地形有所差异，当邻域窗口较小时，破碎点比较多，当邻域窗口比较大时，提取的正负地形又比较过度化，而当邻域窗口适中时，提取的正负地形比较完整。因此本文选取22×22的邻域窗口具有一定的可行性。

图2 不同邻域窗口下的正负地形Fig.2 Positive and negative terrain under different neighborhood windows

3.2 山脊线、山谷线提取

与山脊线不同，在进行山谷线的提取时，首先要提取原始DEM的反地形DEM数据，反地形与正负地形概念不同，它不是客观存在的，只是基于一个基准面将原始地形来进行翻转而得到的地形。由如图3所示，基于ArcMap对原始DEM以及反地形DEM进行洼地填充、水流方向提取、汇流累积量，得到原始DEM汇流累计结果与反地形DEM汇流累计结果，将两个结果分别与正负地形进行相交可以得到研究区域的山脊线与山谷线栅格结果，对栅格结果进行自动矢量化得到的矢量结果图。

图3 山脊线、山谷线矢量化结果Fig.3 Vectorization results of ridge line and valley line

4 结束语

邻域窗口的选择对山脊线与山谷线的提取具有重要的影响，因此本文在基于地表流水物理模拟的基础之上，通过数理统计方法选取最佳邻域窗口来进行山脊线与山谷线的提取。由实验分析可以得到，不同的邻域窗口所得到的正负地形有所差异，通过对均值变点分析法的原理介绍以及应用，得到最佳邻域窗口，该方法提高了地形线提取的准确性，同时这种方法提取的自动化有待进一步研究。