基于随机森林模型的云南元阳梯田地形因子分析

赖自力， 向　杰， 陈建平,2， 王文杰， 田夏一， 胡　彬

(1.中国地质大学(北京)地球科学与资源学院，北京100083； 2.北京市国土资源信息研究开发重点实验室，北京100083)



基于随机森林模型的云南元阳梯田地形因子分析

赖自力1， 向杰1， 陈建平1,2， 王文杰1， 田夏一1， 胡彬1

(1.中国地质大学(北京)地球科学与资源学院，北京100083； 2.北京市国土资源信息研究开发重点实验室，北京100083)

为探究地形因子对梯田空间分布的影响，以云南元阳梯田坝达流域为研究区，运用GeoEye-1遥感卫星数据(2010年，1 m分辨率)和数字高程模型(20 m分辨率)，提取梯田空间分布信息以及海拔高程、地形坡度、地形坡向、地面曲率、高程变异系数、地形起伏度、地面粗糙度、地表切割度、到水系距离9个地形因子。采用随机森林(Random Forest)方法建模，结合ROC曲线和AUC值对模型进行精度评价，根据随机置换残差均方减少量和因子作用曲线，探讨地形因子对梯田分布的影响规律。研究结果表明：该流域内梯田总面积1 158.395 hm2，模型AUC值为0.947，海拔高度因子的随机置换残差均方减少量为388.14，到水系距离因子为199.77，地面坡向因子为80.26，三者占总值的73.45%。由此可见，元阳梯田的空间分布主要受海拔高度、到水系距离、地面坡向3类地形因子的影响，其因子曲线表明梯田分布与地形因子间呈非线性关系。

地形因子；元阳梯田；随机森林；因子评价；云南

0　引　言

地形因子是山地梯田农业发展的重要影响因素之一，对地面物质流动与能量交换有着直接或间接的影响，使得地类的分布呈现一定的规律特征(姚敏等，2006)。针对地类所蕴含的内在规律特征，结合数字地形模型，研究地类分布与地形特征的耦合关系是地学景观研究的一个重要方面。赵松乔(1958)从地形因素(坡度、坡向、高度)出发，对内蒙古自治区土地利用进行了综合评价，分析了地质地貌条件在土地资源调查中的作用；宋乃平等(1993)以宁夏为例，研究了地貌形态与土地利用之间的关系；朱翔(2000)运用GIS技术，在杨系流域采用地形分析方法，研究了农业生产与海拔、坡度、坡向的关系；邱扬等(2003)从地类格局和地形特征入手，探讨了黄土丘陵小流域内土地利用时空分布变化与地形因子间的耦合关系；贾宁凤等(2007)研究了砖瓦窑沟流域土地利用与地形因子的相关性；哈凯等(2015)运用典型相关分析(Canonical Correlation Analysis, CCA)方法，分析了河北怀来县土地类型与地形环境的相关性。可见，基于数字地形模型，运用GIS空间分析技术建立地类分布与地形特征的耦合模型，分析地形因子的影响特性，进而掌握地类的分布规律，能够为优化地类的空间配置提供理论基础。

1　研究区概况与研究数据

1.1研究区概况

流域是天然形成的地理单元，其内部的物质分布及能量传递相对独立(刘世斌，2013)。坝达梯田流域位于云南省东南部、元阳县中部，是元阳梯田文化景观的核心区之一。流域东西长8 km，南北宽8.8 km，总面积47.53 km2。区域内最高海拔2 260 m，最低海拔940 m，高差1 320 m，属于亚热带山地季风气候，年平均气温16.4 ℃，全年相对湿度为85%，年平均降雨量在1 397.6 mm左右(胡文英，2009)。其森林资源丰富，梯田广泛分布，地理位置如图1所示。

图1　 研究区地理位置示意图Fig.1　Sketch map showing the location of the study area

1.2研究数据

本次研究所用数据包括2010年GeoEye-1遥感卫星数据(空间分辨率为1 m，云量为1.2%)、元阳县区域规划图(云南数字乡村网：http://www.ynszxc.gov.cn/S1/default.shtml)和1∶5万地形图(云南省测绘资料，北京54坐标系，1956年黄海高程系)。以地形图为基准数据，对GeoEye-1遥感数据进行几何校正，精度在0.5个像元内。

2　随机森林算法

随机森林算法由Breiman(2001)提出，目前在生物医学、经济管理、地学等众多领域得到了广泛应用。在生物医学研究中，Parkhurst等(2005)基于随机森林算法，分析了沙滩细菌与影响变量的关系；在经济管理研究中，林成德等(2007)运用随机森林算法来筛选企业评估指标，并建立了企业信用评估模型；在地学研究中，李亭等(2014)以深圳为例，分析了不同因子对滑坡的影响规律。

随机森林算法通过自助法(Bootstrap)对样本数据进行抽样，由随机向量(X，Y)构建组合模型(其中X为输入向量，Y为输出向量)，然后将决策树计算结果的平均值作为最终结果，其数学表达式为：

(1)

式(1)中，avk(*)表示取平均值，I(*)为示性函数。由于使用Bootstrap采样，则样本中每个数据未被抽取的概率为(1-1/N)N，当N足够大，收敛于1/e≈0.368，表明接近37%的样本不会用于单次模型构建，由此定义余量函数：

(2)

余量函数可衡量模型可靠度，函数值越大，模型可靠度越高，其泛化误差可表示为：

PE*=Pxy(mg(X,Y)<0)

(3)

(X,Y)表示概率空间，根据幸钦定理，当决策单元增加时，所有序列Θk和PE*几乎处处收敛：

Pxy(PΘ(h(X,θ)=Y)-maxPΘ(h(X,θ)=j)<0

(4)

结果表明，随机森林并不会出现过度拟合现象，处理大数据时稳定性强，便于计算因子非线性作用。

3　地形特征提取与分析

地形是基本的地理组成要素，制约着地球表面物质的分配与能量转移，影响着植被的发育和人类活动(Band,1986)。基于DEM的地形信息，运用严密的数学推导，将地形特征定量化。

根据反映地表信息的不同，将地形特征分为3类：微观地形特征、宏观地形特征、相关地形特征(程维明等，2009；Hutchinson et al., 2000)。微观地形特征包括海拔高程、地形坡度、地形坡向、地面曲率，宏观地形特征包括高程变异系数、地形起伏度、地面粗糙度、地表切割度，相关地形特征包括到水系的距离。

3.1海拔高程

海拔高程(Altitude)对植物的垂直性分带具有一定的影响，研究区最低海拔为940 m，最高海拔为2 260 m，总体呈南高西低。

本次研究使用GeoEye-1遥感卫星数据，经目视解译，共提取流域内梯田1 158.395 hm2。由于梯田大小不一，形态相差较大，仅用梯田中心点信息表示梯田信息，对于大型梯田会造成信息的缺失。因此，根据相同面积内采样数据相等，使用全局采样法，获取流域内7 227个样本数据，其中梯田数据为2 288个，非梯田(除梯田以外的地类)数据4 939个。其海拔高程直方图如图2所示，梯田主要分布在海拔为1 400～1 800 m区域，而非梯田地类在整个研究区均有分布。相对梯田而言，非梯田地类在1 800～2 100 m的范围内较为集中。

图2　研究区海拔高程直方图Fig.2　Altitude histogram of the study area

3.2地形坡度

地形坡度(Slope)为地面某点切平面与过该点水平面的夹角，反映地表面倾斜程度，其函数表达式为：

(5)

式(5)中，s为坡度，p为X方向上的高程变化率，q为Y方向上的高程变化率(贾敦新等，2009)。采用3×3的窗口提取地形坡度信息(图3)，研究区坡度范围在0°～56.8°之间，流域北部和中部坡度较大，西南坡度较缓。梯田集中分布在坡度为10°～20°的范围内，其直方图呈正态分布(图4)。

3.3地形坡向

地形坡向(Aspect)为地面切平面法线与正北方向在顺时针方向上的夹角，其函数表达式为：

(6)

式(6)中，Aspect为坡向，p为X方向上的高程变化率，q为Y方向上的高程变化率(李天文等，2004)。坡向不同则太阳辐射强度不同，导致植被发育程度不同。

一般而言，湿润地区，阳坡植被发育完善，阴坡植被发育相对较差。流域内地形坡向在0°～360°均有分布，即流域内包含了所有方位的地形坡向(图3)。在290°～360°的坡向范围内，即西北—北坡向的地区梯田是非梯田地类的1.6倍(图4)。

3.4地面曲率

地面曲率(SurfaceCurvature,SC)代表地面某点的弯曲变化程度，其函数表达式为：

(7)

式(7)中，SC为地面曲率，s为地形坡度，p为X方向上的高程变化率，q为Y方向上的高程变化率。研究区地面曲率范围在-10.58～10.64之间(图3)，梯田相对非梯田更集中于[-2,2]之间(图4)。

3.5高程变异系数

高程变异系数(Variance Coefficient in Eleva-tion,VCE)反映高程的相对变化，其函数表达式为：

(8)

式(8)中，VCE为高程变异系数，S为高程标准差，Z为高程值。使用3×3的窗口提取研究区高程变异系数(图3d)，其范围在[0.96～1.03]度之间。

3.6地形起伏度

地形起伏度(Undulating Terrain,UT)也称局部地势差，反映地表的相对变化差，函数表达式为：

UT=Hmax-Hmin

(9)

式(9)中，UT为地形起伏度，Hmax为地表分析窗口内的高程最大值，Hmin为分析窗口内的高程最小值。依据坝达梯田流域的地理特征，需用Python进行编程，选取3×3的分析窗口提取研究区地形起伏度。流域内地形起伏度在0～75.37 m之间，其西南起伏较小，东北起伏较大(图5a)。梯田集中分布在起伏度为10～20 m的区间内，相较于非梯田而言更加集中(图4e)。

图3　坝达梯田流域地形因子分布图(1)(a) 地形坡度；(b) 地形坡向；(c) 地面曲率；(d) 高程变异系数Fig.3　Distribution of topographic factors in the study area (1)(a) slope; (b) aspect; (c) surface curvature; (d) variation coefficient of elevation

图4　坝达梯田流域地形因子直方图(a) 地形坡度；(b) 地形坡向；(c) 地面曲率；(d) 高程变异系数；(e)地形起伏度；(f) 地面粗糙度；(g)地表切割度；(h)到水系距离Fig.4　Histograms of topographic factors(a) slope; (b) aspect; (c) surface curvature; (d) elevation variation coefficient; (e) undulating terrain; (f) surface roughness; (g) surface cutting degree; (h) distance to river system

3.7地面粗糙度

地面粗糙度(TerrainRoughness,TR)反映了地表面的起伏以及侵蚀状态，定义为地表面积与其水平面投影面积之比，其函数表达式为：

(10)

式(10)中，TR为地面粗糙度，Slope表示地形坡度(以弧度计算)。研究区地面粗糙度的范围在1～1.83 之间(图5b)，梯田主要分布于1～1.06之间。

3.8地表切割度

地表切割度(TerrainDissectedness,TD)是地面一点领域的平均高程与最小高程之差，反映地表的平均下切度，其函数表达式为：

TD=Hmean-Hmin

(11)

式(11)中，TD为地表切割度，Hmean为地表分析窗口内的高程最大值，Hmin为分析窗口内的高程最小值。研究区地表切割度在0～35.79m之间(图5)，梯田主要分布于5～10m之间(图4)。

3.9到水系的距离

到水系的距离是指地面一点到水系的最短距离。基于数字高程模型，运用水文分析，提取流域内的水系分布(图1)。通过GIS空间分析，求得各点到水系的距离(图5d)。研究区到水系的距离在0～1 100m之间，梯田的数量随着到水系距离的增加而急剧降低(图4h)。

图5　坝达梯田流域地形因子分布图(2)(a) 地形起伏度；(b) 地面粗糙度；(c) 地表切割度；(d) 到水系距离Fig.5　Distribution of topographic factors in the study area (2)(a) undulating terrain; (b) surface roughness; (c) surface cutting degree; (d) distance to river system

4　梯田分布与地形因子分析

4.1建模与精度评价

本次研究使用R语言实现随机森林算法，在树状结构生成时，需要由节点变量数(mtry)对树形生长进行控制。由于随机森林采用的是Bootstrap采样，因此运用10次10折交叉验证进行节点变量数的训练，用均方根误差(RMSE)进行误差评定。10次的训练结果如图6所示，当mtry为2、5、9时，交叉验证RMSE中值分别为0.294 9、0.291 1、0.293 3，即最优mtry值为5。

图6　10次10折交叉验证RMSE箱线图Fig.6　Boxplot graph of ten times of 10-fold cross-validation

接收者操作特征曲线(ReceiverOperatingCharacteristicCurve，ROC曲线)是一种基于混淆矩阵的误差可视化工具，横轴为特异度(Specificity)，纵轴为灵敏度(Sensitivity)，常运用于模型效果评价。其曲线下面积(AreaUnderCurve，AUC)为模型精度评价的定量指标，取值范围为0～1，值越大表明模型精度越高。模型ROC曲线如图7所示，其AUC值为0.947，模型精度高，适用于梯田地形因子分析。

图7　随机森林模型ROC曲线图Fig.7　ROC curve based on random forest model

4.2地形因子重要性分析

各个地形因子对梯田分布的影响程度各有不同，对地形因子的重要性分析可以为梯田管理提供指导。运用随机森林法，基于随机置换残差均方减少量，可进行变量重要性(VIM)量化。过程如下。

(1) 运用Bootstrap采样，对每一组样本进行模型的构建，计算各组对应的袋外数据(OOB)，得到b组袋外误差的均方根误差向量MSEi，其中i的取值范围为1～b。

(2) 将变量X在袋外数据OOB中进行随机的置换，重复步骤(1)，得到随机置换后的均方根误差矩阵MSEpi，其中i的取值范围为1～b，p为变量个数。

(3) 用向量MSEi与矩阵MSEpi对应的行向量相减，平均后再除以标准误差，则得到变量X的重要性评分。

运用随机森林法，计算随机置换残差均方减少量如图8所示。海拔高度的随机置换残差均方减少量的值为388.143，到水系距离(199.772 4)，地面坡向(80.261 33)。可见，研究区内的梯田分布受地形特征的影响各有不同，主要是受海拔高度的影响，其次是到水系距离，再次为地面坡向。

图8　地形因子重要性图Fig.8　Importance of different topographic factors

基于随机森林模型绘制主要因子影响曲线图(图9)。海拔高程在1 000～1 600m的区间内，海拔高程因子对梯田的作用具有相似的促进作用；但当海拔高程超过1 600m以后，随着海拔高程的增加，海拔高程对梯田的影响程度逐渐降低；海拔高程超过1 800m以后，对梯田的分布则产生抑制作用。到水系距离因子对梯田的影响则表现为：随着到水系距离的增加，到水系距离因子对梯田分布的贡献作用急剧下降，当到水系距离在450m的范围内，对梯田的分布呈现促进状态。地形坡度的影响主要表现为坡度在290°～360°之间对梯田的分布具有促进作用，这与地学所观测的基本规律相符。可以看出基于随机森林模型的因子重要性判定具有较高的可信性，对地学景观的定量化研究具有一定意义。

图9　典型地形因子影响曲线图Fig.9　Curves showing the influence of typical topographic factors(a) altitude; (b) distance to river system; (c) aspect

4　结　论

以梯田空间分布与地形因子的关系研究为主线，选取海拔高程、地形坡度、地形坡向、地面曲率、高程变异系数、地形起伏度、地面粗糙度、地表切割度、到水系距离9个地形因子，应用随机森林算法，研究地形因子的重要性，分析因子作用特性。

(1) 运用GeoEye-1遥感卫星数据(2010年，1m分辨率)，成功提取了研究区1 158.395hm2梯田的空间分布特征。基于随机森林法构建地形因子分析模型，选用ROC曲线进行精度评价，测得模型精度为0.947，表明随机森林算法适用于梯田地形因子分析。

(2) 梯田的分布与地形因子存在一定的联系，运用随机森林算法，基于随机置换残差均方减少量进行地形因子重要性评判。海拔高度因子的随机置换残差均方减少量为388.14，到水系距离因子为199.77，地面坡向因子为80.26，占随机置换残差均方减少量总值的73.45%。可见梯田分布主要受海拔高程、到水系距离和地形坡向的影响。运用因子作用曲线，可见地形因子对梯田的影响作用呈非线性。

(3) 地形因子作为影响梯田分布的一个重要因子，在对梯田进行结构优化配置时，需根据地形特征进行梯田的合理管理和协调。

程维明，周成虎，柴慧霞，等，2009. 中国陆地地貌基本形态类型定量提取与分析[J]. 地球信息科学学报，11(6):725-736.

胡文英，2009. 元阳哈尼梯田景观格局及其稳定性研究[D]. 昆明：昆明理工大学.

哈凯，丁庆龙，门明新，等，2015. 山地丘陵区土地利用分布及其与地形因子关系：以河北省怀来县为例[J]. 地理研究，34(5):909-921.

贾宁凤，段建南，乔志敏，2007. 土地利用空间分布与地形因子相关性分析方法[J]. 经济地理，27(2):310-312.

贾敦新，汤国安，王春，等，2009.DEM数据误差与地形描述误差对坡度精度的影响[J]. 地球信息科学学报，11(1):43-49.

李天文，刘学军，汤国安，2004. 地形复杂度对坡度坡向的影响[J]. 山地学报，22(3):272-277.

林成德，彭国兰，2007. 随机森林在企业信用评估指标体系确定中的应用[J]. 厦门大学学报(自然科学版)，46(2):199-203.

刘世斌，2013. 流域土地利用功能分区体系研究[D]. 武汉：中国地质大学(武汉).

李亭，田原，邬伦，等，2014. 基于随机森林方法的滑坡灾害危险性区划[J]. 地理与地理信息科学，30(6):25-30.

邱扬，傅伯杰，王军，等，2003. 黄土丘陵小流域土地利用的时空分布及其与地形因子的关系[J]. 自然资源学报，18(1):20-29.

宋乃平，陈忠祥，1993. 地貌与土地利用关系之探讨[J]. 宁夏大学学报(自然科学版)，14(3):27-31.

姚敏，崔保山，2006. 哈尼梯田湿地生态系统的垂直特征[J]. 生态学报，26(7):2115-2124.

赵松乔，1958. 内蒙古自治区的地形条件在土地利用上的评价[J]. 地理学报，13(3):245-256.

朱翔，2000.GIS地形分析技术在小范围土地利用规划中的应用[J]. 云南环境科学， 19(2):11-13.

BANDLE, 1986.Topographicpartitionofwatershedswithdigitalelevationmodels[J].WaterResourcesResearch， 22(1): 15-24.

BREIMANL, 2001.Randomforests[J].MachineLearning， 45(1): 5-32．

HUTCHINSONMF,GALLANTJC, 2000.Digitalelevationmodelsandrepresentationofterrainshape[M]//WILSONJP,GALLANTJC.TerrainAnalysis:TerrainAnalysis:PrinciplesandApplications.NewYork:Wiley, 29-50.

PARKHURSTDF，BRENNERKP，DUFOURAP, 2005.Indicatorbacteriaatfiveswimmingbeaches:analysisusingrandomforests[J].WaterResearch， 39(7): 1354-60.

Analysis on topographic factors of the Yuanyang terrace in Yunnan Province based on random forest model

LAI Zili1, XIANG Jie1, CHEN Jianping1,2, WANG Wenjie1, TIAN Xiayi1, HU Bin1

(1. School of Earth Sciences and Resources, China University of Geosciences(Beijing), Beijing 100083, China; 2. Key Laboratory of Land and Resources Information Research & Development in Beijing, Beijing 100083, China)

To explore the impact of topography on the spatial distribution of the Yuanyang terrace, this study used the GeoEye-1(year of 2010, 20 m resolution) remote sensing data to extract the spatial distribution of this terrace, and acquired the terrain features based on the digital elevation model(20 m resolution). These terrain features include altitude, slope, aspect, surface curvature, variance coefficient of elevation, undulating terrain, surface roughness, surface cutting degree and the distance to the water system. This research used the factor curves to display the relationship between the terrain factors and terrace based on the random forest modelling. In addition, we employed the ROC curves and AUC to evaluate the accuracy of the model. The result shows that the terraced area is 1 158.395 hm2acres in the study area and the value of AUC is 0.947. It shows that the random forest model is suitable for the terrace research. Furthermore, the value of random permutation residual mean square reduce is 388.14 in altitude and 199.77 in the distance to the water system and 80.26 in slope; these three terrain factors account for 73.45% of the total terrain factors. Thus, the spatial distribution of the Yuanyang terrace is mainly affected by the altitude, the distance to the water system, and the slope. Moreover, the terrain factor has a non-linear relationship with the terrace according to the factor curve based on random forest model.

terrain factor; terrace in Yuanyang County; random forest; factor importance evaluation; Yunnan Province

10.3969/j.issn.1674-3636.2016.03.518

2016-06-06；

2016-06-21；编辑：陆李萍

中国地质调查局项目“矿山环境变化自动监测示范研究”(1212011120029)

赖自力(1993—)，男，硕士研究生，地质工程专业，主要研究方向为矿产资源预测，E-mail：1049962421@qq.com

P237； F301.2

A

1674-3636(2016)03-0518-08