土地利用/土地覆被信息提取在研究区域整体生态环境状况方面起着重要作用^[1]。遥感技术以其快速、同步监测等特点, 成为提取土地利用/植被覆盖信息高效的技术手段。近年来, 随着遥感技术的快速发展, 将多源遥感数据应用于植被覆盖类型和土地利用现状变化监测研究, 以及植被物候监测已经非常普遍。基于决策树法, 结合光谱特征和物候信息, 对植被覆盖类型提取和土地资源监测的技术已经非常成熟^[2—9]。如韩东等基于无人机和决策树算法对浑善达克沙地的榆树疏林草原木本和草本植被类型实现了自动划分和植被覆盖度的估算^[10]。姚永慧等为提取太白山各级植被类型, 通过资源三号、高分二号的多源多时相高分辨率遥感数据与地形约束因子、地面调查数据和以往研究资料相结合, 实现了大中比例尺植被类型图的编制^[11—13]。WANG等采用支持向量机分类方法, 利用北京1号卫星数据提取了鄱阳湖的NDVI指数和NDWI指数, 对鄱阳湖湿地植被功能型进行了分类^[14]。贾明明等基于面向对象分类方法, 通过获取MODIS时序数据的物候特征并提取了HJ-1 CDD的光谱和纹理特征, 将三者结合起来, 实现了对森林植被类型的提取, 结果表明加入物候特征可大幅提高模型的分类精度^[15]。随机森林分类方法(Random forest classification)是由美国科学家Breiman提出的一种基于决策树组合的机器学习新型高效的分类器, 其优越的性能广泛应用于国外诸多领域^[16—26]。然而, 国内学者对于RF模型的研究与应用相对较少。

本研究通过详尽的实地调查, 获取研究区植被覆盖类型多样性方面的基础数据, 建立地形约束因子, 揭示了子午河中下游流域的植被类型、分布格局、群落结构特征, 并基于样方数据和遥感地面调查点数据, 结合影像的光谱特征、地形因子以及纹理特征建立随机森林模型, 对子午河中下游流域的植被覆盖类型和土地利用类型进行分类, 同时, 与传统的SVM模型进行精度对比分析, 评价并分析模型特征变量的重要性, 研究基于多源信息的RF模型在植被覆盖类型和土地利用现状分类上应用的可行性, 基于修正的通用土壤流失方程RUSLE模型, 利用满足要求的分类数据进一步计算研究区土壤侵蚀模数, 并利用ArcGIS软件进行可视化处理, 得到研究区的土壤侵蚀模数分布图, 结合上述土地利用/植被覆盖等信息, 从宏观角度出发对子午河中下游流域进行生态环境质量评价, 量化该地区的生态环境状况。

1 研究区概况

子午河中下游流域位于陕西省汉中市佛坪县、洋县与安康市宁陕县、石泉县四县交界地带, 地理位置为107°46′8″E—108°18′59″E, 33°10′ 38″N—33°41′2″N, 属秦岭南麓, 面积约为2476.9km², 地势上北高南低, 构造形迹以褶皱挤压变形为主, 属相对稳定地区。区内主要河流有子午河、蒲河和椒溪河等；由于冬季受到北方冷气团控制, 夏季受到南方副高压带的影响, 导致气候、地形均呈现差异性变化, 研究区南部地区多年平均降雨量为806mm, 多年平均气温为14.5℃, 多年平均蒸发量为1065.6mm, 而北部地区多年平均降雨量为901mm, 多年平均蒸发量为1213mm, 属于亚热带北缘山地暖温带湿润季风气候。研究区山多地少, 农业基本条件差, 主要农作物为稻谷、小麦、玉米、油料和大豆等^[27]；土地利用类型以草地、农耕地和林地为主；研究区林地以针叶林和混交林为主, 有较少的含常绿阔叶树的落叶阔叶林, 其中针叶林主要为油松林、马尾松林, 混交林主要是油松-栎类混交林, 阔叶林以青冈栎林、栓皮栎林、板栗林、核桃林、槲栎林、化香林等为主^[28—29]。区内海拔一般500—2500m, 植被垂直地带分布较复杂, 自下而上分布有灌草丛带、落叶阔叶林带、针阔混交林带和针叶林带(图 1)。

2 研究方法

本研究基于Landsat 8 OLI_TIRS影像的光谱特征、纹理特征以及地形因子建立随机森林模型, 对研究区土地利用/土地覆被进行分类, 与传统的SVM分类方法进行精度对比分析, 评价并分析模型特征变量的重要性。同时, 基于修正的通用土壤流失方程RUSLE模型, 利用满足要求的分类数据进一步计算研究区土壤侵蚀模数, 绘制该区域的土壤侵蚀分布图, 结合上述土地利用/土地覆被信息, 从宏观上量化子午河中下游流域的生态环境状况, 具体流程如图 2所示。

2.1 遥感数据来源及预处理

本文以Landsat 8 OLI_TIRS影像(15m分辨率)为基础数据源, 下载自地理空间数据云数据共享网站(http://www.gscloud.cn/search)。选取2015—2019年覆盖整个生长季(5—9月)共14景影像, 所选影像云量均少于5%。首先, 利用遥感图像处理平台对遥感图像进行辐射定标处理, 消除传感器本身产生的误差；然后利用ENVI 5.3软件的MODTRAN模型对遥感影像进行Flassh大气校正, 将其保存为TIFF格式, 再利用ArcGIS的数据处理模块对遥感影像进行镶嵌拼接与裁剪。为充分利用landsat8 OLI_TIRS影像丰富的波段光谱信息和高分辨率, 还需利用ENVI 5.3软件将30m分辨率的多光谱影像与15m分辨率的全色波段影像进行图像融合处理, 最终得到覆盖整个研究区既具有多光谱信息且分辨率为15m的完整遥感影像图^[34]。其他辅助数据包括DEM(digital elevation model)高程数据、降雨侵蚀力、土壤可蚀性、水资源量、氨氮、化学需氧量、氮氧化物、二氧化硫、烟尘等数据, 其中降雨侵蚀力和土壤可蚀性来源于《生产建设项目土壤流失量测算导则》(SL 773-2018), 水资源量来源于陕西省水文水资源勘测中心出版的《2019年陕西省水资源公告》, 氨氮、化学需氧量等污染负荷数据来源于长江水资源保护科学研究所提供的《陕西省引汉济渭工程环境影响报告书》。

此外, 本文还运用到了遥感地面样方调查数据。该数据是研究团队于2020年9月7—11日和2021年4月26—30日在调查区展开实地调查, 通过OMAP结合遥感影像与GPS采样, 依据典型性原则, 样方在能够代表整个样地植被、地形及土壤等特征的区域随机布设, 获取了研究区各海拔梯度上共350个植被类型调查点和626个土地类型调查点(图 3)。每个调查点设置1个样方, 样方大小为10m×10m。植被类型调查点涵盖研究区各海拔梯度上主要的4个植被型组, 即灌草丛、落叶阔叶林、针阔混交林和针叶林；土地利用类型调查点覆盖了研究区6种主要的土地利用类型, 即植被、居民地、裸地、建筑物、水域和农田。

2.2 制图方法

首先根据30m空间分辨率的DEM数据利用ArcGIS软件的3D Analyst模块, 计算研究区的坡向并提取相应垂直带的海拔高度数据, 结合研究区植被的垂直带谱信息生成研究区植被类型图的地形约束因子(图 4), 分析并提取区内针叶林、阔叶林、针阔混交林、灌丛草甸等植被型组。然后基于实地样点调查数据和遥感地面调查数据构建RF模型和SVM模型, 对比两种分类结果的总体精度(Overall Accuracy)和Kappa系数, 筛选得到符合精度要求的研究区域的植被和土地利用现状初步分类图。在地形约束因子的辅助下, 通过目视解译对分类结果进行初步校正, 确定植被型组。然后结合专家校正以提高分类精度, 对最终分类结果进行精度评价, 并分析和评价光谱特征变量对模型的重要性和适用性。随机森林模型是一种基于决策树组合的机器学习新型高效的分类器, 通过Weka软件中的随机森林模块和RStudio软件进行预测和计算。利用RF模型进行分类通常只需要设定2个参数, 分别是生长树的数量(ntree)和特征变量的数量(mtry)；本文通过实验发现, 进行植被覆盖和土地利用类型分类, 当ntree≥500, mtry分别为5和4时, 分类结果的袋外误差值(OOB)最小, 其值分别为19.69%和11.4%。SVM模型是基于广义肖像算法的一种新型机器学习算法, 利用有限的样本数量, 求出学习样本的最优超平面, 将不同类型的数据分开。

考虑到Landsat 8 OLI_TIRS遥感影像的特点, 充分利用1—7波段的光谱信息；并且考虑到研究区位于秦巴山区, 地形复杂多样, 本研究选用了研究中常用的相关地形因子, 如坡度(slope)、坡向(aspect)以及高程(elevation)。归一化植被指数(NDVI)、比值植被指数(RVI)、土壤调节植被指数(SAVI)等可以反映不同的植被类型, 且利用多波段信息监测植被比单波段更为稳定。此外, 纹理特征也是区分不同地物的重要变量。最终, 共提取出17项分类特征变量。

根据实地样点调查和遥感地面调查, 共得到350个植被样点数据和626个土地样点调查数据, 对两类调查数据分别随机选取其中2/3作为训练样本, 在Wake软件中设置模型参数, 建立RF模型, 另外1/3未被选中的样点用于分类结果的精度评价。本文选用混淆矩阵(Confusion matrix)的方法对分类结果进行精度评价, 以Kappa系数和总体精度(Overall Accuracy)作为评价指标。初步分类结果经专家校正后, 在Arcgis10.4软件平台并进一步修正以提高分类精度。

结合研究区植被和土地利用现状分类结果, 利用美国农业部提出的修正的通用土壤流失方程RUSLE模型^[36], 基于Arcgis10.4软件平台Spatial Analyst工具模块中的栅格计算器计算出研究区的土壤侵蚀模数。其表达式为：

(1)

式中, A为土壤侵蚀模数(t hm^-2a^-1)；R为降雨和径流因子(MJmm hm^-2h^-1a^-1)；K为土壤可蚀性因子(t hm²MJ^-1hmm^-1 hm^-2)；L为坡长因子；S为坡度因子；C为植被覆盖因子；P为水土保持措施因子, C与P为无量纲单位。f为使A的单位(因R、K使用美习用单位)转换为国际通用单位的常数, 其值为224.2^[36]。

2.3 生态环境状况评价

根据《生态环境评价技术规范》(HJ192—2015), 生态环境状况评价利用一个综合指数(生态环境状况指数, EI)反映区域生态环境整体的状况, 指标体系包括生物丰度指数、植被覆盖指数、水网密度指数、土地胁迫指数、污染负荷指数和一个环境限制指数, 各项评价指标的权重如表 1所示。

生态环境状况计算方法如下：生态环境状况指数(EI)=0.35×生物丰度指数+0.25×植被覆盖指数+0.15×水网密度指数+0.15×(100-土地胁迫指数)+0.1×(100-污染负荷指数)+环境限制指数, 其中各指标计算方法如表 2所示：

根据生态环境状况指数, 将生态环境分为5级, 即优、良、一般、较差、和差, 如表 3所示。

3 实验结果与分析 3.1 分类结果与精度检验

本文分别采用RF模型和SVM模型对秦岭南坡子午河中下游流域的各种地物类型进行识别, 初步得到研究区土地利用/土地覆被类型图；然后在地形约束因子的辅助下, 通过目视解译对分类结果进行初步校正, 确定植被型组；最后结合专家校正得到研究区的土地利用现状和植被类型结果图(图 5, 6)。随机森林的植被分类结果显示, 研究区59.21%的植被类型为针叶林, 25.82%的植被为针阔混交林, 落叶阔叶林和灌草丛仅占14.97%；从土地利用现状分类结果来看, 研究区大部分区域是林地, 占70.05%, 建筑所占的百分比为6.26%, 水域占10.37%。支持向量机模型的植被分类结果表明, 针叶林占林地总面积的61.33%, 针阔混交林占29.22%, 灌草丛和落叶阔叶林分别占5.06和4.39%;土地利用的分类结果表明, 林地面积占研究区总面积的87.59%, 建筑占1.16%, 水域面积占8.99%。

为了定量评价两种分类方法的分类精度, 本文分别计算模型的总体精度(Overall accuracy)、Kappa系数以及林地的使用者精度(User accuracy)和生产者精度(Producer′s accuracy)(表 4, 5)。与传统的SVM方法相比, 利用RF模型对植被覆盖类型进行分类的总体精度为80.31%, 提高了4.37%；Kappa系数为0.734, 提高了4.2%；对于土地利用分类的总体精度为90%, 提高了5.4%, Kappa系数为0.864, 提高了5.9%。通过利用样点对分类结果图进行精度验证, 两种分类方法对于针叶林和灌草丛的识别度较高, 对针阔混交林和落叶阔叶林的分类精度相对较低, 原因在于针阔混交林的构成较为复杂, 与落叶阔叶林的光谱特征易被混淆。土地利用现状分类结果的精度较高, 但两种分类方法对植被和农田的识别度较低, 原因是研究区的主要作物水稻、蔬菜、玉米, 其光谱特征与植被很难区分。对于模型的生产者精度和使用者精度, 随机森林分类法都比支持向量机分类法有了明显提升, 由此可见, 基于光谱特征变量和地形因子以及纹理特征的随机森林分类法, 比支持向量机分类方法在各方面精度都有了明显的提高, 随机森林相对于支持向量机分类法有更好的分类效果。

3.2 特征变量的重要性

对于多源数据的分类, 在各种不同类型特征变量数量较多时, 需要了解这些特征变量对分类结果的贡献, 以便挑选出较重要的特征变量。本文选取的特征变量及其重要性(图 7)。总体上看, 在进行植被分类时, 地形因子中高程因子的贡献率占据主导地位, 说明高程对于研究区的植被类型的影响极高, 与实地调查和文献资料中的结果完全一致, 光谱特征变量中的土壤调节植被指数(SAVI)贡献率较大；进行土地利用现状分类时, 对分类结果起决定性作用的因子是归一化植被指数(NDVI), 原因是NDVI是指示植被生长状态和植被覆盖程度的指示因子, 其次就是纹理特征中的均值、水体归一化植被指数(NDWI)以及高程因子对土地利用现状分类的贡献率较大。综上分析, 在利用随机森林模型分类时应尽量选择使样本差异性较大的特征变量, 要结合研究区的地形特点进行选取。

3.3 模型结果应用 3.3.1 土壤侵蚀模数计算

基于修正的通用土壤流失方程RUSLE模型, 结合上述符合要求的土地利用/土地覆被分类结果, 利用ArcGIS 10.4软件计算研究区的土壤侵蚀模数(图 8)。计算结果显示, 研究区91.56%的区域为微度侵蚀, 土壤侵蚀模数均小于200t km^-2 a^-1, 仅有水源区附近侵蚀相对较严重, 其中输水沿线土壤侵蚀模数为1000—1500t km^-2 a^-1, 属于轻度侵蚀, 水源区土壤侵蚀模数较大, 大部分区域土壤侵蚀模数大于2000t km^-2 a^-1, 主要侵蚀方式为水流侵蚀, 风力侵蚀的可能性较小。只有局部区域由于河道变窄, 水流流速增大导致土壤侵蚀模数为2500—5000t km^-2 a^-1, 属于中度侵蚀。

3.3.2 生态环境状况评价

根据《生态环境评价技术规范》(HJ192—2015), 结合土地利用现状、植被覆盖类型和土壤侵蚀模数以及地面生态环境等相关参数, 计算出研究区的生态环境状况指数, 从而对研究区的生态环境状况进行定量化评价, 得到研究区的生态环境状况分布(图 9)。可以看出, 研究区平均生态环境状况指数(EI)为87.12, 其中74.31%的区域生态环境状况指数都大于75, 生态环境状况为优；仅有15.69%的区域生态环境状况指数为55—75, 生态环境状况为良, 这些区域多分布在水源区附近以及洋县黄金峡水利枢纽附近。经查阅文献和实地调研, 总结出水源区附近和洋县下游与周边的生态环境状况产生差异性的原因主要有以下两点：第一, 水源区附近土壤受到水流侵蚀, 水土流失相对较严重；第二, 由于水库蓄水导致水源区附近植被遭到破坏, 生物丰度下降, 另外由于黄金峡水利枢纽正处于施工期, 有大量污染物排放, 会使下游水质变差, 且洋县上游存在生活污水排放, 致使片区生态环境状况略微受到影响。

4 结论与讨论 4.1 结论

(1) 利用Landsat 8 OLI_TIRS遥感影像对土地利用/土地覆被信息提取的实验中, 随机森林(random forest)分类法能够很好地适用于秦岭南坡子午河中下游流域土地利用/土地覆被信息提取, 分类结果表明研究区59.21%的植被类型为针叶林, 25.82%的植被为针阔混交林, 落叶阔叶林和灌草丛仅占14.97%；从土地利用现状分类结果来看, 研究区大部分区域是林地, 占70.05%, 建筑所占的百分比为6.26%, 水域占10.37%。

(2) 随机森林模型比支持向量机法对土地利用/土地覆被的信息提取更具优势, 分类速度更快、精度更高, 是一种有效的分类方法。模型精度验证结果显示, 随机森林模型对植被覆盖分类的总体精度比支持向量机模型提高了4.37%；Kappa系数提高了4.2%；对于土地利用分类的总体精度前者比后者提高了5.4%, Kappa系数提高了5.9%。

(3) 研究区91.56%的区域土壤侵蚀强度为微度侵蚀, 输水沿线土壤侵蚀强度较高, 属于轻度侵蚀, 主要侵蚀方式为水流侵蚀, 风力侵蚀的可能性较小；局部河道由于断面变窄, 水流流速增大, 导致土壤侵蚀较为严重, 属于中度侵蚀。

(4) 研究区内整体生态环境状况为优, 平均生态环境状况指数为87.12, 仅有15.69%区域生态环境状况为良, 各个区域生态环境状况评价结果与实际调研结果大致相符, 说明随机森林模型的分类结果能够应用于子午河中下游流域的生态环境状况评价。

4.2 讨论

(1) 就土地利用/土地覆被信息提取方法而言, 以往类似的研究多基于朴素贝叶斯、CART决策树、支持向量机和面向对象模型等, 如孙亚楠、刘睿和张露洋等选取不同的特征变量基于决策树、支持向量机和面向对象模型分别对河套灌区永济灌域、新疆石河子市和辽河平原地区的土地利用类型进行了分类, 分类精度相对于传统的像元相似度分析方法有了较大的提升^[38—40]。然而, 这些方法存在耗时长且容易发生过拟合现象。本文以光谱特征、地形特征和纹理特征作为特征变量, 构建的随机森林模型和支持向量机模型对秦岭南坡子午河中下游流域土地利用/土地覆被信息进行了提取, 将两者的分类精度进行对比分析, 结果表明RF模型相对于SVM模型分类速度更快, 且模型的总体精度、生产者精度和使用者精度均优于SVM模型。

(2) 本文的局限性主要体现在以下几个方面：首先, 本文的实验研究仍采用监督分类结合目视解译的方法进行的, 在利用随机森林模型进行监督分类时采用的是Landsat 8 OLI_TIRS遥感影像, 受到其分辨率和波段数量的限制, 并且子午河中下游流域地形复杂, 只是简单地结合光谱特征、纹理特征和地形因子等特征变量构建随机森林模型也具有一定的局限性^[41]。其次, 本文虽然分析了特征变量对随机森林模型的重要性, 但没有具体分析和评价各特征变量对随机森林模型分类精度会产生多大的影响^[42]。因此, 后续的研究中, 需要使用高分辨率的高光谱影像, 并结合研究区的水文气象因子以及影像的几何特征等特征变量构建更加精确的随机森林模型, 并对模型的特征参数进行更加深入的研究, 实现更高精度的遥感信息提取, 从而使生态环境评价更加切合实际, 为有关部门提供更加准确的数据参考。