丹江口库区产水量时空动态与情景模拟*

张 建 雷 刚 漆良华,

(1.国际竹藤中心 国家林业和草原局/北京市共建竹藤科学与技术重点实验室 北京 100102； 2.国际竹藤中心安徽太平试验中心 太平 245700)

产水量是流域重要的生态系统服务功能，可表征一定时空范围内生态系统水分保持的过程与能力，既影响区域水资源的整体水平，又反映区域自然环境和人为活动的关系(陈姗姗, 2016； Langetal., 2018)。因此，深入了解区域产水量时空动态变化特征及其主要驱动因素具有重要意义。

随着3S技术与水文生态模型的发展，越来越多的学者通过水文模型分析和评估区域生态系统产水量(Lehetal., 2013)。其中InVEST模型可以较好地把握总体格局，将量化的生态系统服务功能和情景以地图的形式表达，可设置不同情景来模拟未来土地利用条件下生态系统服务功能和多种服务之间的权衡关系(刘智方等, 2016； 唐尧等, 2015)。该模型由美国斯坦福大学、大自然保护协会、世界自然资金会共同开发，输出结果的空间化表达使生态系统服务的重要区域易于识别(Rimaletal., 2019)。InVEST模型产水量模块被Leh等(2013)、Neson等(2009)应用于西非加纳和科特迪瓦、美国俄勒冈州的威拉米特河流域； 顾晋饴等( 2018)利用InVEST模型发现太湖流域水源涵养空间分布差异性和不均衡性不断增大，降水量对水源涵养能力的影响大于土地利用/覆被类型的影响； 程一凡(2019)利用模型评估了三江源国家公园近十年水源涵养量变化，得到不同时期不同区域水源涵养量变化特性，但对驱动因素的探讨依然空缺； 邱问心等(2018)以浙江临安区水涛庄水库集水区为例，基于水量平衡法、综合蓄水法和InVEST模型法评价研究区水源涵养功能，结果表明模型结果精度可达83.36%，InVEST模型测算水源涵养功能是可行的且具有较高推广价值。

丹江口市为南水北调中线工程水源区，近年来城乡发展、人口聚集和土地利用/覆被类型变化等影响了库区水源涵养服务能力，因此为学者们所关注。丁霞等( 2019)对库区马尾松(Pinusmassoniana)林水源涵养价值进行了研究，徐慧( 2017)利用遥感技术探究了南水北调中线工程的初始运行对水源地带来的影响，这些研究多局限于某一种生态系统类型的单一生态服务价值(李亦秋等, 2011； 柳晶辉等, 2008)，缺乏时空动态变化和驱动因素的研究。因此，本研究以南水北调中线工程水源区丹江口市为研究区域，采用InVEST模型研究库区2003、2013和2018年产水量时空动态，并利用情景模拟方法探讨年降水量和土地利用/覆被类型变化对产水量的影响，以期为库区水源涵养能力提升与调控提供科学依据。

1 研究区概况

丹江口市地处湖北省西北部，十堰市东部，汉江中上游(110°47′53″—110°34′47″E，32°14′10″—32°58′10″N)，是南水北调中线工程重要水源地，东西最大横距73 km，南北最大纵距81 km。南水北调中线工程于2003年开工，同年丹江口市开展移民搬迁工程； 2013年南水北调中线水源地工程完工，开展调水试验，同年移民搬迁工程圆满完成； 截止2018年底丹江口库区累计供水2.22亿m3。丹江口市属亚热带半湿润季风气候区，气候较温和，四季分明，光照时间长，雨热同季且无霜期较长，年均气温15.9 ℃，年均降水量约800 mm，年蒸发量1 979.1 mm。土壤以黄棕壤和黄壤为主，成土母质由石灰岩、片麻岩等发育而成，质地疏松。马尾松、栓皮栎(Quercusvariabilis)、柏木(Cupressusfunebris)、柑橘(Citrusreticulata)等在库区均有分布。

2 研究方法

2.1 数据来源及处理

InVEST模型所需数据包括研究区土地利用/覆被类型图、年均降水量、年均潜在蒸散量、土层厚度、植物可利用水含量、集水区范围、季节性因子Z值等。

1) 土地利用/覆被类型图 通过遥感影像图进行解译分类得到。原始遥感影像数据源于美国地质调查局(http:∥www.usgs.gov/)，选择丹江口市2003年Landsat 7 ETM+、2013年和2018年Landsat 8 OLI_TIRS 影像，分辨率为30 m，云量均3%以下，影像质量良好，在ENVI软件中通过几何校正、融合镶嵌、裁剪处理，采用监督分类方法进行解译，将研究区土地利用类型分为林地、耕地、园地、建设用地、水域和未利用地6类(图1)。

图1 2003—2018年丹江口市土地利用/覆被类型Fig. 1 Land use/cover types of the Danjiangkou City in 2003—2018

2) 年均降水量 年降水数据下载于中国气象数据网(http:∥data.cma.cn/)，为规避单一年份数据代表性低，选取3个时期(2002—2004、2012—2014、2017—2019 年)气象数据的均值，整理研究区周边气象站点及监测设备降水量数据，通过克里格插值得到年均降水量栅格图。

3) 年均潜在蒸散量 利用中国地面气象观测报告，收集研究区周边4个气象站数据。利用Hargreaves方程(Droogersetal., 2002)计算年均潜在蒸散量ETo，通过克里格插值得到年均潜在蒸散量的栅格图。

ETo=0.001 3×0.408×RA×(T+17)×

(TD-0.0123P)0.76。

式中： RA为太阳大气顶层辐射量；T为年最高温均值和最低温均值的平均值；TD为年最高气温均值和年最低气温均值之差。

4) 土层厚度 根据森林资源二类调查数据结果中的属性数据，采用克里格插值，得到模型所需图层厚度的栅格图。

5) 植物可利用水含量 表征土壤为植被生长所储蓄的总水量，参考周文佐等(2003)提出的计算公式，将计算结果进行克里格插值得到模型所需栅格图。

PAWC=(54.509-0.132SAN-0.003SAN2-0.055SIL-0.006SIL2-0.738CLA+0.007CLA2-2.688C+

0.501C2)×100%。

式中： PAWC为植物可利用水含量；SAN为土壤的沙粒含量(%)； SIL为土壤粉粒含量(%)； CLA为土壤黏粒含量(%)；C为土壤有机质含量(%)。

6) 集水区范围 指被分水岭所包围的集水流域，内部包含若干个互不嵌套的子流域。从地理空间数据云网站(http:∥www.gscloud.cn)下载研究区30 m分辨率高程数据，经镶嵌、裁剪、填洼后生成DEM数据； 基于DEM使用ArcGIS水文分析工具和矢量数据生成流域和子流域，每一级流域赋予唯一数值。

7) 季节性因子Z值 是降水特征的常数，通常用Zhang系数表征(Zhangetal., 2001)。对于降水总量相等的区域，降水次数越多该值越大。本研究经多次调节Z值，确定季节性因子取值为11时，模型评估结果与实测值最为接近。

2.2 年产水量模型原理

InVEST模型产水量模块是基于Budyko水热耦合平衡假设和年均降水的水量平衡估算方法建立的(潘韬等, 2013)。模型不区分地表水、地下水和基流，某栅格单元的年降水量减去年蒸发量即年产水量，参考InVEST模型手册(Redheadetal., 2019)，其算法为：

式中:Yxj为栅格x中土地利用/覆被类型j的年产水量； AETxj为栅格x中土地覆被类型j的年蒸散量；Px为栅格x的年降水量。

式中：Rxj为为第j类土地利用/覆被类型在栅格单元x的干燥度指数，为年蒸散量与年降水量的比值，无量纲；wx为植物年需水量和年降水量的比值，无量纲。

式中：Kxj为栅格单元x中土地利用/覆被类型j的植被蒸散量的系数； ETox为栅格单元x中的年均潜在蒸散量。

式中： AWCx为土壤有效含水量，由土壤质地和有效土壤层厚度决定；Z值为季节性因子，代表季节性降水分布和降水深度，不同研究区Z值存在差异。

2.3 情景模拟法

由于InVEST模型产水量模块是基于水量平衡法原理开发的，年降水量和年蒸散量是影响结果的主要因素，而年蒸散量主要受植被和土地利用/覆被类型影响(潘韬等, 2013； Terradoetal.， 2014)。为了进一步探究年降水量和土地利用/覆被类型变化对产水量的影响，本研究设计2种情景，情景1： 模型均输入2003年降水数据，土地利用/覆被数据为2013年和2018年实际数据，与2003年结果相比，即年降水量数据保持不变而土地利用/覆被类型变化，可研究2013年和2018年土地利用/覆被对产水量的影响； 情景2： 模型均输入2003年土地利用/覆被数据， 2013年和2018年降水量数据为相应的年降水量，与2003年结果相比，可探究年降水量变化对产水量的影响。引入实际产水量总变化量WT(赵亚茹等， 2019； 杨洁等， 2020)，其正负值表示年降水量和土地利用/覆被类型变化使产水量增加或减少。

ΔWC=WO-WC;

ΔWL=WO-WL;

ΔWT=ΔWL-ΔWC。

式中：WO为实际产水量；WC为年降水量不变而土地利用/覆被类型变化的产水量；WL为土地利用/覆被类型不变而年降水量变化的产水量；WT为各模拟的总变化量，即在土地利用/覆被类型变化和年降水量变化双重影响下产水量的变化量； ΔWC为实际情况与年降水量不变情景差值； ΔWL为实际情况与土地利用/覆被类型不变情景的差值。

2.4 模型检验与校正

土地利用/覆被图是模型运行的数据基础，其精度直接影响结果的输出，采用Kappa分析方法分别对3期土地利用/覆被图结果进行精度验证。

在模型其他参数确定的情况下，季节性因子Z值是影响模型准确性的重要参数。将产水量模块运行结果与实际监测数据进行拟合，通过调整季节性因子Z值大小，确定最佳值。选择《十堰市水资源公报(2013年)》中丹江口市2013年产水量为标准数据，与模型计算的产水量进行检验校正。

3 结果与分析

3.1 模型校验

在对土地利用/覆被类型的校验中，利用谷歌地图获取3个时期的高清历史影像图，结合野外调查数据和历史影像图建立样本区，对监督分类解译的数据修正，经Kappa系数检验，3期土地利用/覆被图精度达到0.84、0.86、0.85，解译结果符合精度要求。

将所需数据按要求输入模型后，在其他数据已确定情况下，其运行结果受季节性因子Z值的影响。将年产水量模型运行结果与监测数据进行拟合，结果发现随着季节性因子Z值增加，产水量呈减少趋势(图2)，确定季节性因子Z值为11时，模型结果与监测值最为接近，误差为0.72%。

图2 模型校验结果Fig. 2 Model verification result

3.2 产水量时空动态

整体上看(图3)，丹江口市近15年产水量呈先减弱后增强的趋势， 2003年单位面积产水量最深可达646.96 mm，远高于2013年(451.25 mm)和2018年(527.53 mm)。2003年丹江口市产水深度310.09 mm，产水量为13.6亿m3，产流系数为0.324； 2013年产水能力最弱，产水深度146.67 mm，产水量为5.25亿m3，产流系数0.281； 2018年产水量相较2013年增加了27.88%，产水量为7.28亿m3，产水深度为166.06 mm，产流系数0.296。

图3 2003、2013和2018年丹江口市产水深度时空分布Fig. 3 Spatial and temporal distribution of water yield depth in Danjiangkou in 2003, 2013， and 2018

同一年内，丹江口市产水能力在空间上分布不均匀，北部丘陵山区产水深度高出南部39.84%～57.79%； 时间序列下， 2003—2013年，全市面积的85.10%为减弱区，13.50%的区域表现为增加趋势(主要分布在林地建设区)； 2013—2018年，丹江口市89.60%的区域产水深度有所增加，而10.20%的区域表现出减少的趋势(图4)。

图4 2003—2018 年丹江口市产水深度变化Fig. 4 Variation of water yield depth in Danjiangkou from 2003 to 2018

下垫面是影响产水量变化的重要因素，利用分区统计功能得到不同土地利用/覆被类型下的产水深度变化情况(图5)。丹江口市4种土地利用/覆被类型中，产水深最高的是建设用地(530.35 mm)，其次分别为耕地(292.85 mm)、园地(284.07 mm)和林地(273.76 mm)； 相较于其他土地利用/覆被类型，建设用地无植被截留降水，蒸散量小，因此产水深度最好。从变化趋势上可看出， 2003—2018年各土地利用/覆被类型的产水深度呈先减小后缓升的趋势； 总体而言，耕地保持增长而其他土地利用/覆被类型则降低。

3.3 情景模拟分析

总体上看(图6)，在年降水量不变而土地利用/覆被类型变化情景下，情景1中2013和2018年产水量分别为10.5亿m3和8.6亿m3，表明土地利用/覆被类型变化是导致产水量减少的重要因素； 在土地利用/覆被类型不变而年降水量变化情景下，情景2中2013年和2018年产水量变幅较小，分别为5.68亿m3和9.7亿m3，说明年降水量是影响产水量变化的重要因素。

图5 不同土地利用/覆被类型产水情况Fig. 5 Water yield of different land use types

表1表明了年降水量变化和土地利用/覆被类型变化对产水量的影响程度。2013年，两种情景模拟结果显示，土地利用/覆被类型变化影响下实际产水量减少了0.43亿m3，年降水量变化使实际产水量减少了5.25亿m3，二者的贡献率分别为7.57%和92.43%，说明2003—2013年降水量与土地利用/覆被类型变化均使得产水量减少，且年降水量是影响产水量的主要因素。而2018年，土地利用/覆被类型变化导致实际产水量减少了2.42亿m3，贡献率为64.71%，年降水量变化使实际产水量减少了1.32亿m3，说明相较于年降水量，土地利用/覆被类型变化对2013—2018年产水量的影响更大。

图6 模拟不同情景下产水量变化Fig. 6 Simulate the change of water yield under different scenarios

表1 年降水量与土地利用/覆被类型变化对产水量的影响Tab.1 Annual precipitation and land use/cover type change influence water yield

在仅考虑年降水量不变而土地利用/覆被类型变化情况下， 2013—2018年水域和林地面积持续增加，特别是林地面积增加了10.03%，耕地和园地面积不断减少，建设用地先减少后增加，未利用地变化不明显(表2)。由于建设用地基本为不透水面，耕地和园地的植物密度、根系深度及人为劳作带来的地表变化等都有助于产水量的增加，而林地通过冠层截留降水、枯落物层吸收降水、土壤层蓄渗降水以及其较深的根系可有效拦截降水，且有强大蒸腾作用，从而使得产水量较少。

表2 丹江口市2003-2018年不同土地利用/覆被类型面积变化Tab.2 Area changes of different land use/cover types in Danjiangkou from 2003 to 2018

4 讨论

Lü等(2013)研究了中国东南部南北样带1981—2000年产水量和气候及土地利用/覆被类型变化之间的动态响应，其中产水量的时空变化受年降水量影响最显著，与本研究结果相似； 王耕等( 2018)利用年产水量数据进一步对大凌河上游汇水区水源涵养功能定量评价，发现产水量与水源涵养均有明显的空间异质性； Wang等(2020)评估了时空变化下丹江口库区景观格局和土地利用/覆被对水资源的影响，对未来土地利用管理和政策的实施具有重要意义。当前产水量研究方法中均存在一定的优势和局限性，而InVEST模型考虑了蒸散发、降水量和多年径流系数等数据，研究结果更为准确。

本研究利用InVEST模型对丹江口市2003—2018年产水量开展了研究，结果表明近15年产水量时空变化较大，产水量分布不均，且各土地利用/覆被类型的产水深度存在差异。研究结果与太湖流域(顾晋饴等， 2018； 陈骏宇等， 2016)、三江源流域(程一凡， 2019； 吕乐婷等， 2020)大凌河流域(王耕等， 2018)、石羊河流域(赵亚茹等， 2019)、大连市(吕乐婷等， 2019)等流域或地区的结果相似，其中杨洁等(2020)研究了黄河流域不同土地利用/覆被类型的产水深度，发现永久性冰川雪地、裸岩石质地、城镇建设用地的产水深最高，与本研究结果中建设用地产水量最高相同。

产水量大小受该区域年降水量、蒸散量以及二者之间平衡关系的影响，而实际蒸散量除了受气象因素影响外，还直接受制于土地利用/覆被类型的影响。Lang等(2017)量化了气候变化和土地利用/覆被对产水量的影响，发现年降水量对年产水量变化的贡献率高达97.44%，与本研究情景1中的模拟结果相似。土地利用/覆被类型变化通过改变下垫面结构和类型直接或间接影响产水量，其次是土壤孔隙度、土壤质地、结构均对产水量带来间接影响(Jatinetal., 2018； Pamukcuetal., 2016)。因此产水量变化是多因子耦合驱动的过程(Lüetal., 2017； Suetal., 2013)，集合了包括气候变化与土地利用/覆被类型变化在内的自然、社会、经济等多方面因素的共同作用。

本研究忽略了上游对该区域的水源补给，可能会导致InVEST模型拟合结果偏大； 季节性因子Z值最终设定为11，但由于小流域气候等因素也会存在一定误差； 模型结构的本土化不足以及未能获得长期野外观测数据，会对研究结果的准确性产生影响。在今后研究中，应加强对野外数据的监测收集、调整模型参数使其本土化，确保评估结果的可靠性。

5 结论

本研究基于InVEST模型，以丹江口市域范围为研究区，定量研究了2003、2013和2018年产水量的动态变化特征，并运用情景模拟法探讨了年降水量和土地利用/覆被类型变化对产水量的影响。结果表明： 1) 丹江口市2003—2018年产水量表现出先减弱后增强的趋势；2) 单位面积产水量在空间上分布不均，平均而言北部丘陵山区高出南部39.84%～57.79%；3) 各土地利用/覆被类型产水深度均有差异，其中产水深度最高的是建设用地(530.35 mm)；4) 情景模拟表明年降水量和土地利用/覆被类型是影响产水量的主要因素， 2013年产水量中年降水量贡献率为92.43%，土地利用/覆被类型变化对2018年产水量贡献率为64.71%。