沂河流域植被覆盖时空演变及其与径流的关系研究

石智宇, 赵 清, 王雅婷, 张连蓬

(江苏师范大学 地理测绘与城乡规划学院, 江苏 徐州 221116)

植被是联结地表与大气之间物质、能量交换的关键环节，在陆地表面能量交换、水分循环和生物地球化学循环过程中起着至关重要的作用[1]。植被是地理环境的重要组成部分，具有涵养水源、保持水土、调节气候等作用，植被变化与气候、土壤、地形条件、人类活动等要素密切相关。径流是流域水循环的重要环节，径流量的变化会影响水文系统的演化和水资源的开发利用，进而影响生态环境及社会经济的发展[2-3]。植被变化可以通过改变冠层截留、植被蒸腾作用和土壤入渗特性等水文过程，显著地改变地表径流，进而影响流域水资源可利用量[4-7]。因此，基于植被覆盖时空演变过程，逐像元分析流域植被覆盖变化对径流变化的影响，对揭示流域生态环境演变过程和改善流域水资源开发利用有重要意义。

近些年来，国内外学者围绕植被覆盖变化及其对径流的影响方面开展了大量研究。毕早莹等[5]基于遥感植被指数和水文、气象数据，通过基于Budyko假设的弹性系数法，分析了1982—2015年窟野河流域植被变化及其对径流变化的影响，发现窟野河流域内植被得到了极大的恢复，植被变化成为决定径流变化的主导因素，其次为气候变化和人类活动；Kim等[8]选取4个不同冠层的径流区研究植被冠层对地表径流的影响，结果表明枯落物层的蓄水能力和叶片形状影响地表径流的产生，除此之外，土壤水分条件也是影响地表径流的重要因素；支童等[9]以陕北秃尾河流域为研究对象，分析流域植被变化特征，探讨了降水与植被覆盖变化对径流的影响并指出影响径流缩减的物理机制；李天生等[10]基于Budyko理论具体分析了珠江流域中上游地区气候和植被变化对该区域径流的影响，结果发现流域降雨减少是径流减少的主导因素，而植被在研究时段没有显著的变化趋势，其对径流减少的贡献不显著。

遥感技术因能提供反映不同时空尺度上的植被覆盖信息及其动态变化，方便快捷，逐渐成为监测植被覆盖时空演变的主要手段，其中归一化植被指数(Normalized Difference Vegetation Index，NDVI)对绿色植被信息敏感，可以很好地反映地表植被生长状况，其变化在一定程度上能代表地表植被变化[11-12]，是目前应用最为广泛的植被指数之一[13-14]。SCS(Soil Conservation Service)径流模型是美国农业部水土保持局开发的流域水文模型，由于该模型简单易行，所需参数较少，对观测数据要求不高且有效考虑到流域下垫面的特点，被广泛采用[15-17]。使用SCS模型并结合地理信息技术可以获取地表径流的空间分布，从而实现在像元尺度上研究植被覆盖变化与径流的关系[15,18]。

沂河流域地处鲁中南山地丘陵区，是沂蒙山区核心组成部分和重要的水源涵养区，流域内地形起伏大，生态脆弱，水土流失较严重[19]。近年来沂河流域内经济社会持续发展，生态环境发生较大变化，人水关系趋于紧张，但针对沂河流域的研究多集中在水文变化、水土流失、土地利用变化等方面[19-22]，对于植被覆盖变化及其水文效应研究相对较少。在人类活动日益加剧和经济社会快速发展的背景下，沂河流域植被覆盖时空格局是否发生变化？如果沂河流域植被覆盖呈现变化趋势，那么其与地表径流变化的相关关系如何？这些问题都有待解答。因此本文选取沂河流域作为研究区，基于2000—2020年MOD13Q1 NDVI遥感数据，使用趋势分析法，分析沂河流域20年间植被覆盖时空演变特征，在此基础上，结合SCS模型，使用相关性分析法，在像元尺度上探讨植被覆盖变化与径流的相关关系，以期为流域生态环境保护和水资源规划利用提供理论依据，同时为中小型流域植被覆盖变化和径流关系研究提供一定参考。

1 研究区概况

沂河是淮河流域沂沭泗水系中主要河流之一，发源于山东省沂源县西部，流经沂源、沂水、沂南、临沂、蒙阴、平邑、郯城等县市，在江苏省新沂市汇入骆马湖，河道全长333 km，流域总面积为11 820 km2，是鲁南地区和苏北地区重要的山洪河道，有东汶河、蒙河、祊河等诸多支流汇入。

由于沂河下游开辟新沂河、“分沂入沭”等河道整治工程，流域界限被破坏，因此本文选取沂河临沂市以上相对完整的流域作为研究区域。沂河临沂以上河道长223 km，流域面积约9 950 km2，属暖温带半湿润季风气候，年均气温13～14.3℃，流域地势西北高、东南低，每年6—9月为汛期，多年平均年降水量830 mm[21]，土壤类型较多，主要有棕壤、褐土、潮土、砂姜黑土等，流域处于暖温带落叶阔叶林带，主要树种有刺槐、油松、赤松、加拿大杨、旱柳等[19]。降雨是沂河流域径流补给的主要来源且径流年内变化十分剧烈，汛期多山洪暴雨，枯水期径流量很小甚至断流，导致沂河流域成为为易旱易涝地区，对当地经济发展带来不利影响[21,23]。

2 数据与研究方法

2.1 数据来源与处理

本文遥感影像数据采用美国国家航空航天局(NASA)MODIS数据中的MOD13Q1 NDVI数据(http:∥ladsweb.modaps.eosdis.nasa.com)，文件格式为HDF，空间分辨率为250 m，时间分辨率为16 d，数据时间范围为2000年2月至2020年12月，共480景影像。利用MRT(MODIS Reprojection Tool)软件对获取的数据进行投影转换、裁剪和格式转换等处理，获取质量较为可靠的NDVI数据集。采用最大值合成法(Maximum Value Composites，MVC)获取年最大NDVI作为年NDVI值，能较好地代表当年植被生长最好的状况[24]，并可进一步消除云、大气、太阳高度角等的部分干扰[25]。

研究区土地利用数据使用来自于中国科学院资源环境科学与数据中心(http:∥www.resdc.cn)提供的1 km栅格数据，将土地利用类型分为耕地、林地、草地、水域、建设和未利用土地，包括2000年、2010年、2020年三期土地利用数据。

土壤数据来源于联合国粮农组织(Food and Agriculture Organization of the United Nations，FAO)和维也纳国际应用系统研究所(International Institute for Applied System Analysis，IIASA)所构建的世界土壤数据库(Harmonized World Soil Database，HWSD)，采用FAO 90土壤分类系统，包括上层土壤(0—30 cm土层)和下层土壤(30—100 cm土层)的砂粒、黏粒、碎石等百分比含量及USDA土壤质地分类数据。

DEM数据使用地理空间数据云(http:∥www.gscloud.cn/)提供的SRTMSLOPE 90 m分辨率坡度数据产品，经过ArcGIS裁剪得到研究区坡度数据。降水数据来源于中国气象数据网(http:∥data.cma.cn)，包含研究区内沂源、沂水、沂南、蒙阴、平邑、费县和临沂7个雨量站点2000—2020年降水量数据(表1)，经反距离加权插值法(IDW)获取与NDVI数据具有相同投影和空间分辨率的栅格数据。

表1 沂河流域雨量站位置信息及数据时间

2.2 研究方法

2.2.1 趋势分析 本文采用一元线性回归分析逐像元模拟NDVI年际空间变化趋势。一元线性回归分析可以模拟每个栅格的变化趋势，以单个像元时间变化特征反映整个区域时空格局演变[26]，计算公式为：

(1)

式中:slope为像元NDVI回归方程的斜率；i为年序号；n代表时间序列长度；NDVIi代表第i年的NDVI值。若slope>0,则表示NDVI随时间变化呈增加趋势,若slpoe<0,则表示NDVI随时间变化呈减少趋势,拟合斜率绝对值越大,表示变化越明显。

2.2.2 径流模拟 集总式水文模型具有结构简单、参数较少、复杂程度较低等优点，SCS(Soil Conservation Service)水文模型是美国农业部水土保持局于1954年开发的集总式水文模型，实用性强且改进空间大，近年来得到了比较广泛的应用[27]。该模型综合考虑了流域降雨、地形、土壤类型、土地利用方式、前期土壤湿润状况与径流的关系，SCS水文模型基本产流方程[28]为：

(2)

(3)

式中:P为降雨量(mm);Q为径流深(mm);S为最大可能滞留量(mm);CN值是反映降雨前流域特征的一个综合参数,用于描述降雨—径流关系。CN值通常由土壤性质、土地利用方式和降雨事件前期土壤湿润情况等数据来确定[28]。

根据研究区各类型土壤上层与下层土壤质地数据，基于SPAW软件计算研究区上层土壤与下层土壤最小下渗率并选取较小一组作为该土壤类型最小下渗率，参考SCS径流模型的水文土壤分类标准[18](表2)划分流域水文土壤组，在ArcGIS软件下重分类，得到研究区水文土壤组数据(图1)。

表2 SCS模型水文土壤组定义指标

图1 沂河流域水文土壤组

SCS径流模型考虑前期降水对径流的影响，将前期土壤湿润程度根据径流事件发生前5天的降雨总量(即前期降雨指数API)划分为干旱(AMCⅠ)、平均(AMCⅡ)和湿润(AMCⅢ)3种状态，本文研究中假设前期土壤湿润程度为平均状态(AMCⅡ)。

SCS径流模型中不同的下垫面类型对产汇流的影响不同，本文选取沂河流域2000年、2010年、2020年三期土地利用数据，通过GIS软件分析三期土地利用(图2)，得到研究区三期土地利用类型变化，其中2000—2005年、2006—2015年、2016—2020年分别基于2000年、2010年、2020年土地利用数据生成的CN值模拟地表径流深。

结合研究区概况，在SCS模型提供的CN值查算表基础上，参考国内有关学者的研究成果[29]，确定正常条件下的CN值。考虑到研究区地形坡度对CN值的影响，本文采用Huang坡度修正公式[30]对CN值进行修正，公式如下：

(4)

式中:CN2为AMCⅡ条件下CN值；CNa为坡度修正后的CN2值;α为多边形坡度值,用百分比(%)表示。

基于GIS软件，将研究区水文土壤、土地利用数据两层数据交叉，得到包含土壤、土地利用双重信息的数据图层，输入修正后CN值，生成AMCⅡ条件下CN值图层，结合研究区降水数据，利用ArcGIS栅格运算功能得到逐年地表径流深图层。

图2 沂河流域土地利用变化

2.2.3 植被覆盖变化与径流的相关分析 为了分析研究区植被覆盖变化与径流的相关关系，本文采用相关分析方法，基于像元尺度对2000—2020年沂河流域年最大NDVI与径流深进行逐像元相关分析，使用相关系数来反映植被NDVI和地表径流深的相关程度，在p<0.05置信水平对相关性结果进行显著性分析。相关分析计算公式如下：

(5)

3 结果与分析

3.1 NDVI年际变化特征

提取沂河流域年最大NDVI代表当年植被覆盖最好情况，并将其作为该年的NDVI值(图3)，可以看出，2000—2020年沂河流域NDVI整体呈现波动增加趋势，20年间流域年最大NDVI介于0.67～0.79之间，增长速率为1.46%/10 a(R2=0.118)。植被覆盖时间变化大致可以分为3个阶段：(1) 2000—2002年植被覆盖处于较低水平，其中2002年出现20 a间NDVI最低值；(2) 2003—2013年植被覆盖较上一阶段出现较大幅度提升，NDVI值整体保持小幅波动上升趋势；(3) 2014—2020年植被NDVI值较上一阶段出现较大程度波动，2014年出现该阶段NDVI最低值，该年植被覆盖情况几乎退化至2000—2001年水平，但总体上该阶段植被覆盖依然呈现增加趋势。总体而言，2000—2020年沂河流域NDVI呈现增加趋势，流域植被覆盖程度得到了改善，植被生长情况总体向较好方向发展。

图3 2000-2020年沂河流域NDVI年际变化趋势

3.2 NDVI空间分布特征及其变化趋势

在空间分布上，2000—2020年沂河流域年均NDVI空间分布见图4A。可以看出，沂河流域NDVI多年均值范围在0.16～0.91之间，流域内植被覆盖空间分异较大。植被NDVI高值区多集中在流域南部与东部，其中蒙阴县、平邑县与费县三县交界处的蒙山山区NDVI高值最为集中，部分地区NDVI值在0.9以上，该区位于沂蒙山区腹地，生态环境较好，同时建有国家森林公园，植被覆盖度较高。植被NDVI低值区主要集中在流域北部以及流域内各市(县)域城镇建成区，其中流域北部沂山山区和鲁山山区NDVI值除少数地区外整体较南部蒙山山区低，表明流域内沂山山区和鲁山山区多年来植被覆盖情况较蒙山山区差。此外，流域内各市(县)域城镇建成区NDVI值均较低，其中临沂市兰山区低植被覆盖区面积最大，表明城镇用地的增加导致以上区域地表植被覆盖情况较差。

在空间变化趋势上，基于一元线性回归分析在像元尺度上计算2000—2020年沂河流域NDVI空间变化趋势，参考史晓亮等[31]对淮河流域植被NDVI划分的变化趋势等级，将沂河流域NDVI变化趋势划分为3个等级，并统计各个等级的面积及所占比例(图4B，表3)。可以看出，植被NDVI基本不变的区域约占研究区总面积的78.35%，表明流域大部分地区植被覆盖保持相对稳定，不存在明显的增减趋势。植被NDVI呈增加趋势的区域所占比例为15.05%，主要分布在流域北部，以沂源县和蒙阴县境内最为集中，分析认为近些年来对沂蒙山区进行的一系列生态保护与修复措施，使得该区域植被NDVI增加面积持续增加。植被NDVI出现下降的区域面积相对较小，约占研究区总面积的6.60%，主要分布在流域东部和南部，其中流域南部平邑县至兰山区境内出现植被NDVI呈连续片状减少区域，该区域主要集中在祊河流域，结合图2可以看出，该区域20 a间草地面积不断减少，同时耕地和建设用地面积不断增加，分析认为该区域人类活动的增加导致植被NDVI的减少。此外，研究还发现在流域各城镇建成区及其周边区域内NDVI值呈减少趋势，表明城镇化的发展导致NDVI下降，植被覆盖水平降低。

3.3 沂河流域地表径流变化特征

基于GIS软件结合SCS模型，在AMCⅡ条件下生成沂河流域20年逐年地表径流深图层，并依此生成2000—2005年，2006—2010年，2011—2015年，2016—2020年4个时期地表径流深变化图(图5)。结果表明：(1) 流域不同时期径流深均呈现出不同程度的由西北向东南逐渐增加趋势，径流深低值主要集中在流域北部及南部蒙山地区，高值主要集中在流域东部及南部。(2) 不同时期，沂河流域径流深出现较大幅度波动，径流存在明显的丰枯变化，径流深高值区和低值区空间范围和面积均发生了不同程度变化。

和继军等[32]提出，在临界值下，坡面产流量随坡度的增加而增加。随着坡度的增大，水的势能变大，在地表中的流速增大，缩短了入渗时间，径流深减小，从而使地表流量增大，沂河流域北部鲁山山区、沂山山区及南部蒙山山区地形坡度较大，坡面产流能力较强，坡面产流主要向东南部地形较平坦的平原地区汇集。同时结合沂河流域土地利用图发现，20年间流域南部耕地及建设用地面积不断扩大，而北部土地利用结构相对稳定，林地与草地面积比例较南部高。建设用地所造成的不透水面导致地表水渗透比例减小，而耕地由于农业活动频繁，促使土壤压实和结皮，入渗速率和土壤蓄水含量降低[33]，从而导致地表径流深增加。因此结合研究区高程及土地利用结构，分析认为沂河流域地表径流深空间分布趋势主要是地形坡度和土地利用结构共同作用的结果。此外，沂河流域径流变化与降水变化具有显著相关性，流域内季风气候显著，降水年内分配不均且年际变化大，导致沂河流域径流出现较大幅度波动[23]。

图4 2000-2020年沂河流域植被NDVI空间分布及变化趋势分布

表3 2000-2020年沂河流域NDVI变化趋势统计

3.4 植被NDVI与地表径流的相关性分析

基于SCS模型模拟沂河流域20年间逐年地表径流深，逐像元计算20年间年最大NDVI和径流深之间的相关性，并对相关性进行显著性检验(图6)。结果显示，沂河流域植被NDVI和地表径流之间正相关和负相关共存，统计表明，流域内88.79%的区域植被NDVI与径流呈正相关，11.21%的区域植被NDVI与径流呈负相关。由图6A可知，流域北部呈正相关区域面积远大于呈负相关区域，呈负相关的区域则主要集中在流域南部，尤其以兰山区为最。此外，流域内植被NDVI与径流变化相关系数有较大分异，流域正相关系数较大，部分地区可达0.8以上，而在流域内部分呈负相关的区域，相关系数则可接近-0.7。

图5 2000-2020年沂河流域不同时间段地表径流深

对相关性结果采用p<0.05置信水平进行显著性检验(图6B，表4)，可以看出研究区内以不显著正相关为主，研究区内正、负相关达到显著水平的面积分别占统计总面积的34.07%和4.72%，其中呈显著正相关的区域主要集中在流域北部，对比图4B，发现呈显著正相关的区域与NDVI呈改善趋势的区域基本一致，表明随着该区域植被NDVI的增加，地表径流也增加。呈显著负相关的区域主要集中在研究区南部，尤其以临沂市兰山区最为集中，结合图2和图4A，分析认为城镇建设用地的增加导致城市地表不透水面增加，同时造成植被覆盖的减少，导致地表径流的增加。此外，植被NDVI与地表径流呈显著正相关与显著负相关的像元在流域内相互交错分布，表明流域内植被NDVI与径流变化关系在小空间尺度下存在明显差异。

图6 沂河流域植被NDVI与径流变化相关系数及显著性

表4 沂河流域植被NDVI与径流变化相关系数显著性统计

4 讨论与结论

4.1 讨 论

植被覆盖变化受多种因素耦合作用，在空间上呈现差异性，由其引发的水文效应则更加复杂。由于区域气候和地理的差异性，研究尺度与方法的不同，以及植被本身的复杂性，特定流域内植被变化对水文过程的影响结果并不一致。近几十年以来，沂河流域经济社会快速发展，流域内植被覆盖情况的改变及其引发的水文效应势必对现在及未来沂河流域生态环境产生影响。

本文基于遥感影像，在像元尺度分析植被变化及其对径流的影响，实现在较高精度下分析沂河流域植被变化与径流的相关关系。但需要指出的是本文在利用SCS模型模拟地表径流时，只基于前期土壤水分处于平均状态并结合三期土地利用模拟地表径流，且沂河流域降雨主要集中在主汛期的7—9月份，多数降雨为集中性暴雨，因此基于SCS模型的沂河流域径流模拟精度依然尚待改善。此外，本文只是定量描述了植被变化与径流的关系，而受植被变化时空特征影响冠层截留、蒸散发等水文过程，需要进一步利用分布式水文模型，在小尺度背景下，从产汇流过程上具体分析植被空间格局变化对水文过程的影响机制。

沂河是鲁南与苏北地区重要的山洪河道，同时沂河流域也是鲁南山地生态环境的重要组成部分，因此沂河流域植被覆盖变化及其对径流的影响对鲁南与苏北地区经济社会发展意义重大。未来沂河流域经济社会发展应该更加注重土地利用合理规划，加强生态环境保护，而未来沂河流域生态环境变化研究也应该将目光放在自然因素与人为因素相互耦合作用研究上。

4.2 结 论

(1) 2000—2020年沂河流域植被NDVI呈波动增加趋势，增长速率为1.46%/10 a。流域植被NDVI空间分异较大，高值区主要分布在流域南部与东部，低值区主要集中在流域北部以及流域内城镇建成区。流域植被NDVI空间变化趋势整体较为稳定，以基本不变为主，NDVI呈增加趋势的区域主要集中在流域北部，占研究区面积的15.05%，NDVI呈减少趋势的区域主要分布在流域东部和南部及各城镇建成区内，占比为6.60%。

(2) 2000—2020年沂河流域不同时期径流深均呈现出不同程度的由西北向东南逐渐增加趋势，径流深出现较大幅度波动，径流存在明显的丰枯变化。径流深低值主要集中在流域北部及南部蒙山地区，高值主要集中在流域东部及南部，径流深高值区和低值区空间范围和面积均发生了不同程度变化。

(3) 沂河流域植被NDVI和地表径流深的相关关系呈正相关与负相关共存的状态，流域不同空间位置的植被变化与径流变化的关系存在较大分异。88.79%的区域植被NDVI与径流呈正相关，以不显著正相关为主，呈显著正相关的区域主要集中在流域北部，植被NDVI与径流呈负相关仅占比11.21%，其中显著负相关区域主要集中在流域南部兰山区境内。