近30年千岛湖流域产水量时空变化及其影响因子分析

朱志洪,周本智,王懿祥,祁 军,李爱博,黄润霞

(1.中国林业科学研究院亚热带林业研究所,浙江 杭州 311400;2.浙江农林大学环境与资源学院,浙江 杭州 311300;3.上海春沁生态园林建设股份有限公司,上海 201210)

生态系统服务是指生态系统直接或间接地为人类生存提供各种服务资源[1]。流域是地表上的特殊生态系统,是水资源在空间上的集合。产水量的实质是降水量与实际蒸散量的差值,表征着流域的水供给能力[2],是重要的生态系统服务功能之一,对调节和改善流域水分循环十分重要,同时也与流域的其他生态系统服务功能如土壤保持、农作物生产及生物多样性保持等息息相关[3]。因此,探究流域产水服务时空变化及其影响机制对于流域水资源合理开发和实现生态治理具有重要的意义。

InVEST模型因其强大的空间表达能力、数据易获取性、结果可信度高和便于情景模拟等优点,成为目前国内外应用较为广泛的生态系统服务评估模型[4-5]。国内外学者通过模型的产水模块对不同地区的产水情况进行了量化和分析,如Redhead等[6]采用不同数据源对英国不同流域产水量进行量化,与真实资料对比分析,验证了InVEST模型结果的可靠性;Yohannes等[7]量化了贝蕾萨流域的产水量和输沙量两个指标的时空变化,发现两者受景观组成的强烈影响;Marquès等[8]对西班牙的Francoli流域产水量进行量化并预测了未来气候变化对该流域产水量的影响。国内学者也在不同空间尺度进行了产水服务时空变化特征分析,并探究气候和土地利用变化对产水的影响。如杨洁等[9]通过情景分析法量化了气候和土地利用变化对黄河流域产水服务的影响,指出1995—2019年降水变化和土地利用变化对流域产水深度的贡献率分别为84%和16%;戴尔阜等[10]采用地理探测器探究横断山区产水服务的空间异质性,发现在不同地貌和气候分区,产水服务的影响因素具有显著差异;王晓峰等[11]采用相关性分析和地理加权回归对秦岭生态屏障区产水服务驱动要素进行分析,发现降水的主导范围最大。目前对于流域源供给能力研究大部分集中在我国大中型流域以及干旱、半干旱地区,主要通过情景分析法、相关分析、地理加权回归等方法探究气候和土地利用变化对产水的响应[2],但对于东南沿海地区中小型流域的研究较少。

千岛湖流域位于浙江省西北部,是国家一级水资源保护区,也是整个长江三角洲的重要战略饮用水源地[12]。近年来,随着“两山理论”“山水林田湖草沙”一体化保护和系统治理的实践,生态保护和修复已成为千岛湖流域规划发展的重要组成部分,流域的生态系统服务功能逐渐受到研究人员的关注。此前研究主要针对千岛湖流域的水质情况、生物多样性和生态系统服务价值评估等[13-16],对单一水源供给能力的研究较为缺乏。流域的水源供给能力对其生态系统服务功能和社会经济发展都有着直接或间接的影响,不同学者由于研究目的不同对流域范围界定各不相同,相晨等[17]对千岛湖流域主体即淳安县进行生态系统服务价值评估,汤旭光等[18]对千岛湖流域多年植被覆盖度进行观测,流域范围包括淳安县以及安徽省的新安江流域。深入了解流域产水量时空变化特征及其空间格局的影响因子,对于流域的水资源利用和水生态保护有着重要的作用,同时为退耕还林还草、荒山荒地造林等相关环境保护措施提供参考。因此,本研究基于时空尺度探究千岛湖流域近30年产水量时空动态变化特征,通过情景分析法,设计4种不同情景定量讨论降水和土地利用对流域产水服务的影响,并分析单一土地利用类型变化情景对流域产水服务的影响,以期为千岛湖流域实现水资源可持续管理和生态系统可持续发展提供参考。

1 材料与方法

1.1 研究区概况

本研究界定千岛湖流域为整个淳安县县域以及建德市的新安江街道(118°33′～119°34′E,29°18′～30°03′N),流域面积4 526 km2,地势四周高、中间低,属浙西山地丘陵区,上游的新安江及其支流是流域主要汇入河流。流域位于亚热带季风气候北缘,雨量充沛,四季分明,年平均气温17.2 ℃,年平均降水量约为1 515.1 mm,平均日照1 850.3 h,是我国典型的气候湿润区。植被覆盖度较高,以亚热带常绿阔叶林为主,动植物资源丰富,土壤以红壤、黄壤和岩性土为主。流域内的新安江水库是20世纪50年代为建新安江发电站而形成的大型水库,湖区面积573 km2,相应库容178.40亿m3,是钱塘江中下游的主要水源。

底图审图号:GS(2016)2556号。下同。

1.2 研究方法

1.2.1 产水模块原理

InVEST模型产水模块是基于水量平衡原理,通过降水、实际蒸散等参数计算区域产水量,公式如下:

(1)

式中:Yx为栅格单元x的产水量;Aetx为栅格单元x的实际蒸散量;Px为栅格单元x的降水量。两者比值计算公式如下:

(2)

式中:Rx为栅格单元x的无量纲干燥指数;wx为物理参数。

Rx=Kc×Eet0x/Px;

(3)

wx=Z×Aawcx/Px。

(4)

式中:Kc为表征植被蒸散的系数;Eet0x为栅格单元x的参考蒸散量,采用FAO修正的Penman-Monteith方程[19]进行计算;Aawcx为植被含水量;Z为表征季节降水分布的参数,与降水次数正相关。

Aawcx=min(Sdepth,Rdepth)×Ppawc。

(5)

式中:Sdepth为土壤深度;Rdepth为根系深度;Ppawc为植被有效含水量[20]。

(6)

式中:Ssan、Ssil、Scla分别为土壤砂砾、粉粒、黏粒含量(质量分数),%;CSOM为土壤有机质质量分数,%。

1.2.2 情景分析法

1) 气候和土地利用变化情景。InVEST模型产水模块实质是将研究区的降水量与实际蒸散量的差值作为研究区的产水量,其中实际蒸散量与地表的土地类型相关。为更好地研究气候变化和土地利用变化对产水量的影响,本研究设计4种情景(情景1～4),深入研究不同时间段气候和土地利用变化对产水量的影响,并通过贡献度进行量化。将1995—2019年划分为1995—2005和2005—2019两个时间段,对应气候变化中的情景1和2。以情景1为例,在模型输入1995年的土地利用数据和2005年的降水量和潜在蒸散量数据,与1995年真实情景相比,土地利用数据不变,气候要素发生变化,探究1995—2005年气候变化对产水服务的影响。同样以这两个时间段对应土地利用变化中的情景3和4,并以情景3为例,在模型输入1995年降水量和潜在蒸散量数据,以及2005年土地利用变化数据,与1995年真实情景相比,气候要素不变,土地利用数据发生变化,探究1995—2005年土地利用变化对产水量的影响,并采用贡献度进行量化。

R1=Δ1/(Δ1+Δp);

(7)

Rp=Δp/(Δ1+Δp)。

(8)

式中:Rl和Rp分别为土地利用变化和降水量变化对产水量的贡献度,Δl和Δp分别为土地利用变化和降水变化导致的产水量的变化量。

2) 极端土地利用情景。极端土地利用情景以2019年的土地利用数据和气候数据为基准,由于千岛湖流域主要土地利用类型为林地、水体、农田和建设用地,考虑到流域的土地利用现状以及未来的发展规划,设置4种土地利用变化情景(情景 5～8),依次为农田全部转化为建设用地、农田全部转化为林地、建设用地全部转化为林地和建设用地全部转化为农田,其他用地类型保持不变,以分析单一土地利用类型变化对流域产水服务的影响。

1.2.3 地理探测器

地理探测器是探测空间分异特征及其主要驱动力的一种方法,其基本思想是:如果自变量与因变量的空间分布特征具有相似性,那么自变量对因变量就有重要的影响[21]。它包含风险探测、因子探测、生态探测和交互探测4个模块。本研究通过因子探测和交互探测两个模块分析千岛湖流域产水服务空间格局的主要影响因子,计算公式如下:

(9)

式中:q表示影响因子对产水量的解释度,范围为[0,1],值越接近1说明该类因子对产水量空间分布影响越高;h为分区数目,L为影响因子的样本数量;N和Nh分别是全区和层内单元数;σ和σh为全区和层h的方差。交互作用探测可探测两个影响因子之间的交互作用关系,评估其交互作用是增强还是减弱对产水量的解释力。

1.3 数据来源及处理

选取研究区周边气象站的降水、气温等气候资料,采用克里金插值法进行插值得到空间分辨率为30 m的降水和蒸散数据,数据来源于中国气象数据网(http://data.cma.cn)。土地利用数据是基于Landsat TM/OLI遥感影像,结合NDVI数据,MODIS土地利用数据通过随机森林算法分类获得,包括水体、建设用地、农田、未利用地、草地和林地6类。土壤砂砾、粉粒、黏粒数据来自中国科学院资源环境科学数据中心(http://www.resdc.cn),空间分辨率为1 km;土壤有机质数据来源于国家地球系统科学数据中心(http://www.geodata.cn),空间分辨率为250 m;土壤的根系最大埋藏深度数据来源于国家青藏高原科学数据中心(http://data.tpdc.ac.cn)的基于世界土壤数据库(HWSD)的中国土壤数据集,空间分辨率为1 km。数字高程模型(DEM)数据来源于地理空间数据云(http://www.gscloud.cn),空间分辨率为30 m。道路和行政区划数据来源于全国地理信息资源目录服务系统的1∶25万全国基础地理数据库。季节常数采用Donohue等[22]方法通过0.2倍的降水次数进行计算,并在此基础上多次模拟运行,当Z取29.6时,模拟产水量与流域水资源量最接近。植被蒸散系数(Kc),最大根系深度数据通过参考文献[23-24]、联合国粮农组织植被蒸散系数参考指南和InVEST模型推荐的参数获得。

上述所有模型输入数据经重采样后,空间分辨率为30 m,坐标与投影保持一致。由于地理探测器的输入自变量必须为类型量,因此本研究采用ArcGIS10.2通过自然断点分级法对降雨等气候因子和地形因子进行划分。

2 结果与分析

2.1 千岛湖流域土地利用变化及对产水量的影响

千岛湖流域土地利用类型以林地、水体、农田和建筑用地为主,草地和未利用地的面积占比较小(表1)。由表1分析可知,1995—2019年林地面积占比由83.65%上升到85.60%,面积增加88.70 km2,主要由农田和水体转化而来,建筑用地面积占比由1995年的0.51%上升到2019年的1.12%,增幅为121%,主要由农田和林地转化而来,可以看出建筑用地有明显的扩张趋势;农田面积占比由1995年的5.26%下降到2019年的3.13%,且主要向林地和建筑用地转移,说明1995—2019年千岛湖流域在快速城镇化的同时,退耕还林等生态修复政策的施行也获得了较好的成效。水体面积占比分别由1995年的10.56%下降到2019年的9.96%,未利用地和草地面积有明显增加,但整体占比仍较低。

表1 1995—2019年千岛湖流域土地利用变化转移矩阵

不同的土地利用类型由于蒸散能力、土壤导水能力、植被冠层截流能力等的差异,其产水能力也各不相同。从平均产水量来看,建筑用地、草地、耕地和未利用地的产水量较为接近,占所有地类产水量的16.14%～21.32%,其中建筑的产水量最高,草地、耕地、未利用地次之,林地和水体的平均产水量最低,占产水量的5.95%～16.37%。建筑用地地表多为不透水层,降水直接转化为地表径流,因此单位面积建筑用地产水量最高;林地由于其内部较高的植被蒸散和冠层截留降水,植被根系吸收降水等原因,使得降水更少地转化为产水量,故产水能力低于建筑;草地等低植被覆盖地类,水体的蒸散能力最高且降水量相对较低,因此产水量最低。从产水总量来看,林地的产水总量最高,水体次之。2019年千岛湖流域林地占整个流域面积的85.84%,占流域产水总量的87.63%,是千岛湖流域产水总量的主要贡献地类。水体产水总量次之,两者产水总量之和占流域产水总量的94.84%(图2)。

图2 不同土地利用类型产水情况Fig.2 The water yield of different land use types

2.2 千岛湖流域产水量时空变化特征

基于InVEST模型的产水模块得到1995—2019年平均产水总量为55.48×108m3,与相关研究成果[25]和《安徽省水资源公报》发布的千岛湖流域多年平均水资源总量进行对比,当Z系数为29.6时,两者最为接近。

1995—2019年千岛湖流域产水总量年际变化十分明显(表2)。1995—2005年,产水总量呈现出明显的下降趋势,整体下降52.95%;2005—2015年,产水总量逐步上升,从2005年的29.54×108m3升到2015年的78.67×108m3,升幅166.32%;2015—2019年,产水总量再次降低,下降28.22%。从总体趋势来看,1995—2019年间,流域产水总量变化趋势与年降水量变化趋势相同,采用平均产水量表征单位面积产水量,其变化趋势与降水量一致,因此降水量是产水量年际变化的主要驱动因子。

表2 千岛湖流域各年份产水情况

分析可知,千岛湖流域产水量空间分布格局变化不大,整体表现出较一致的规律性(图3)。

图3 千岛湖流域产水量空间分布格局Fig.3 The spatial distribution patterns of the water yield in the Thousand-Island Lake basin

高值区连片分布在中西部千岛湖两岸人口密集区,多年平均产水量在1 400～1 600 mm之间,低值区分布在千岛湖东部湖区水面之上,多年平均产水量在800～900 mm之间。这种分布格局主要与流域年平均降水量和下垫面蒸散能力相关,降水量高,蒸散能力弱的地区产水量高;反之,降水量较低,蒸散能力强的地区产水量较低。

2.3 不同情景下产水量模拟分析

通过模拟气候与土地利用变化可知,在气候变化情景下,与1995年基准年相比,2005年产水量减少739.17 mm(情景1),2019年较2005年增加597.53 mm(情景2),整体变化较为明显;在土地利用变化情景下,与基准年相比,各年份模拟产水量变化较小,情景3产水量增加了3.64 mm,情景4则减少2.66 mm。采用贡献度进行量化可知,1995—2005年气候变化贡献度和土地利用变化贡献度分别为100.50%和-0.50%,2005—2019年,贡献度分别为100.45%和-0.45%,气候变化对产水量的贡献度均远超于土地利用变化,流域的水供给能力主要受气候条件影响(表3)。

表3 不同情景下产水量的变化情况和贡献度

情景3和情景4的不同变化情况则说明了土地利用变化对流域产水量影响的不确定性。这可能与千岛湖流域土地利用变化多表现在耕地的减少和建筑用地的增加,林地、水体、未利用地和草地的面积变化较小且存在互相转化有关,建筑用地面积的增加会增大产水量,林地和水体面积的增加则会降低产水量,不同年份土地利用类型的相互转移导致了产水量总体变化的不明显。

而极端土地利用变化情景下,产水量的变化情况见表4。由表4可知,与2019年真实情景对比,只有当农田全部转化为建设用地(情景5)时,千岛湖流域产水总量增加了0.13×108m3,增幅为0.23%,其余情景产水总量呈减少趋势。各情景变化量大小依次为情景6(0.17×108m3,减幅0.31%)>情景7(0.11×108m3,减幅0.19%)>情景8(0.04×108m3,减幅0.08%)。主要原因在于,建设用地是流域平均产水量最高的地类,其面积的扩张使得流域下垫面不透水层增加,降水量更多地转化为产水量。因此其余土地利用类型向建设用地的转化,均会导致流域产水总量的增加。同理可得,建设用地转化为林地(情景7)和农田(情景8),均使得流域产水总量减少,其中情景7的变化量显著高于情景8,主要原因在于农田的平均产水量较高于林地,因此情景8的产水总量变化量更低。此结论也可由情景6得出,农田向林地的转化使得流域产水总量降低了0.31%,是所有情景中的最大变化量,说明林地面积的增加会减少流域产水总量。因此,流域建设用地的扩张有利于提高产水总量,退耕还林等植被恢复措施会减少流域产水总量,但有利于流域水土保持,在未来土地利用规划中,适当规划林地、农田和林地的面积有利于实现流域水资源的合理配置。

表4 千岛湖流域极端土地利用情景下的产水总量模拟结果

2.4 千岛湖流域产水格局影响因子分析

千岛湖流域产水格局影响因子探测结果见表5。

表5 产水量空间格局影响因子解释度(q)探测结果

由因子探测结果(表5)可知,各环境因子的解释度q值为实际蒸散量>土地利用>降水量>潜在蒸散量>坡向>坡度>海拔,影响因子在不同程度上解释了产水量的空间分异,实际蒸散量的q值最高,为0.87,是影响产水量空间分异的第一主导因素,其次为土地利用和降雨,q值分别为0.75和0.30。海拔和坡度的q值最小,说明海拔和坡度对产水量空间格局分布的影响最低。

从因子交互作用探测结果(表6)可知,千岛湖流域各因子间交互作用解释均大于单因子,有双因子增强和非线性增强两种组合结果。其中实际蒸散量与其他6个因子的交互作用对产水量空间分布格局影响较为显著,q值均高于0.8,降水量和实际蒸散量的交互作用最高,q值为0.995 5,土地利用和气候因子的交互作用高于其同地形因子的交互,但均为双因子增强。海拔、坡度和坡向与土地利用及实际蒸散量的交互作用均显著高于与其他因子的交互,表明地形因子通过影响气候因子和人类活动间接对产水量的空间格局产生影响。因此在对流域进行生态环境治理和水资源规划时需要考虑多种环境因子的复合作用影响。

表6 产水量空间格局影响因子交互作用q值探测结果

3 讨 论

1995—2019年,气候变化是千岛湖流域产水量年际变化的主要因素,土地利用变化的影响较低。这与相关研究结论一致[9,26-27],杨洁等[9]和窦攀烽等[26]通过情景分析量化降水对黄河流域和宁波地区产水量的影响,发现降水对产水量的贡献度分别为84%和97.56%,吴健等[28]对东北地区产水量与降水进行相关分析,发现两者呈显著正相关关系。一方面,降水是流域水资源的主要来源;白路遥等[29]在通过多年气候数据探究长江黄河流域气候变化对水资源的影响,发现长江流域降水与水资源量的相关系数达到0.9以上,流域水资源变化与降水量密切相关;另一方面,InVEST模型产水模块不分地表、地下径流和基流,通过在栅格尺度考虑可利用水,土地类型等因素,将栅格单元降水与实际蒸散量的差值作为流域产水量。因此本研究模拟的产水量变化主要受到降水变化的影响。空间上千岛湖流域的分布格局变化不大,千岛湖湖区两岸以西人口密集区是产水高值区,湖区东部水面则是产水的低值区。这种分布格局与太湖流域和东江湖流域等大型湖泊所在流域的空间分布相似[23,30],可能与降水的空间异质性和下垫面的蒸散能力相关,如水体的蒸散能力通常高于陆地和植被[31],因此在湖区降水相对较低的情况下其产水量也低于其他地类,而人口密集区域植被蒸散能力较弱且降水相对较高,降水多转化为地表径流,故成为产水量的高值区。但这种空间分布格局又与部分研究如白龙江流域[32]较为不同,这可能因为白龙江流域内湖泊面积占比较小且分布较为破碎,地形气候复杂,故产水量空间分布主要受下垫面覆盖的影响,气候影响相对较低。

定量分析气候与土地利用变化对产水量的影响可知,1995—2005年气候变化对产水量的贡献度为100.50%,土地利用变化的影响为-0.50%,这个结果与窦攀烽等[26]在宁波地区的研究结论相似,气候变化对流域产水服务的影响显著高于土地利用变化。进一步采用极端土地利用变化情景模拟分析发现,建设用地面积的增加有利于提高流域产水总量,而农田转化为林地会导致流域产水总量减少,说明植被恢复会降低流域产水总量,森林面积的增加一定程度上会导致流域产水量的减少,但是从水土保持角度而言,植被恢复使得更多的水留在了森林,有利于流域水源涵养能力的提升。不同土地利用类型的相互转化会导致产水量增加或减少,使得最终量化结果不显著。因此在建设用地产水量难以利用的前提下[27],需综合考虑流域农田、林地及其他土地利用类型的转化,以实现流域水资源可持续发展。

产水量是区域多种因子共同作用的结果,本研究通过地理探测器探究影响产水量空间格局分布的因子。分析表明,实际蒸散量和土地利用类型是产水空间异质性的主要驱动因素,多种因子两两交互作用对产水量的空间分布格局均有显著影响。戴尔阜等[10]通过地理探测器探究横断山区产水空间异质性的驱动因子,发现在海拔较低的中低起伏山地,实际蒸散量和土地利用是主要驱动因子,这与本研究地理探测器结果相似。黄治化等[33]和郑续等[34]分别对石羊河流域和疏勒河流域进行地理探测,结果表明不同地理分区产水驱动因子各不相同。产水量的空间分布格局和不同地区的地表覆盖、植被蒸散能力、降水分布以及地形地貌均有关联。千岛湖流域面积相对较小,地势平缓,降水量的空间分布较为均一,故实际蒸散量成为产水空间异质性的主要驱动因子。相关研究表明,实际蒸散量不仅受到降水、气温等气候因子的影响,而且与地表的土地利用覆盖息息相关[35]。因此在未来的城市经济发展中,要注意土地利用的合理规划,充分考虑地形、气候等因素,以减少对区域水供给服务的负面影响。

另外,InVEST模型评价区域产水量具有一定的不确定性,气象数据的插值未考虑到海拔的影响,且研究采用的土壤数据是基于全国尺度,分辨率不高,因此在模型数据获取时进行研究区的野外样品采样和实验获取数据,实现参数的本地化。

综上,1995—2019年,千岛湖流域产水量年际变化与降水量特征高度一致,空间格局变化不大。气候变化对流域产水量的影响显著高于土地利用变化。 水体和林地的平均产水量显著低于建筑用地、农田、草地和未利用地等,森林植被恢复会降低流域产水量,但有利于提升流域水源涵养能力。 实际蒸散量和土地利用变化是产水量空间分布的主要驱动因子;土地利用变化通过改变下垫面蒸散能力影响实际蒸散量,从而影响产水量;地形因子与其他因子交互作用对产水量空间分异有明显作用。