基于InVEST模型的疏勒河上游产水量时空变化特征分析

魏培洁,吴明辉,贾映兰,高雅月,徐浩杰,刘章文,陈生云1,,6,*

1 中国科学院西北生态环境资源研究院 冰冻圈科学国家重点实验室 疏勒河源冰冻圈与生态环境综合监测研究站,兰州 730000 2 中国科学院大学,北京 100049 3 祁连山国家公园国家长期科研基地,西宁 810000 4 兰州大学草地农业科技学院草地农业生态系统国家重点实验室,兰州 730020 5 中国科学院西北生态环境资源研究院 黑河上游生态水文试验研究站,兰州 730000 6 青海师范大学高原科学与可持续发展研究院 地理科学学院,西宁 810008

生态系统服务是指人类从各种生态系统中获取的所有惠益,由生态系统的支持功能、供给功能、文化功能和调节功能及其相互作用形成[1—2]。生态系统功能是生态系统服务的基础和物质保障,其对人类生存和社会发展至关重要,特别是与水资源相关的功能[3]。随着全球气候变化、地表水资源短缺和水环境问题加剧,空间定量化、可视化评估区域产水量已成为当前研究的热点议题之一[4—8]。产水量是流域特别是干旱地区重要的调节功能,它对工业、农业、渔业、人类消费和娱乐活动至关重要[9—10]。一方面,产水量是人类生存和区域经济可持续发展的重要物质基础；另一方面,产水量与区域自然地理条件和人类活动密切相关[10]。因此,研究区域产水量变化的驱动因素变得十分必要,其可为区域水资源的有效管理和可持续发展提供科学依据[10—11]。

自1970年以来,遥感技术和水文模型发展迅速,越来越多的学者尝试通过模型模拟方法来定量化、可视化和精细化地评估和分析区域产水量,如MIKESHE(MIKE System Hydrological European)模型、TOPMODEL模型(TOPgraphy based hydrological Model)、SWAT模型(Soil and Water Assessment Tool)和InVEST模型等[12—13]。InVEST是一款由美国斯坦福大学、世界自然基金会、大自然保护协会和其他相关机构联合开发用于生态系统服务功能评估的模型方法[14]。InVEST模型产水量模块基于水量平衡原理,充分考虑不同土地利用类型下土壤渗透性的空间差异及地形、地表粗糙程度等因素对径流影响,以栅格为单元定量估算不同土地利用类型的产水量[12—13]。其模型由于参数调整灵活、评估结果空间表达性强,且能够实现时空间和多目标权衡[15],被国外许多国家和地区广泛应用。如：Leh等[16]、Marquèsa等[17]和Redhead等[18]利用InVEST模型分别对西非科特迪瓦、西班牙东北部弗兰科利河流域和英国多条河流的产水量进行了评估,并取得了很好的应用效果。此外,中国学者利用InVEST模型的产水量模块分别对北京山区[19]和贵水流域[20]、陕西省商洛市五大水系流域[21]、三江并流区(云南、四川和西藏的交界处)[22]以及高原湖泊[23]等区域的产水量及其时空分布特征进行了科学评价,模型在各研究区均取得了很好的应用效果。同时,InVEST模型已成功用于青藏高原地区。其中,潘涛等[7]和吕乐婷等[24]分别评估了三江源区1981—2010、1980—2016年产水量的时间变化特征与空间分布状况。王雅等[25]和赵亚茹等[10]分别对位于祁连山区中部的黑河流域、东部的石羊河流域的产水量空间分异特征进行了分析,这些研究对西北地区水资源规划和管理都具有重要意义。

目前,气候变化和土地利用变化是影响产水量的关键因素[26]。其中,气候变化在宏观尺度上可以通过改变流域的降水量和蒸散量来影响产水量[27—28],土地利用可通过改变流域下垫面间接影响产水量[29—30]。一些研究已探究了气候变化和土地变化对流域产水量的影响。例如：Pessacg等[9]利用降水数据集评估了楚布特河流域的产水量,表明降水量与产水量之间存在正相关关系；Stone等[31]重点研究了气候变化对密苏里河流域产水量的影响,证实气候变化对产水量影响显著；Zhang等[32]和Li等[33]通过分析土地利用变化对产水量影响,得出不同土地利用方式的改变对产水量的影响不同。目前,全球气候变化对中国水资源的影响日益加剧[34]。青藏高原作为对全球环境变化特别是气候变化响应最敏感的地区之一,近年来受气候变暖影响,50年来气温显著升高0.3 ℃/(10a),部分区域降水波动较大[35],并且气候暖湿化促进了植被生长,NDVI显著增加[36]。然而,青藏高原地区气候变化和土地利用变化对产水量的相关研究仍然缺乏,尤其植被根系深度和NDVI对产水量的影响更是鲜有研究。

水在维持流域生态系统平衡和生态承载力可持续性方面扮演着极其重要的角色[37]。干旱区水资源空间分布不均衡,降水多集中于上游山区,上游的出山口径流量基本代表着流域的总水资源[38]。疏勒河流域地处青藏高原东北缘祁连山西段,是我国河西走廊内陆干旱区三大内陆河流域之一,也是河西农牧民群众的“生命线”和“天然水塔区”[39]。然而,疏勒河流域产水量相关的研究报道甚少。开展上游产水量研究,对整个流域水资源的动态评估、生态平衡的维系以及气候调节具有十分重要的意义[10,15]。鉴于此,基于InVEST模型产水量模块,本研究对疏勒河上游高寒山区产水量时空变化特征进行了分析,并进一步探讨其对环境要素变化的响应,旨在为流域水资源的有效管理和可持续发展提供科学参考。

1 研究区概况

疏勒河流域地处我国西北干旱区腹地,是河西走廊仅次于黑河的第二大内陆河[40—41],其水源补给主要源于冰川和冻土融水以及山区降水[37,42]。疏勒河上游(96°37′12″—98°59′24″ E,38°13′12″—39°52′12″ N)是指以昌马堡水文站(96°51′ E,39°49′ N)为出水控制点的上部区域,为疏勒河干流的形成区[42],流域面积约为10973.9 km2,海拔2031—5763 m(图1),多年平均径流量为13.81×108m3。地形上主要由疏勒南山、托来南山和疏勒河谷地组成,山地地形陡峭,河谷地形相对低缓[38]。气候类型为典型的大陆性气候,按陆地表层干湿指数区划属于半干旱区—干旱区边缘地带,呈现干寒气候条件下的草原—荒漠景观[38—40]。土壤类型主要为高山寒漠土、高山草甸草原土、栗钙土、淡栗钙土和山地灰钙土[40],土壤质地主要以壤土为主,砂土和粘土占比均较小。疏勒河上游是整个流域的水源涵养区和产流区,为中下游山前平原和绿洲地区赖以生存和发展提供了水源保障。

图1 研究区位置和水文站点分布图Fig.1 Location of study area and hydrological station

2 研究方法与数据处理

2.1 产水量模型

InVEST产水量模块基于水量平衡原理,用各栅格单元的降水量减去实际蒸散后的水量,进而得到该栅格产水量。该模型假设每个栅格单元所有的产水量通过地下径流或地表径流的方式到达流域出水口[14]。模型主要算法如下[7—8,43]：

(1)

式中,Yxj为第j土地利用类型栅格x的产水量(mm),AETxj为第j土地利用类型栅格x的年实际蒸散发量(mm),Px为栅格x的多年平均降水量(mm)。

(2)

式中,Rxj为第j土地利用类型栅格x的Budyko干燥指数,ωx表示植被有效含水量与年均降水量的比值：

(3)

(4)

AWCx=Min(DS,Dr)×PAWCx

(5)

式中,AWCx为栅格x的植被有效含水量(mm)(由有效土壤深度和土壤质地决定)；Z为季节参数(即Zhang系数)用于表征降水的季节性特征,kxj为第j土地利用类型栅格x的植被蒸散系数,ET0x表示栅格x的参考蒸散量(mm);Min为minimum的缩写,DS为土壤深度(mm),Dr为根系深度(mm),PAWCx为栅格单元x的植被可利用水分(mm)。

2.2 参数来源及分析

2.2.1模型参数

产水量模块所需输入参数有气象、土地利用、土壤深度、植被可利用水、生物物理参数表、流域以及Zhang系数。另外,相关分析涉及NDVI。本研究所有数据的空间分辨率均重采样到500 m,且坐标系均使用世界大地坐标系统(WGS1984)。

(1)降水量和参考蒸散发

利用中国气象科学数据中心(http://data.cma.cn/user/toLogin.html)收集周边气象站2001—2019年的降水量和气温数据。用ANUSPLIN模型进行空间插值,得到降水量和气温的空间栅格数据。本文选用Modified—Hargreaves法计算参考蒸散发,与其它方法相比,该方法只需最高气温、最低气温和太阳大气顶层辐射便可计算[45]。公式如下：

(6)

式中,Tmax和Tmin分别为月最高气温和月最低气温(℃),Ra是月太阳大气顶层辐射(MJ m-2d-1),C、E和T是常数,其值参考胡庆芳等[46]推荐青藏高原C、E和T值。

(2)土地利用

土地利用数据源于美国宇航局地球科学数据系统(https://search.earthdata.nasa.gov/search)中的MCD12Q1数据产品。它是是根据一年的Terra和Aqua观测所得的数据经过处理,描述土地覆盖的类型。MCD12Q1采用了五种不同的土地覆盖分类方案,信息提取的主要技术是监督决策树分类。本文采用国际地圈—生物圈计划(IGBP)的分类方法。研究区土地利用类型主要包括草地、永久性湿地、农田、建设用地、农作物—自然植被交错带、永久性冰雪和裸地。

(3)土壤数据和PAWC

基于兰州市寒区旱区科学数据中心下载的第二次土壤调查数据库生成土壤深度和土壤质地栅格图。PAWC是指土壤土层中为植物生长提供的水量所占比例,为0—1的小数。周文佐等[47]提出了一种根据物理和化学性质计算PAWC的方法,具体算法如下;

PAWC=54.509-0.132Sand-0.003(Sand)2-0.055Silt-0.006(Silt)2-0.738Clay+0.007(Clay)2-2.668C+

0.501(C)2

(7)

式中,Sand为土壤砂粒含量(%),Silt为土壤粉粒含量(%),Clay为土壤粘粒含量(%)；C为土壤有机质含量(%)。

(4)流域和生物物理参数表

数字高程模型(DEM)由资源环境科学与数据中心(http://www.resdc.cn/)提供的雷达地形测绘SRTM(90 m)数据。基于数字高程模型数据通过ArcGIS 10.4中的水文分析工具生成的子流域。生物物理表主要用来反映土壤覆盖度和LULC的属性,包括LULC编码、植物蒸散系数(Kc)和根系深度等,生物物理表中每种LULC类型的系数来源于文献[45,48]和《InVEST用户指南》[49]。

(5)Zhang系数

Zhang系数是反映当地降水格局和水文地质特征的气候季节因子,其值在1—30之间变化。根据昌马铺水文站多年平均径流量观测资料,采用水量平衡法对Zhang系数进行校验[7],当Zhang系数为3.33时,InVEST产水模块的模拟效果最好。

(6)NDVI与植被根系深度

NDVI可以准确反映地表植被覆盖状况,是各种植被指数中最常用的表征地表植被状况的指标[36]。本研究中所涉及的2001—2019年NDVI数据源于资源环境科学与数据中心(http://www.resdc.cn/)。该数据集是基于SPOT/VEGETATION NDVI卫星遥感数据,在月数据的基础上采用最大值合成法生成年度NDVI数据集。植被根系深度栅格数据是InVEST模型基于生物物理参数表生成的。

2.2.2时空变化趋势分析

本研究主要通过OriginPro 9.1中的线性拟合方法分析研究目标(气温、降水量和实际蒸散量、NDVI和产水量)的时间变化趋势。利用Matlab—R2015B对研究目标的空间变化趋势进行线性回归分析,该方法广泛应用于时序数据动态检测,线性回归的斜率是量化动态变化趋势的最佳指标。具体算法如下[50]：

(8)

式中,n为研究年限(本研究n=19),i为年份序号,Mi为第i年的研究目标值,Slope>0表示19年来研究目标变化趋势为正,Slope<0则表示研究目标变化趋势为负。

利用线性回归的相关系数对研究目标动态趋势进行显著性检验,相关系数的计算公式如下：

(9)

式中,i为研究期间的年份序号,Mi为第i年的研究目标值,cov为协方差函数,var为方差函数。

2.2.3产水功能重要性分级

合理分配、利用和保护水资源的前提是需对区域产水功能重要性进行分级。本研究基于InVEST模型对研究区产水量进行了定量评估,并依据《国家生态保护红线—生态功能红线划定技术指南(试行)》中涵养水源重要生态功能区的分级方法[51],利用ArcGIS 10.4中的Quantile(分位数)功能对2001—2019年的年均产水量进行分类。按产水量值将其重要性由低到高依次划分为4个重要性级别,即一般重要[0.00—98.38)、中等重要[98.38—227.33)、高度重要[227.33—344.76)和极重要[344.76—488.72)。

2.2.4统计分析

利用MODIS Reprojection Tools与ArcGIS对所获取的栅格数据进行拼接、提取、裁剪等操作,以得到本研究所需的数据集。其次,将2001—2019年昌马堡水文站实测地表年径流量与模拟年总产水量进行线性拟合,并计算调整R2、P以校验模型结果；拟合直线斜率与1∶1直线斜率间的差异通过SPSS 22.0中的一般线性模型进行分析。另外,利用R v4.0.4对环境要素与产水量之间相关性进行Pearson相关矩阵分析。

3 结果与分析

3.1 模型校验

将模型预测结果与实测数据进行比较,通过不断调整参数,得到最适于研究区的预测模型则为模型校验[14]。本研究利用昌马堡水文站地表径流量实测数据对InVEST模型产水量的预测结果进行了校验。结果显示(图2)：模型模拟的多年平均总产水量为13.66×108m3,水文站实测多年平均总径流量为13.81×108m3,其相对误差为1.11%。同时,模型模拟的总产水量与实测地表总径流量之间存在较强的线性关系(R2= 0.986P<0.01)。虽然,模拟总产水量相较实测地表总径流量略低估了0.37×108m3/a,但两者所拟直线斜率与1没有显著差异(P>0.05)。可见,InVEST模型在本研究区适用性很好。

图2 InVEST 模拟总产水量与实测地表总径流量对比 Fig.2 InVEST modelled total water yield vs measured total surface runoff红线代表斜率为0,截距为1的直线

3.2 产水量时空分布及变化特征

2001—2019年疏勒河上游产水量显著增加,增加速率为0.26×108m3/a(图3)。2003年总产水量最低(9.76×108m3),2017年总产水量最高(17.25×108m3)。从空间分布格局来看,山地高海拔区产水量较高,河谷低海拔区产水量较低,尤其是东南部低海拔区与西北部东侧最低(图4)。2001—2019年,所有栅格的产水变化量空间趋势如图4所示。产水量显著增加区域约占整个研究区的32.01%,主要分布在山地高海拔地区；仅有8.39%的区域显著下降且集中分布在研究区的西北部。

图3 2001—2019年产水量时间变化Fig.3 Temporal variation of water yield in 2001—2019

图4 产水量空间分布及变化Fig.4 Spatial distribution and its variation of water yield

3.3 产水功能重要性分级

研究区产水功能重要性分级总体呈现出由河谷低海拔区向山地高海拔区递增的阶梯分布趋势(图5)。其中,一般重要区域所占面积最大,比例为47.26%,主要分布在东南部低海拔区与西北部东侧；中等重要区域所占比例为24.48%,主要分布在西北部；高度重要和极重要区域所占比例分别为14.51%和13.74%,集中分布于山地高海拔区。

图5 产水功能重要性空间分级Fig.5 Importance spatial classification of water yield function

3.4 影响流域产水量的因素分析

疏勒河上游的平均海拔高于3945 m,导致该地区气候寒冷干燥,年均气温为-6.06℃,年降水量325.96 mm。2001—2019年研究区年均气温、年降水量和年实际蒸散发的时间变化趋势如图6所示。年均气温与年降水量分别以0.02 ℃/a (P>0.05)和5.92 mm/a(P<0.05)的速率增加,研究区气候呈现暖湿化趋势。在此影响下,年实际蒸散量和NDVI分别以3.00 mm/a (P>0.05)和0.003/a (P<0.05)的速率递增。

图6 年均气温、年降水量、年实际蒸散量和年NDVI的时间变化Fig.6 Temporal variations of the mean annual air temperature, annual precipitation, annual actual evapotranspiration and annual NDVI

研究区气温、降水量、实际蒸散量和NDVI的空间分布均具有明显的空间异质性。具体而言,气温的空间变化范围为-16.92—5.33℃,自西北向东南逐渐下降且河谷低海拔区高于山地高海拔区,2001—2019年研究区大部分区域(76.08%)的气温呈明显上升趋势,但未达到统计学显著水平(P>0.05)。空间上,降水量在127.42—527.13 mm之间变化,总体呈西北低、东南高且河谷低海拔区低于山地高海拔区的分布格局(图7),与气温的分布情况恰好相反。研究期间,降水量整体显著增加(P<0.05)且东南局部地区(24.41%)极显著增加(P<0.01)。实际蒸散量的空间变化范围为0.00—357.54 mm, 其高值区主要位于东南部低海拔区与西北部东侧；2001—2019年实际蒸散量显著增加区域约为33.39%(P<0.05),主要分布河谷低海拔区,而显著下降区域约为4.01%(P<0.05),散落分布在整个研究区。NDVI空间分布格局与实际蒸散量的分布情况一致,研究期间,NDVI显著增加区域约为56.92%(P<0.05),分布在河谷低海拔局部地区,极显著下降区域约占11.27%(P<0.05),主要分布在山地高海拔区。此外,植被根系深度的空间分布与NDVI一致。

图7 环境要素空间分布及其变化Fig.7 Spatial distributions and the variations of environmental elements

Pearson相关分析表明(图8),产水量与降水量显著正相关,却与气温、实际蒸散量、NDVI和根系深度显著负相关(P<0.001)。不同环境因素之间的相关性表明,实际蒸散量与根系深度、NDVI和气温显著正相关(P<0.001)。气温升高或NDVI、根系深度增加均可能会促进实际蒸散量增加,从而对产水量产生负向影响。

图8 产水量与影响因子的相关矩阵Fig.8 Correlation matrix of water yield and driving factors***, P<0.001

4 讨论

产水量是一项重要的生态功能指标,直接影响着社会经济和生态系统的可持续发展[34]。本研究着重对疏勒河上游产水量的时空变化特征以及影响因素进行了探讨分析。结果表明,2001—2019年疏勒河上游产水量显著上升,这同吕乐婷等[24]对1980—2016年三江源区的研究结果类似,却同王雅等[25]对2000—2014年祁连山黑河中游的研究结果相反。这可能与各研究区的气候条件以及研究的时间尺度不同有关。产水量是区域水循环过程中综合考虑收(降水量)支(实际蒸散量)平衡的产物[10],因此,区域降水量和实际蒸散量以及两者之间平衡关系是影响产水量和产水能力的主导要素,这是许多产水量研究中的一致结论[7,52—53]。潘韬等[7]和吕乐婷等[24]对三江源产水量的研究得出,降水和实际蒸散发是决定生态系统产水量的两个关键环节,其中降水量与产水量显著正相关,实际蒸散量与其显著负相关。刘军会等[52]和刘娅等[53]分别对青藏高原和京津冀地区的产水量进行了研究,结果表明降水量和实际蒸散量综合影响区域产水变化量。然而,气温、NDVI和植被根系深度会通过影响实际蒸散量来间接影响产水量。其中,NDVI和植被根系深度主要通过下垫面植被变化来影响实际蒸散量,进而影响产水量。Lu等[35]对青藏高原拉萨河流域产水量的研究发现,降水量和NDVI变化是影响区域产水量的关键要素,特别是通过植被变化来推动产水量变化。19年来,疏勒河上游实际蒸散增加趋势(3.00 mm/a)低于降水增加趋势(5.92 mm/a),降水量对产水量的补给远超过实际蒸散对产水量的消耗,因此,疏勒河上游产水量整体呈显著增加趋势。

疏勒河上游发育大面积多年冻土,其中山地高海拔区的土层存在冻结面,这使得降水无法下渗而在冻结面汇流形成径流。与此同时,山地高海拔区的降水量高,气温和实际蒸散量低且土地利用类型主要为冰川和裸地。冰川属于水域,最易融化形成径流；裸地的植被覆盖度低且根系浅,截留能力弱,更多降水也直接形成径流[24,54—55]。相比之下,河谷低海拔区由于冻土分布面积和深度受限,降水或冰水可以渗入土壤深层,径流形成的几率大大降低。另外,河谷低海拔区的主要土地利用类型为草地,其植被覆盖度和实际蒸散量都很高,并且植被根系深而密集,可以有效的截留降水,使得产水量偏低。2001—2019年山地高海拔区产水量显著增加,主要是因为降水量补给显著增加而蒸散量无显著变化,致使生态产水以增加为主。与其相反,西北地区的裸地被草地和农田大面积取代,相比裸地,草地与农田的根系较深且植被生长消耗的水量也比较大,因此产水量显著下降。

疏勒河上游产水重要性分级与产水量空间分布格局基本一致,这与孙艳伟等[14]和王玉纯等[51]研究结论一致。具体而言,产水极重要和高度重要区降水量较多、蒸散量较小,并且多分布冰川和裸地,是疏勒河流域水资源集成区,而产水一般重要和中等重要地区为相对干旱的堆积平原和盆地,为水资源消耗区。为使研究区实现健康可持续发展,极重要和高度重要区需作为重点和优先保护的水资源地而纳入政策和管理决策之中[56]。

本研究对疏勒河上游产水量的模拟结果较为理想,基本符合区域的实际情况。但因模型自身设定、参数算法以及数据精度等原因,使研究仍存在一些不确定性。首先,InVEST模型简化了汇流过程,没有区分地表径流、土壤中径流和基流,并假定所有栅格的产水量都到达同一出水口[14],并且产水量常与社会经济和人类活动密切相关,而模型的输入数据只有自然数据。特别是,疏勒河上游是一个典型的高寒山区,InVEST模型未考虑冰川和冻土融水,这会增加水文模拟的不确定性。因此,InVEST模型在这些方面都有待进一步改进。其次,降水量和参考蒸散量是精确模拟产水量的重要参数,而西北地区水文站点较少且气候变化具有复杂性,利用ANUSPLIN模型空间插值生成的降水数据以及通过Modified—Hargreaves法计算的参考蒸散发的准确性都会受到气象数据精度的影响[57],从而影响模型的模拟精度。尽管存在这些不确定性,但研究结果仍可以有效反映产水量的时间变化和空间分布,揭示产水量与环境要素之间的关系,为疏勒河流域水资源的有效管理、合理利用和保护提供科学参考。

5 结论

基于InVEST模型,结合环境要素对疏勒河上游产水量的时空分布格局及影响因素进行分析,得出：

(1) 模型中的Zhang系数为3.3时,研究区产水量模拟效果最佳；研究区年均总产水量为13.66×108m3,2001—2019年研究区产水量以0.26×108m3/a速率显著增加。从空间分布来看,山地高海拔区产水量较高,河谷低海拔区产水量较低,尤其东南部低海拔区和西北部东侧最低。研究期间,产水量显著增加区约占总面积的32.01%,主要分布在山地高海拔区；仅有8.39%的研究区显著下降且集中分布于西北部。

(2) 产水量与降水量显著正相关,而与气温、实际蒸散量、NDVI和根系深度显著负相关；实际蒸散量与气温、NDVI和植被根系深度显著正相关。

(3) 在对研究区产水量进行定量评估的基础上,进行了产水功能重要性空间分级。得出：一般重要区域所占面积最大,约占整个研究区的47.26%,主要位于东南部低海拔区与西北部东侧；高度重要和极重要区域共占总面积的28.25%,集中分布于山地高海拔区。