地形起伏度最佳统计单元算法的比较研究

陈学兄, 常庆瑞, 毕如田, 刘正春, 张小军

(1.山西农业大学 资源环境学院, 山西 太谷 030801; 2.西北农林科技大学资源环境学院, 陕西 杨凌 712100; 3.山西农业大学 园艺学院, 山西 太谷 030801)

地形起伏度是定量描述地貌形态，宏观描述地形变化整体趋势的主要指标[1-2]，也是导致水土流失的最直接因素，它不仅具有数学意义且具备土壤侵蚀和地貌学的意义，已被广泛应用于各个研究领域，如水土流失定量评价[3]、潜在水土流失评价[4]、地质环境质量定量评价[5]、区域滑坡灾害评价[6]、生态敏感性评价[7]、水土流失敏感性评价[8-9]等。在此基础上，国内学者对地形起伏度最佳统计单元的计算方法有了一定的研究[10-12]，最佳统计单元的确定属于典型的变点分析问题，其计算方法有：人工作图法、最大高差法、模糊数学法、均值变点分析法和累积和CUSUM分析算法，除人工作图法外，其他方法都以定量的形式寻找突变点，哪种方法更适合山西省地形起伏度最佳统计单元的计算还需研究。鉴于此，本研究使用SRTM DEM数据，分析并探讨最佳统计单元各计算方法的准确性及适用性，确定山西省地形起伏度提取的最佳统计单元大小，以期为生态敏感性评价、水土流失定量评价等相关研究提供数据资料，也为山西省地形起伏度的进一步应用研究提供参考。

1 研究区概括

山西省，位于我国华北西部的黄土高原地带，东隔太行山与河北为邻，西隔黄河与陕西相望，南与河南省接壤，北与内蒙古自治区相连，疆域轮廓呈由东北斜向西南的平行四边形，全省总面积约16万km2，地貌类型复杂多样，有平原、台地、山地、丘陵、高原等，其中山地丘陵约占全省总面积的80%之余。山西属温带暖温带、半湿润半干旱大陆性季风气候，降水偏少，年降水量为400～600 mm，自西北向东南递增。境内矿产资源丰富多样，且分布广泛，尤以煤炭资源储量大，质量优，品种齐全，易开采而著称于世。

2 数据与研究方法

2.1 基础数据

本研究使用分辨率为90 m的SRTM DEM数据，数据采样格网大小为90 m，格式为ArcGIS的Grid格式，地理坐标系为：WGS_1984_Albers。

2.2 研究方法

2.2.1 地形起伏度的提取 使用邻域分析法提取山西省地形起伏度，所使用矩形窗口大小为n×n像元(n=2，3，…，30)。地形起伏度提取的步骤如下：(1) 统计n×n(n=2，3，…，30)窗口内像元的最大和最小值；(2) 计算各窗口内最大值和最小值的差值；(3) 统计不同窗口的平均起伏度值(表1)。

2.2.2 最佳统计单元的计算 运用人工作图法[10-11]、最大高差法[10]、模糊数学法[10]、均值变点分析法[1,11,13]和累积和CUSUM分析算法[12]分别计算最佳统计单元，研究各计算方法在确定最佳统计单元上的区别与适用性。

(1) 人工作图法。利用统计软件对平均起伏度进行回归分析，提取拟合度最好的拟合曲线，从拟合曲线上人工判断拐点。

(2) 最大高差法。在(i，j)样本集中，样点i在邻域面积j处的高差变化值(dH)为：dHj=Hj-H(j-1)，比较i点所有的高差变化值dHj，取最大的dHj所对应的j点为拐点，则拐点所对应的面积为最佳统计单元的面积。

(3) 模糊数学法。样本集(i，j)中，样点i在j邻域面积处，有：

(1)

(2)

(3)

式中：模糊系数(ak)均匀分布于(0,1]区间，其中2≥k≤N，N为邻域面积个数(N=29)，a2=1；定义拐点为最大Rj值对应的那一点。

(4) 均值变点分析法。均值变点分析法对恰有一个变点的检验最有效，该方法已被广泛应用于计算地形起伏度最佳统计单元，根据均值变点分析法的原理构建样本序列，具体过程为[11]：首先，计算各统计单元下的单位地势度(T)，单位地势度(Ti)=平均起伏度(ti)/邻域面积(si)，(i=2,3,…,30)；其次，对单位地势度(T)取对数，得到序列X，X为{xk,k=1，2，3，…，29}；最后，根据下面公式对表1数据进行计算处理，计算出统计量S和Si的值。

(4)

(5)

(6)

3 结果与分析

3.1 地形起伏度与邻域面积关系分析

计算不同统计单元的邻域面积(表1)。由表1可知，平均起伏度值随邻域面积的增大逐渐增大，但当邻域面积达到一定的阈值时，其变化趋势逐渐平缓。

3.2 山西地形起伏度最佳统计单元的探讨

3.2.1 最佳统计单元的确定

(1) 人工作图法。利用统计软件对表1中的平均起伏度值进行回归分析，提取拟合度最好的2条拟合曲线，见图1。

表1 统计单元与地形起伏度对应关系

图1 邻域面积与平均地形起伏度对应关系拟合曲线

由图1可知，使用SRTM DEM数据分析研究山西省地形起伏度时，采用幂模型和对数模型对平均起伏度值和邻域面积进行拟合，幂模型的拟合效果最好，而对数模型次之，拟合方程分别为y=22.777x0.3853，y=46.95ln(x)-60.119；平均起伏度值随邻域面积的增大逐渐增大，但其增幅一直在降低并渐渐趋于平缓，由于从拟合曲线上人工判断拐点容易受主观因素的影响，所以此方法不参与比对分析。

(2) 最大高差法。根据最大高差法的含义，由表1数据得知平均起伏度的最大变化值在3×3处，即3×3(0.072 9 km2)为使用此方法计算得出的最佳统计单元。

(3) 模糊数学法。根据模糊数学法的概念，各变量的取值为ak=1，27/28，…，1/28(k=2，3，4，…，N)，N=29。计算结果见表2。

表2 模糊数学法统计结果

由表2可知，在j=3时，Rj的值为最大，即使用该方法计算得到的最佳统计单元为3×3(0.072 9 km2)。

(4) 均值变点分析法。根据均值变点分析法的计算过程，计算得出统计量S=22.23，Si的值见表3。

由表3可以得知，在第10个点时S与Si的差值达到最大，这一点就是所要求的点，即利用该方法计算得到的最佳统计单元为11×11(0.980 1 km2)。

(5) 累积和(CUSUM)分析算法。 为确定其突变情况，依据累积和算法分析突变点的原理对表1中的数据序列进行了分析，分析原始数据序列中累加和Si与统计单元窗口的关系，见图2。由图2可知，随着邻域面积的增大，CUSUM值逐渐减小，最小值(-827.34)出现在14×14处，然后随着邻域面积的增大，累积和(CUSUM)的值也逐渐增大，因此可确定14×14(1.587 6 km2)为原始数据序列发生突变的位置。

表3 均值变点分析法的统计结果

图2 平均起伏度的累积和CUSUM图

3.2.2 计算结果对比分析 不同计算方法计算所得的最佳统计单元面积见表4。与已有研究成果[13]进行比对分析，以验证不同计算方法其结果的准确性。在ArcGIS支持下，首先将不同计算方法所得的最佳统计单元下的山西省地形起伏度分为5个等级：平坦(0～30 m)、微起伏(30～70 m)、小起伏(70～200 m)、中起伏(200～500 m)和大起伏(>500 m)，统计各级所对应的像元数，并计算各级所占的面积百分比(表5)，然后依据《中国1∶100万地貌制图规范(试行)中的地貌基本形态划分体系》及其他相关研究成果[14-15]，对基本地貌类型进行了简单的划分，并计算各类地貌所占的面积百分比，结果见表5。

表4 不同方法得到的统计单元面积结果对比

表5 不同方法计算结果分类

为比对不同方法其计算结果的准确性，使用SPSS统计软件分别分析最大高差法(或模糊数学法)、均值变点分析法和累积和算法与文献[13]数据的相关性，分析得出均值变点分析法与文献[13]中的起伏度最接近，相关系数为0.89，累积和分析算法次之，相关系数为0.757，最大高差法与模糊数学法相关性最差，相关系数为0.39，说明最大高差法与模糊数学法已不适用，这与前人研究结果相一致[10]，均值变点分析法和累积和分析算法是相对比较有效的方法，而且均值变点分析法更适合山西省地形起伏度最佳统计单元的确定。

通过上述分析可以得出：不同方法计算所得的最佳统计单元面积有所不同(表4)，人工作图法易受主观因素的影响，很难确定最佳统计单元，最大高差法和模糊数学法所得的最佳统计单元均为3×3(0.072 9 km2)，均值变点分析法所得的最佳统计单元为11×11(0.980 1 km2)，而累积和分析算法的最佳统计单元为14×14(1.587 6 km2)；由表5得知，均值变点分析法和累积和分析算法所得的最佳统计单元下的起伏度与文献[13]中的起伏度较接近，以小起伏为主，其面积占研究区总面积的1/2左右，且这两种方法划分的基本地貌类型所占面积比例与实际较相符，山地丘陵占80%之余；均值变点分析法是计算山西省地形起伏度最佳统计单元的一种最为理想的方法。

4 结论与讨论

(1) 基于90 m分辨率SRTM DEM数据，利用邻域分析法在不同大小窗口下(2×2，3×3，4×4，…，30×30)对山西省地形起伏度进行了提取(表1)，平均起伏度值起初随邻域面积的增大逐渐增大，但邻域面积达到一定的阈值时，其变化趋势逐渐平缓，采用幂模型对平均起伏度值和邻域面积的拟合效果最好，对数模型次之；

(2) 运用人工作图法、最大高差法、模糊数学法、均值变点分析法和累积和(CUSUM)分析算法分别计算最佳统计单元。结果表明：人工作图法易受主观因素的影响，最大高差法和模糊数学法计算所得最佳统计单元下的地形起伏度分级结果与实际很不相符，运用均值变点分析法与累积和分析算法计算所得的最佳统计单元分别为11×11(0.980 1 km2)和14×14(1.587 6 km2)，最佳统计单元下的地形起伏度分级结果与实际较相符，且使用均值变点分析法和累积和分析算法划分的基本地貌类型与实际较相符，均值变点分析法和累积和分析算法是相对比较有效的方法，其他方法已不适用，这与前人研究结果相一致；

(3) 均值变点分析法的结果与已有研究结果[13]最接近，研究区地形起伏度以小起伏为主，说明该方法是计算最佳统计单元的一种最为理想的方法；

(4) 已有研究[13]所用数据为30 m分辨率的ASTERG DEM数据，与本文所用数据类型及数据分辨率有所不同，而文献[13]使用均值变点分析法计算所得的地形起伏度最佳统计单元为17×17(0.260 1 km2)，可见对于同一地区而言，不同DEM数据类型、不同分辨率DEM，所使用的最佳统计单元大小是不一样的。

本研究分析研究出了适合山西省地形起伏度最佳统计单元的最为理想的计算方法，对均值变点分析法和累积和分析算法计算所得最佳统计单元大小存在差异还需进一步的研究。

参考文献:

[1] 刘育成,赵廷宁.基于变点分析法提取废弃采石场地形起伏度的方法[J].水土保持研究,2016,23(3):269-273.

[2] 涂汉明,刘振东.中国地势起伏度研究[J].测绘学报,1991,20(4):311-319.

[3] 刘新华,杨勤科,汤国安.中国地形起伏度的提取及在水土流失定量评价中的应用[J].水土保持通报,2001,21(1):57-59,62.

[4] 马晓微,杨勤科,刘宝元.基于GIS的中国潜在水土流失评价研究[J].水土保持学报,2002,16(4):49-53.

[5] 闫满存,李华梅,王光谦.广东沿海陆地地质环境质量定量评价研究[J].工程地质学报,2000,8(4):416-425.

[6] 祁元,刘勇,杨正华,等.基于GIS的兰州滑坡与泥石流灾害危险性分析[J].冰川冻土,2012,34(1):96-104.

[7] 刘康,欧阳志云,王效科,等.甘肃省生态环境敏感性评价及其空间分布[J].生态学报,2003,23(12):2711-2718.

[8] 凡非得,王克林,熊鹰,等.西南喀斯特区域水土流失敏感性评价及其空间分异特征[J].生态学报,2011,31(21):6353-6362.

[9] 徐天献,王玉宽,傅斌.汶川地震重灾区土壤侵蚀敏感性评价[J].中国水土保持,2011,18(1):39-42.

[10] 赵斌滨,程永锋,丁士君,等.基于SRTM-DEM的我国地势起伏度统计单元研究[J].水利学报,2015,46(S1):284-290.

[11] 陈学兄,常庆瑞,郭碧云,等.基于SRTM DEM数据的中国地形起伏度分析研究[J].应用基础与工程科学学报,2013,21(4):670-678.

[12] 王志恒,胡卓玮,赵文吉,等.应用累积和分析算法的地形起伏度最佳统计单元确定[J].测绘科学,2014,39(6):59-64.

[13] 陈学兄,毕如田,刘正春,等.基于ASTER GDEM数据的山西地形起伏度分析研究[J].山西农业大学学报:自然科学版,2016,36(6):417-421.

[14] 中国科学院地理研究所.中国1:1000000地貌图制图规范(试行)[M].北京:科学出版社,1987:33-34.

[15] 李炳元,潘保田,韩嘉福.中国陆地基本地貌类型及其划分指标探讨[J].第四纪研究,2008,28(4):535-543.