实际蒸散发包括植被蒸腾和水分蒸发, 是指土壤-植被-大气系统中的水汽交换量与过程, 与生态系统服务、生态安全格局息息相关。在全球陆地降水中实际蒸散发的占比高达60%—65%^[1], 其量的变化会直接影响到地表生态系统的水资源供给服务能力；同时实际蒸散发作为连接地表植被气孔行为、水碳交换的纽带, 控制着地表植被生态系统的格局以及净初级生产力等关键生态过程^[2]。在气候变化逐渐加剧和人类活动不断加强的背景下, 评估区域实际蒸散发的响应变化与成因, 对于区域社会—经济—自然复合生态系统的可持续发展有着重要的意义。目前大尺度的实际蒸散发并不能通过观测直接获取, 模型模拟是实际蒸散发数据的主要获取手段。Yan基于冠层导度过程模型分析出了全球陆地实际蒸散发整体上升的线性趋势^[3]；李修仓等使用平流—干旱模型发现在气温日较差和日照时数的下降是珠江流域实际蒸散发减少的主要原因^[4]；刘健通过互补相关模型发现由于太阳辐射和风速的下降, 鄱阳湖流域早期的蒸散发呈现出下降趋势^[5]；刘冲使用WaSSI-C生态水文模型计算出了鄱阳湖流域实际蒸散发1983—2011年间的上升趋势, 并认为植被恢复是主要原因^[6]。综合已有的文献可以看出, 气象要素和地表植被是实际蒸散发的重要影响因素, 在全球气候变化逐渐加强, 生态修复工程大力推行的背景下, 我国各区域实际蒸散发的时空格局也可能会发生改变。

长江流域是我国重要的水资源供给区域, 随着南水北调工程的建设, 其重要性得到了进一步的提升。由于流域内水资源时空分配不均, 长江流域旱涝灾害频发^[7-8], 水资源供给服务能力不稳定。实际蒸散发是流域生态系统水资源供给服务能力的直接影响因素^[1], 分析长江流域实际蒸散发的变化趋势与成因对于流域的水资源管理有着重要的意义。Gao^[9]和Wang^[10]发现20世纪60年代至20世纪末期间, 受到日照时数下降的影响, 长江流域实际蒸散发呈现出下降的趋势。早期关于长江流域实际蒸散发的研究在趋势分析方面较少考虑实际蒸散发演变趋势的阶段特性, 在影响因素分析方面仅分析了气象因素的影响, 忽略了地表环境变化的作用。随着我国退耕还林工程和农业效率的提高, 21世纪以来长江流域地表植被覆盖度快速增长^[11], 同时流域内气候变化仍在进一步加剧, 20世纪90年代流域内出现了气温突增的状况^[12]。在气候和地表环境变化的综合影响下, 近几十年来长江流域实际蒸散发的响应状况目前仍然未知。本文基于水量平衡原理对GLEAM_v3.2a^[13]、MOD16^[14]、GLDAS_Noah2^[15]、ERA_Interim^[16]四种已公开发布的实际蒸散发模型计算结果数据集在长江流域的精度进行了对比验证, 选取精度最高的实际蒸散发数据集, 使用分段线性拟合、逐步回归、相关性分析等方法分析了1981—2017年间长江流域实际蒸散发演变趋势及主导因素, 研究为长江流域的生态服务评估、生态安全保障的相关研究提供科学的依据与参考。

1 研究区域与数据 1.1 研究区域概况

长江发源于青藏高原唐古拉山脉各拉丹冬峰, 全长6397 km, 流域面积180万km², 涉及了中国19个省、自治区、直辖市, 涵盖人口约为4亿。流域内地形地貌复杂多样, 整体上呈现出西北高东南低的三级阶梯状分布态势(图 1)。气候上, 流域由西向东横跨了青藏高寒区、西南热带季风气候区、华中亚热带季风气候区, 多年平均气温约为15℃, 东高西低、南高北低的；多年平均降雨量约为1120 mm, 降水量的空间分布极不均匀, 自东南向西北呈减小趋势。作为中国第一长河, 长江流域内部水系发达, 湖泊众多, 水资源丰富, 流域内的江都水利枢纽和丹江口水库被选为南水北调水利工程东线和中线的水源地, 分别于2013年末和2014年末开始向中国北方水资源匮乏地区调水, 截止2019年年末累计调水300亿m^3[17]。由于水资源分布的极不均匀和极端气候的频繁出现, 长江流域也是我国旱涝灾害频发的地区^[7-8], 有研究显示气候的变化对流域内旱涝灾害的发生有驱动的作用^[18-19]。

1.2 研究数据

本次研究的数据来源分别是英国布里斯托大学(https://www.gleam.eu/)GLEAM_v3.2a、美国蒙大拿大学(http://www.ntsg.umt.edu/project/modis/mod16.php)MOD16、美国国家航空航天局(https://search.earthdata.nasa.gov/search)GLDAS_Noah2、和欧洲中期天气预报中心(https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-interim)ERA_Interim四种实际蒸散发数据集, 所选数据集概况如表 1所示。用于蒸散发数据集精度评价的径流量数据来源于长江泥沙公报(http://www.cjw.gov.cn/z-wzc/bmgb/)。气象站点监测数据下载自中国气象数据网(http://www.nmic.cn/site/article/id/29086.html), 具体包括降水量、气温、风速、相对湿度以及日照时数监测数据, 植被归一化指数数据取自美国国家航空航天局提供的VIP数据集(1981—2014年)和MOD13数据集(2015—2017年)。在进行数据分析之前, 首先在ArcGIS平台上, 对气象站点监测数据进行空间插值, 对实际蒸散发数据和植被归一化指数数据进行拼接、裁剪、重投影、重采样、加和等预处理, 最终将所有数据的坐标与投影均转化为WGS84与Albert等积投影, 时空间分辨率统一为1 a与1 km, 用于进一步的数据分析。

2 研究方法 2.1 水量平衡原理

本文中水量平衡原理主要应用于蒸散发数据集的精度验证和讨论部分长江流域生态系统水资源供给服务的计算, 水量平衡原理的计算公式如下:

(1)

式中, P为流域年降水量, mm；R为流域年径流量, 同时也可以表示流域生态系统的水资源供给服务量^[20], mm；ET为年实际蒸散发, mm；C为流域地面及地下储水量的变化量, mm, 在年度或多年尺度上, C可以忽略不计^[1]。

在蒸散发数据集的对比验证过程中, 本文基于水量平衡原理, 计算了验证时段内四个数据集的平均绝对误差MAE、相对绝对误差MRAE、均方根误差RMSE三项误差指标以及相关系数指标r。误差指标可以描述各蒸散发模拟数据集相对的误差大小, 而r则可表示各数据集所描绘的长江流域实际蒸散发变化趋势与现实状况的吻合程度。四项指标的计算方法在相关文献^[21]中有详细说明, 在此不做赘述。

2.2 流域实际蒸散发演变趋势的分析方法

分段线性回归是Shao等提出的一种线性回归模型, 它能检测出时间序列数据趋势变化的m个突变点, 识别出研究时间范围内m+1个有着不同线性变化趋势的时间段^[22]。本研究同时选用线性回归和分段线性回归法, 分析长江流域实际蒸散发及其影响因素的整体与阶段变化趋势状况。分段线性回归方程如下:

(2)

回归方程中, t表示年份, 取值区间为研究的时间范围；ET_t表示年份为t时流域ET的拟合值；CP_j(Change Point)表示分段线性回归法中检测出的第j个趋势突变点, 0 < j < m；a_j和b_j分别表示第j个分段中, 回归方程的斜率和常数项, a_j可以表征第j个时间段内, 流域实际蒸散发的线性变化趋势。拟合优度R²和显著性检验指标P是两个常用的统计检验指标, 本文采用以上两个指标对趋势分析结果进行检验。

2.3 实际蒸散发变化的成因分析方法

本文采用逐步回归分析和相关性分析的方法, 探究气象与植被因素对流域实际蒸散发趋势的影响。逐步回归分析法可以保证最后所得的回归模型中所有的影响因子都具有显著性, 以建立最优回归方程, 对实际蒸散发的变化进行归因。最优回归方程如下:

(3)

式中, a₀为常数项；a₁, …, a_k为各入选自变量的标准化回归系数, 能够表示了各自变量对因变量变化的相对影响大小；X₁, …, X_k为入选最终回归方程的自变量；Y为因变量。

3 研究结果及分析 3.1 长江流域实际蒸散发数据集的精度对比验证

在区域尺度实际蒸散发的研究中, 时空连续的实际蒸散发数据必不可少, 为了满足这一需求, 许多学者开展了大尺度实际蒸散发模型模拟的研究, 目前为止, 已有多种长时间序列的实际蒸散发数据集被发布, 例如GLEAM_v3.2a、MOD16、GLDAS_Noah2、ERA_Interim、JRA55^[23]等。在众多的数据集中, MOD16数据集因其易获取、空间分辨率高、精度较好的优点, 在国内各流域实际蒸散发的研究中应用较多^[24-27]。考虑到实际蒸散发数据的可获取性和研究结果的可靠性, 本文基于水量平衡原理, 使用2005—2013年间的降水量空间插值数据与长江泥沙公报中的径流量数据, 计算了金沙江流域、嘉陵江流域、乌江流域、汉江流域、洞庭湖流域、鄱阳湖流域的年实际蒸散量, 基于此评价了GLEAM_v3.2a、MOD16、GLDAS_Noah2、ERA_Interim四种蒸散发数据集的精度, 并在长江流域尺度上对比了年实际蒸散量模型模拟值与水量平衡计算值的年际变化状况, 结果如图 2所示。

在4种蒸散发数据集中, ERA-Interim三项误差指标值最大, 其均方根误差RMSE高达201.24 mm, 高于MOD16在全球26个流域的检验结果(RMSE=194 mm)^[28]；GLDAS-Noah2、GLEAM-v3.2a、MOD16数据集的误差较小, RMSE均小于MOD16在全球26个流域的检验结果, 其中MOD16蒸散发数据集的误差最小, 但其对应的相关系数也最低, GLEAM-V3.2a数据集的误差虽然略高于MOD16数据集, 但其相关系数r最高。综合考虑的MAE、MRAE、RMSE以及r四项精度评价指标, 认为4种蒸散发数据集中, GLEAM-v3.2a数据集最能反映出长江流域ET的实际时空格局状况。因此, 本文最终选用GLEAM-v3.2a数据集, 用于长江流域ET时空格局与成因的分析。

3.2 长江流域实际蒸散发的演变趋势

基于GLEAM-V3.2a实际蒸散发数据, 使用线性趋势分析与分段线性趋势分析1981—2017年间长江流域实际蒸散发的变化趋势, 结果如图 3所示。线性趋势分析的结果显示, 在1981—2017年间, 长江流域实际蒸散发呈现出1.16 mm/a的显著上升趋势(P < 0.01, R²=0.39), 而分段线性趋势分析的结果(R²=0.47)表明, 在1981—2017年间, 长江流域实际蒸散发的线性趋势在1997年左右存在突变点, 1981—1997年间长江流域实际蒸散发以0.02 mm/a线性趋势不显著下降(P>0.05), 1998—2017年间长江流域实际蒸散发则以1.94 mm/a的线性趋势显著上升(P < 0.01)。对比线性趋势分析和分段线性趋势分析的拟合优度R², 可以发现分段线性趋势分析的结果能够更加精确的描绘出长江流域实际蒸散发的变化趋势。Gao和Wang的研究^[9-10]也发现早期长江流域内实际蒸散发呈现出下降的线性趋势, 这点与本文的分段线性趋势分析中1981—1997年间的趋势分析结果一致。此外, 为了进一步验证分段线性趋势分析中1998—2017年间的显著上升趋势结果, 本文使用2000—2017年中国水资源公报和中国泥沙公报中的降水量、径流量数据, 基于水量平衡原理计算了长江流域ET的变化趋势(2013年年末南水北调工程正式通水, 但年均约60亿m³的跨流域调水量在长江径流量多年均值中占比不到1%^[17], 故忽略其对ET计算的影的分段线性趋势分析结果(1.94 mm/a)非常接近, 证明了本文趋势分析结果的可靠性。

实际蒸散发的变化趋势在长江流域内部具有一定的空间差异性。以突变点1997年为时间间隔, 计算1981—2017年、1981—1997年、1998—2017年三个时间段内长江子流域实际蒸散发的变化趋势, 结果如图 4、表 2。在1981—1997年间, 11个长江子流域的ET变化趋势均不显著, ET不显著下降的区域主要集中在长江流域西北部(金沙江流域、岷沱江流域、嘉陵江流域、汉江流域)和东部地区(太湖流域、鄱阳湖流域), 下速率在-0.77—-0.03 mm/a之间, ET不显著上升区域主要分布在干流区域、乌江流域以及洞庭湖流域, 上升速率在0.38—1.23 mm/a之间；1998—2017年间, 所有子流域ET均呈现出增长的趋势, 显著增长的区域主要集中在西北部分子流域(岷沱江流域、嘉陵江流域、汉江流域)以及上游干流区和鄱阳湖流域, 其中嘉陵江增长速率最大, 为3.39 mm/a。从整个研究时间范围(1981—2017年)来看, 在11个子流域中, 除了西北部的金沙江流域、岷沱江流域不显著外, 其他子流域实际蒸散发均呈显著增长的趋势, 增长速率的高值区主要集中在中、下游干流区和鄱阳湖流域, 增长速率均超过了1.8 mm/a。

3.3 长江流域气象与植被因素演变趋势

气候因素和下垫面植被的状态是影响实际蒸散发过程的重要因素, 在探究长江流域蒸散发演变趋势成因之前, 有必要预先分析气象要素与植被状况的变化趋势, 以便于理解各影响因素作用于实际蒸散发的过程与机制。基于气象监测站点气象数据的空间插值结果和植被归一化指数数据, 使用线性趋势分析的方法, 分析了长江流域气象与植被要素在1981—2017年、1981—1997年、1998—2017年3个时间段内的变化趋势, 结果如表 3所示。在1981—2017年间, 长江流域内植被归一化指数呈现出显著上升的趋势, 且1998—2017年内的上升速率要显著高于1981—1997年；1981—2017年间气温也显著上升, 但在1981—1997年、1998—2017年两个时段内, 气温上升的趋势并不显著, 这间接说明气温在1997年左右存在着突增的状况, 钟晨晨在研究中也发现了这一现象^[12]；日照时数在整个研究时间范围内保持着下降趋势, 但下降趋势有所减缓；而风速、降水、湿度的变化趋势都具有阶段性, 前期(1981—1997年)长江流域风速和降水量有所下降, 相对湿度有所升高, 但后期(1998—2017年)风速和降水量则呈现出上升的趋势, 相对湿度有所下降。

3.4 长江流域实际蒸散发影响因素分析

气候要素能够改变水汽传输环境、影响能量供给的条件, 进而决定着实际蒸散发可能达到的最大量, 而下垫面上生长的植被则分别作为截流蒸发和蒸腾效应的过程载体, 会进一步影响实际蒸散发量的大小^[1]。近三十多年以来, 长江流域实际蒸散发及其气象、植被影响要素呈现出了明显的阶段性变化, 为了探究不同阶段长江流域实际蒸散发演变主要影响因素的时空演替状况, 本文使用逐步回归分析和皮尔森相关性分析的方法分别在长江流域和子流域的尺度上探究了不同阶段气象与植被要素在蒸散发演变过程中的重要性。

根据逐步回归分析结果(表 4)可以发现, 在长江流域尺度上, 实际蒸散发的主要影响因素具也有阶段性变化。在1981—1997年间, 6个影响要素中仅有日照时数入选最终回归方程, 标准化系数为0.470(P < 0.05), 日照时数的下降是长江流域实际蒸散发下降的主要原因, 日照时数的下降会减少地表蒸散发的能量来源, 进而削弱地表蒸发蒸腾的过程, Wang^[10]在研究中也得出了相似的结论；1998—2017年间, 植被归一化指数NDVI和日照时数入选了逐步回归方程, 标准化回归系数分别为0.768、0.551(P < 0.01), 说明在这一阶段的植被恢复对长江流域的蒸散发过程有明显的正向促进作用。植被的增长会增加蒸散发过程的降雨截流蒸发量与植被蒸腾量, 虽然日照时数的下降对实际蒸散发有抑制效应, 但在植被恢复的主导下, 长江流域实际蒸散发仍呈现了出显著上升的趋势。整体上来看, 1981—2017年间, 植被恢复是长江流域实际蒸散发显著增长的主导因素, 其标准化回归系数最高, 为0.805(P>0.01)。

在子流域尺度上, 本文基于皮尔森相关系数, 分析了各子流域内气象与植被要素在不同阶段与实际蒸散发的相关性, 结果如图 5—7、表 5—7所示。在1981—1997年间, 日照时数在长江大部分子流域内与实际蒸散发呈显著正相关关系, 西北部的岷沱江、嘉陵江、汉江流域以及东南区域的鄱阳湖流域主要受到日照时数下降的影响, 实际蒸散发量有所减少。而在干流区、乌江流域、洞庭湖流域内, 实际蒸散发在日照时数、植被、降水、气温等因素的综合作用下, 呈现出上升的趋势。太湖流域与其他子流域有所区别, 其实际蒸散发仅与植被归一化指数显著相关, 通过分区统计发现太湖流域的植被归一化指数在这一阶段有所下降(张亮在分析长江流域植被覆盖度时空变化时也得出了一致的结论^[29]), 地表植被退化是太湖流域实际蒸散发下降的主要原因。在1998—2017年间, 除了洞庭湖、鄱阳湖、太湖流域外, 长江大部分子流域中植被归一化指数与实际蒸散发显著正相关或正相关系数最大, 植被恢复是长江流域这些地区实际蒸散发上升的重要原因。洞庭湖、鄱阳湖流域实际蒸散发在这一时期仅与相对湿度显著负相关, 太湖流域实际蒸散发则仅与气温显著正相关。此外风速、日照时数在嘉陵江、乌江、上游干流区对蒸散发也有较大的影响。在整个研究时间范围(1981—2017年)上来看, 植被恢复是所有子流域实际蒸散发上升的重要原因, 风速的回升对西部部分子流域实际蒸散发的升高也有较大的正向作用。在大部分子流域中, 气温在1981—1997年、1998—2017年两个时段内与实际蒸散发的相关性并不明显, 但在更长的时段(1981—2017年)上与实际蒸散发显著相关, 这可能是由于气温在1981—1997年、1998—2017年两个时段内变化并不显著, 在1997年左右突增的状况导致的, 说明气温的突增现象在长江流域实际蒸散发的长期升高趋势中仍起到较大的促进作用。

4 讨论

实际蒸散发作为地表水能循环的关键环节已得到了广泛的关注, 但早期的相关研究中^[9-10]较少考虑区域实际蒸散发变化趋势的阶段性。结合早期Gao^[9]和Wang^[10]的研究结果可以发现, 长江流域实际蒸散发的演变趋势具有明显的阶段性, 从20世纪60年代开始的显著下降^[9], 到20世纪末下降趋势有所减缓, 再到进入21世纪后流域实际蒸散发显著上升(图 3)。20世纪60年代以来, 长江流域气温一直具有上升的趋势, 实际蒸散发量在气温升高的促进下理应有所增加, 但日照时数的下降减少了地表蒸散过程的能量来源, 导致20世纪末之前, 长江流域实际蒸散发在气温升高的情况下仍呈现下降趋势^[10]。随着时间的推移, 长江流域日照时数下降趋势有所减缓, 而长江流域植被生态系统由于农业效率的提升和植树造林工程的推行^{[11, 29]}, 质量有了显著的提高(表 2), 受到植被变化的主要影响, 自20世纪末以来, 长江流域实际蒸散发显著增加, 茂盛的植被使得植被蒸腾量和降水截留蒸发量有所增加, 同时气温持续升高、风速由降转升(表 2)也在一定程度上促进了流域实际蒸散发的升高。影响因素对实际蒸散发的作用存在着空间差异性, 这不仅体现在长江流域内部(图 5—7), 同时也存在于不同流域之间。曹文旭等^[30]基于Budyko假设发现气候变化在潮河流域实际蒸散发1961—2009年间的增长趋势中起到主导作用, 贡献率高达80%, 而植被变化的贡献率仅为20%, 与本文相比, 实际蒸散发影响因素分析结果存在明显差异。本文认为这一状况可能与流域间干湿状况的差异以及影响因素变化程度的大小有关, zhang^[31]在研究中发现黄流域的蒸散发大小主要受到水分的限制, 而珠江流域的实际蒸散发则对太阳辐射量变化更为敏感；而张丹等^[32]基于Budyko假设发现, 在降水、潜在蒸散发、地表环境三个要素中, 我国干旱流域对降水最为敏感, 湿润流域对地表环境参数敏感性最强。

实际蒸散发的持续升高会直接影响到长江流域生态系统的水资源供给服务能力, 基于水量平衡原理, 结合GLEAM实际蒸散发和降水量插值数据, 可以计算出长江流域生态系统水资源供给服务量在实际蒸散发显著增加背景下的变化趋势, 结果如表 8所示。21世纪以来, 长江流域的年降水量变化趋势并不显著(表 2), 年实际蒸散发量显著增加(图 3), 导致了流域生态系统水资源供给服务总量呈现出下降的趋势。流域东部地区受到降水量增加的影响, 水资源供给量有所增加, 气候有暖湿化的趋势, 而西部地区水资源供给量的变化趋势与流域整体一致, 有所减少, 区域旱情加剧, 这与黄涛^[33]的研究结果一致。长江流域西部地区是我国生态修复工程实施的主要区域之一, 1997—2006年间云南地区由于植树造林工程的实施, 森林覆盖率由34%增长到了50%^[34], 植树造林能够增强区域生态系统固碳释氧、水土保持等生态服务功能, 但同时也会削弱区域生态系统的水资源供给能力。相较于本地的原生林, 快速生长的人类造林有着更强的蒸散能力^[35], 在区域旱情加剧的背景下, 建议在长江流域西部干旱地区的生态保护与修复过程中, 优先加强对当地原生植被的保护, 同时在植树造林、退耕还林工程中, 选择针叶乔木等蒸腾能力较弱的树种, 以减少植被耗水量, 缓解旱情加剧的趋势。

5 结论

本文对从多个实际蒸散发模型模拟数据集进行了精度对比验证, 然后基于精度最优的数据集, 分析了1981—2017年间长江流域实际蒸散发的趋势与成因。结果表明:

(1) 在GLEAM_v3.2a、MOD16、GLDAS_Noah2、ERA_Interim四种蒸散发数据集中, GlEAM_v3.2a的综合精度最高；

(2) 长江流域实际蒸散发具有先减后增的阶段性变化趋势。在1981—1997年间, 长江流域实际蒸散发以0.02 mm/a的速率不显著减少, 减少的区域主要集中在长江流域西北部和长江流域东部；在1998—2017年间长江流域实际蒸散发以1.94 mm/a的速率显著增加, 所有子流域实际蒸散发均有所增加；整体上, 在1981—2017年间, 长江流域实际蒸散发以1.16 mm/a的速率显著增加；

(3) 1981—1997年间, 长江流域西北部和东部区域的实际蒸散发主要受到太阳辐射量下降的影响, 蒸散发有所下降, 而干流区、乌江以及洞庭湖区域由于受到日照时数、降水量、植被状况等因素的综合影响, 在日照时数下降的情况下, 实际蒸散发仍有所上升；1998—2017年间, 植被恢复是长江流域大部分区域实际蒸散发上升的重要原因, 东南地区的洞庭湖、鄱阳湖流域实际蒸散发主要受到相对湿度的影响, 此外西部部分子流域中, 风速的回升对实际蒸散发的增加也起到了促进作用。气温在长江流域实际蒸散发的短期(1981—1997年、1998—2017年)阶段性变化中与实际蒸散发的相关性并不明显, 但在其长期(1981—2017年)的变化趋势中, 气温仍起到非常重要的作用；

(4) 受到实际蒸散发显著增加的影响, 长江流域生态系统水资源供给服务总量在1998—2017年间有所下降, 流域内部生态系统水资源供给服务量呈现出西部减少、东部增加的趋势；

(5) 生态系统水资源供给服务量的减少会使得长江流域西部缺水地区旱情加剧, 建议在长江流域西部地区加强对本地原生林的保护, 同时在植树造林、退耕还林等生态修复工程中优先选择针叶乔木等蒸腾能力较弱的树种, 以降低植被的耗水量, 缓解区域干旱加剧的趋势。