河北坝上地区绿水时空变化及其驱动因素研究*

2021-06-09石嘉丽张晓龙李红军沈彦军

中国生态农业学报(中英文) 2021年6期

石嘉丽,张晓龙,李红军,沈彦军**

(1.中国科学院遗传与发育生物学研究所农业资源研究中心/中国科学院农业水资源重点实验室/河北省节水农业重点实验室石家庄 050022;2.中国科学院大学北京 100049)

“绿水”的概念最早由瑞典水文学家Falkenmark[1]提出,是指在热力作用和植物根系吸水作用下由降水转化而来的水资源[2],包括绿水流和绿水蓄量。绿水流主要为地表实际蒸散发,满足了世界83%的雨养农业用水需求[3],是植物生长的必要条件。作为水量平衡和热量平衡的纽带,绿水流的时空变化深刻地影响着生态系统的协调均衡。绿水蓄量主要为土壤水分变化量[4],直接影响土壤水分传输、土壤蒸发、植物根系吸水等过程。一般情况下,绿水蓄量的年际变化可忽略不计。开展绿水时空动态变化研究可以为理解植被耗水、水文过程变化及水资源合理利用提供新的思路。

绿水水量计算方法主要包括水量平衡法和水文模型模拟法[5-6]。水量平衡法是通过计算区域降水量与径流量(包括地表径流与地下径流)之差得到区域绿水量[5-7],因其原理易于理解,可操作性强而广受青睐。水文模型法对水文过程的反映详尽,但对驱动数据要求较高,机理相对复杂,在实际应用中存在一定局限性。应用较广的水文模型有LPJ(Lund-Potsdam-Jena)[6]、SWAT(Soil and Water Assessment Tool)[8]、HIMS-VIH(Hydro-Informatic Modeling System-Vegetation Impacts on Hydrology)[9]、MIKE SHE[10]等模型。众多学者已在多种尺度下开展了绿水资源的研究,如从小流域尺度[5]到区域[8]再到全球尺度[6],日尺度[10]到年际尺度[8]再到年代尺度[7],这些研究从不同角度揭示了绿水多时空尺度的变化过程。在气候变化和人类活动对水循环影响日益深刻的背景下,探究绿水变化的驱动因素成为绿水研究的热点之一。Zang等[11]定量衡量了土地利用变化对黑河流域蓝绿水转化的影响,发现土地利用变化使2.06亿m3的蓝水转化为绿水;赵安周等[12]探讨了气候变化和人类活动对黄河支流渭河流域蓝绿水资源的影响,发现人类活动对绿水量变化影响较小,气候变化是导致研究区绿水量变化的主要原因;冯畅等[13]评估了变化环境对涟水流域蓝绿水的影响,发现在时间变化特征上气候变化使绿水流呈一定的上升趋势而绿水蓄量呈持续下降趋势,在空间分布上绿水资源与土地利用空间分布格局表现出较显著的相关性。这些研究对了解当地绿水变化机理具有重要贡献,但普遍对蒸散发过程考虑不足,气象站点驱动数据在从站点到区域的尺度转换中也存在一定不确定性。

不断发展的遥感技术为绿水研究提供了更加便利的条件,可以避免数据资料不足、尺度转化造成的问题,已成为当前研究绿水流的主要手段。其中MODIS(Moderate Resolution Imaging Spectroradiometer)蒸散产品在全球各区域得到了广泛应用[14-17],由北京师范大学开发的GLASS(Global Land Surface Satellite)蒸散产品[18]在MODIS数据集的基础上,利用贝叶斯方法改进的PM、PT-JPL、MS-PT及半经验彭曼算法计算得到,以通量站观测数据为参考值确定每种算法的权重,其数据精度得到较大提升,具有良好的应用前景。

河北张家口坝上地区为典型的农牧交错地带,属于内流区,肩负着构建京津冀绿色生态屏障的任务,具有防风固沙和涵养水源等生态功能[19]。但该地区生态环境脆弱,处理好水资源开发与环境保护的关系,对实现区域经济、社会、环境协调发展和保障京津冀地区生态环境安全具有重要意义。气候变化背景下,以退耕还林、还草为主体的国家生态建设工程极大地改变了下垫面状况,进而直接影响该地区绿水的时空演变过程。作为水热平衡的重要组分,绿水资源的时空变化规律及其驱动因素尚不明确,开展该研究是科学进行区域水资源管理及合理规划的前提,也是“首都水源涵养功能区和生态环境支撑区”建设的迫切需求。因此,本研究选择河北张家口坝上地区为研究区,基于广泛使用的蒸散产品、气象再分析数据产品、土地利用、地表净辐射等遥感产品,分析该地区绿水时空变化特征及其主要驱动因素,以期为坝上地区“两区”建设、水资源高效利用提供科学依据。

1 研究区概况与数据来源

1.1 研究区概况

“坝上”由地形陡然升高而形成,位于华北平原向内蒙古高原过渡的地区。张家口坝上地区包括张北县、沽源县、康保县和尚义县,位于113.81°～116.06°E,40.73°～42.17°N,海拔多在1200～1800 m,南高北低(图1)。该区域为东亚大陆性季风气候,是半干旱区向干旱区的过渡带,年平均气温−0.3～3.5 ℃,年均降水量400 mm左右[20],集中在6—8月。特殊的地形和气候条件使该地区生态系统相对脆弱,自我修复能力差。主要土地利用类型为耕地和草地,产流微弱,绿水是降水量消耗的主要途径。近20年来高耗水作物种植规模不断扩大,过度的灌溉活动使地表水断流、湿地湖泊萎缩,水资源形势愈加严峻。

1.2 数据来源与处理

实际蒸散发(ET)数据与总第一性生产力(GPP)数据来自于北京师范大学开发完成的GLASS数据集[18]。气温、降雨数据来源于中国科学院青藏高原研究所开发的中国区域地面气象要素驱动数据集[21-22]。净辐射(Rn)数据来源于张晓龙[23]估算的1982—2017年中国大陆区域逐日Rn。土壤水分数据来自中国科学院青藏高原研究所开发的基于微波数据同化的中国土壤水分数据集[24-27]。土地利用数据来源于中国科学院资源环境与数据中心1∶100万土地利用数据产品。数据详细信息见表1。

ET数据为8 d数据,对其进行格式、投影转化,影像裁切处理,再将8 d数据线性插值得到每日数据,累加求和得到多年月均数据及年尺度数据。年均降雨量、Rn及气温由日值数据累加及求均值得到,最后将以上数据重采样统一为0.01°。将2000年、2010年、2015年3期土地利用数据重分类为耕地、林地、草地、建设用地、水体、未利用地6类,并结合野外实测数据对错分类别进行修正,将栅格格式转化为矢量后分别对绿水量进行提取。本研究除土壤水分数据以外,其他所有研究数据研究时段均为2001—2015年,受数据获取限制,土壤水分数据研究时段仅为2002—2011年。

表1 研究所用数据类型、来源及时空分辨率Table 1 Data source and spatio-temporal resolution used in the study

2 研究方法

2.1 绿水蓄量计算方法

绿水蓄量即土壤水分含量变化量,计算公式如下[28]:

式中:ΔW为土壤水分含量变化量即绿水蓄量(mm);h为土层厚度(cm);Δθ为土壤体积含水量变化量(%)。土壤水分含量变化量分为0～5 cm、5～20 cm、20～100 cm 3层计算。

2.2 贡献率计算方法

采用多元回归分析法计算各因子对绿水变化的贡献率。为衡量各因子对绿水量相对影响程度,需对研究数据进行标准化处理,使之介于0～1。对标准化后的数据进行多元回归分析,计算公式如下[29]:

式中:Ys为因变量变化标准值,X1s、X2s、X3s…为自变量标准化值,a、b、c为各因子标准化后回归系数,η1为X1变化对于Y变化的相对贡献率。

2.3 评估指标

用最小二乘法基于像元拟合得到线性变化趋势线,趋势线的斜率即为线性变化系数θslope,以此评估绿水量及各影响因子的变化趋势。在IDL的环境实现计算,计算公式为[30]:

式中:θslope为线性变化系数,i为从2001年开始的第i年,n为研究年数15,Yi为第i年的变量值。利用绿水量及各影响因子年值序列和时间序列的相关关系来判断年际间变化的趋势,θslope＜0表示变化趋势为减少,θslope＞0表示变化趋势为增加。

利用皮尔逊相关系数衡量绿水与降水、温度、Rn的变化相关性,计算公式为[31]:

式中:r为X、Y的相关系数,Xi为第i年的影响因子值,Yi为第i年的绿水量(mm),N为研究年数15。r取值在−1到1之间:r＞0,表示两要素正相关;r＜0,表示两要素负相关。r的绝对值越大,表明二者的相关性越强。

3 结果与分析

3.1 绿水量时间变化特征

本项研究时间尺度较短,土壤水分变化量微乎其微,故年均绿水量计算时仅将绿水流即ET纳入考虑,在多年月均绿水量计算时则同时考虑了多年月均绿水流和绿水蓄量。坝上4县绿水量的年际变化趋势如图2所示,在2001—2015年间,坝上4县绿水量的变化总体呈不显著(P=0.47)的减小趋势,多年平均值为371.11 mm,其中2003—2006年、2008年、2010年、2012年、2013年共8个年份高于平均水平,其余7个年份绿水量低于平均水平,在2003—2007年间连续下降,2007—2015年间绿水量波动起伏较大。波动最为突出的年份为2008年和2009年,相对变化率分别为16.5%和−16.5%。最大值出现在2003年,绿水量为415.34 mm;最小值出现在2009年,绿水量为322.35 mm,年最大绿水量与最小绿水量差值为92.99 mm。

对多年月均绿水量进行空间统计得到多年月均绿水量分布图(图3),绿水量的多年月均值介于2.5～121.8 mm,月份之间差异较大,呈现周期性的单峰变化趋势,峰值出现在7、8月。多年月均绿水量的季节性差异很强,春季(3—5月)绿水总量介于70.64～121.54 mm,夏季(6—8月)为141.15～324.94 mm;秋季(9—11月)绿水总量为50.88～96.56 mm,冬季(12—2月)仅为14.42～27.52 mm,整体表现为夏季＞春季＞秋季＞冬季,坝上4县绿水量的变化基本与水热同期。

3.2 绿水量空间变化特征

张家口坝上4县2001—2015年平均绿水量为279.02～540.50 mm,且具有显著的空间异质性(图4a),由西北至东南呈现很明显的递增规律。其中沽源县的绿水量最高,多年平均值为442.5 mm,该地降水量较大,植被覆盖度高且多为林地,绿水量相应较大;康保县的绿水量最低,多年平均值为343.7 mm,该地降水量低且空气干燥,绿水量较小。像元尺度绿水量变化趋势介于−10.60～2.59 mm·a−1(图4b),仅有0.2%的区域变化显著(P＜0.05)。总体上呈减小趋势,其中有21.2%的区域绿水量变化呈增加趋势,有78.8%的区域呈减少趋势。

3.3 绿水量变化驱动因素分析

3.3.1 气象因子时空变化特征分析

降水是蒸散发水分的主要来源,降水量直接影响ET。2001—2015年坝上4县的平均降水量为348.84～435.38 mm,年均降水量空间上由西北向东南逐渐增加(图5a),与多年平均绿水量的空间分布有明显的递变一致性。其中康保县的年均降雨量最小,沽源县最大。2001—2015年降水量整体呈不显著的增加趋势(图5d),研究区东部降水增加较为缓慢,西部增加较快,高值区集中分布于尚义县和张北县西部。气温是反映气候变化的另一关键因子,也是影响绿水量的另一重要气象要素。坝上4县2001—2015年的年均温为0.85～7.67 ℃(图5b),由东北向西南逐渐减小,尚义县气温最低,沽源县气温最高。气温的年际变化幅度较小(图5e),θslope的绝对值最大为0.11,整体呈不显著的增加趋势(P=0.62)。其中有68.5%的区域呈增加趋势,31.5%呈减小趋势。研究区Rn的多年平均值为2712.30～3097.22 MJ·m−2,随海拔升高由北向南不断增加(图5c)。Rn年际波动较大,如图5f所示,变化趋势自西向东不断增加,在西部呈减少趋势,在东部呈增加趋势。其中有68.9%的区域呈减小趋势,31.1%的区域呈增加趋势,区域整体上呈不显著的下降趋势(P=0.65)。

3.3.2 GPP、土地利用时空变化分析

水循环、碳循环是陆地表层耦合的两个基本生物物理过程[32]。GPP与绿水分别作为水、碳循环的重要产物,二者的变化密切相关。如图6a所示,GPP的多年平均值为140.29 g(C)·m−2,2001—2015年间变化幅度较大,最小值为61.16 g(C)·m−2,最大值为255.37 g(C)·m−2。空间上由西北至东南逐渐增加,高值区集中在沽源县。变化趋势(图6b)也由西北向东南逐渐增加,整体呈不显著的增加趋势(P=0.42)。其中31.3%呈减少趋势,68.7%呈增加趋势,增加最快的区域多位于沽源县中部和尚义县南部。尚义县南部年均GPP较低但增长速度较快,这与该区域林地面积扩大有密切关系。

土地利用变化直观地反映了人类活动对下垫面的作用,是影响绿水变化的重要原因。因土地利用在短时间内变化较小且土地利用产品为多年合成产品,故此处基于2000年、2010年、2015年土地利用数据进行分析,并将原始产品土地利用类型合并为耕地、林地、草地、水体、建设用地和未利用地6类,对原始数据中错分的水体部分基于野外调查数据进行更正(图7)。研究区内不同土地利用类型面积为耕地＞草地＞林地＞未利用地＞建设用地＞水体(表2),2001—2015年耕地、草地、水体面积有所缩减,林地面积连续增加,建设用地面积不断扩张,未利用地面积未发生较大变化。林地的增长面积小于耕地、草地、水体的缩减面积,耕地、草地被建设用地侵占,湖泊干涸演变为未利用地,这在一定程度上使得ET总体上呈减少趋势。

表2 2001年、2010年、2015年张家口市坝上4县各土地利用类型面积占比Table 2 Proportions of different land use types of the four countiesof Zhangjiakou in Bashang area in 2001,2010 and 2015 %

3.3.3 绿水量对主要驱动因子的响应关系

水热条件的变化极大地影响着蒸散,结合绿水量年内变化与水热同期的现象,本文选取了降水、温度、Rn3个气象因子探讨绿水量对气候变化的响应,以降水表征水分条件、以气温和Rn表征热力条件。同时将GPP纳入主要驱动因子考虑范畴,以GPP反映研究区内绿色植物蒸散能力。

基于IDL分别计算各因子与绿水量的相关系数,在像元尺度上分别分析绿水和降水、温度、Rn、GPP的相关性。2001—2015年各影响因子与绿水的相关系数分布如图8所示。坝上4县降水量与绿水量相关系数(图8a)空间分布不均,介于−0.04到0.68之间,显著相关区域仅占13.9%(P＜0.05)。其中呈正相关关系的区域占研究区面积的98%以上,尤其以康保县西北部及尚义县南部的小部分区域最为明显,低值区多位于张北县西部和沽源县东南部。研究区内温度与绿水整体上呈负相关关系(图8b),相关系数最大值为0.18,最小值为−0.76,28.1%的区域呈显著相关(P＜0.05)。有大于99%的区域呈负相关,仅有极少区域呈正相关。相关性较强的区域集中于沽源县。Rn与绿水相关性较弱,仅有0.02%的像元呈显著相关(P＜0.05),相关系数介于−0.62到0.43之间(图8c),总体上与绿水量呈现负相关的关系,只存在15.9%的正相关区域。负相关相关性最强的区域分布于研究区南部边界,正相关区域集中分布于康保县西部、尚义县北部以及沽源县东北部的小部分区域。如图8d所示,GPP与绿水的正相关性很强,平均相关系数为0.63,其中75.0%的区域呈显著正相关(P＜0.05),相关系数最大值高达0.96,只有极小部分区域呈负相关关系。空间分布自西向东递减,其中尚义县的GPP与绿水的正相关关系最为明显,平均相关系数达0.74。

研究区内绿水变化受多种因素共同作用,故以标准化后的气温、降水、Rn和GPP为自变量,以绿水量为因变量进行多元回归分析,定量计算各影响因子对绿水变化的贡献。结果表明绿水量与气温、Rn呈负相关关系,与降水、GPP呈正相关关系,与相关分析结果一致。对绿水的影响程度为GPP＞气温＞降水＞Rn,其中GPP对绿水量变化的贡献率最大,为45%,气温的贡献率为28%,降水的贡献率为16%,Rn对绿水量变化的贡献率最小(11%)。

3.3.4 绿水量与土地利用的关系

提取2001年、2010年、2015年各土地利用类型实际绿水量,结果如图9所示。各年份绿水量均表现为林地＞草地＞耕地＞建设用地。未利用地在研究区内多为干涸湖泊、沼泽等地类,且分布较为分散,受来自周围林地、水体等高绿水量区水汽的影响,绿水量也较高。绿水量高值区的土地利用类型多为草地,草地的植被覆盖度较高且该区域降水量较大,ET相应较大;林地集中分布于绿水量的增长区,这可能由植树造林工程和总初级生产力增加造成。

4 讨论

坝上地区绿水量整体呈减小趋势,年内变化呈现较强的季节性,基本与区域水热同期。吴桂平等[33]在鄱阳湖流域、杨秀芹等[34]在淮河流域对实际蒸散发的同期研究中也得到了相同的季节变化特征。绿水的时间变化主要受气温、降水、净辐射等水热条件的影响。坝上4县降水、温度、GPP整体上呈不显著增加趋势,Rn呈不显著下降趋势,但变化趋势存在显著的空间差异,这与其他学者的研究结果一致[23,35-37]。研究发现绿水与降水、GPP呈正相关关系,与气温、Rn呈负相关关系。坝上4县地处干旱与半干旱区,夏季贡献了主要的蒸散量,夏季平均气温约18 ℃,降水量仅为200 mm左右,热量充足而降水不足,温度升高会使植被受水分胁迫,降水较气温对该地区绿水变化的限制性更强,这可能是造成绿水与降水呈正相关而与温度呈负相关的原因。绿水与降水呈正相关关系得到了广泛认同[38-39],温度与绿水的负相关关系在其他流域也有类似表现。如在干旱区石羊河流域,Xiong等[40]研究发现地表温度与ET呈负相关关系,水分限制是造成这种温度与绿水量显著负相关关系的主要因素,这和本研究结论大体一致。Goroshi等[41]在印度干旱半干旱区、张明明[42]在中国干旱半干旱区的坝上地区、马梓策等[38]在锡林河流域、杨宇娜等[43]在吉林泰周边等的研究中也得到绿水与气温呈负相关的结论。此外,坝上4县绿水量与GPP呈现显著的正相关关系,气温的变化与GPP的年际波动均相反(图10),可能也是气温与绿水量负相关关系形成的重要原因。

气温、降水等因素在时间尺度上影响绿水,土地利用则更直接地影响绿水的空间分布。2000年以来,国家大力实施退耕还林还草等生态建设工程,研究区内植被覆盖变化明显。各土地利用类型绿水量多年平均值差异较小,但绿水量变化趋势受土地利用变化影响显著。研究区内绿水增加区集中于尚义县南部和沽源县,2001—2015年间两地林地面积增加显著,绿水增加区与林地面积增加区相对一致;绿水量剧烈减少区(θslope≤−5)位于张北县中部,安固里淖干涸可能是导致绿水量剧烈减少的重要原因。安塞[44]在京津冀地表蒸散发的研究中发现了ET的空间分布与土地利用类型分布非常一致,与植被覆盖度有良好的相关性;李发鹏等[45]在黄河三角洲源区的研究中强调了陆面蒸散发的空间分布主要受人类活动的影响。这与本研究得出土地利用影响绿水空间分布的结果相对一致。

本研究基于遥感数据产品进行绿水的时空特征及其影响分析,遥感数据产品虽在许多学者的研究中得到验证且具有一定区域适用性,但仍存在一定不确定性。其中,绿水数据分辨率为0.01°,土地利用数据为30 m,二者空间分辨率的差异可能是造成各土地利用类型年均绿水量差异不明显的重要原因。此外,绿水变化还受其他气候因子(如风速、日照时数等)作用的影响,且各气候因子之间也存在复杂的相互作用,多因素水热耦合对绿水时空演变的影响机理还有待进一步研究。