生态学报  2018, Vol. 38 Issue (24): 8835-8848

文章信息

李汇文, 王世杰, 白晓永, 唐红, 操玥, 王明明, 吴路华.
LI Huiwen, WANG Shijie, BAI Xiaoyong, TANG Hong, CAO Yue, WANG Mingming, WU Luhua.
西南近50年实际蒸散发反演及其时空演变
Inversion and spatiotemporal evolution of actual evapotranspiration in southwest China for the past 50 years
生态学报. 2018, 38(24): 8835-8848
Acta Ecologica Sinica. 2018, 38(24): 8835-8848
http://dx.doi.org/10.5846/stxb201805231129

文章历史

收稿日期: 2018-05-23
修订日期: 2018-09-19
西南近50年实际蒸散发反演及其时空演变
李汇文1,3,4 , 王世杰1,2 , 白晓永1,2 , 唐红3 , 操玥1,4 , 王明明1,3,4 , 吴路华1     
1. 中国科学院地球化学研究所, 环境地球化学国家重点实验室, 贵阳 550081;
2. 中国科学院地球化学研究所, 贵州省科技厅普定喀斯特研究综合试验站, 安顺 562100;
3. 中国科学院地球化学研究所, 月球与行星科学研究中心, 贵阳 550081;
4. 中国科学院大学, 北京 100049
摘要: 利用CRU 4.0及GLDAS Noah 2.1数据集,采用随机森林算法对1966年至2016年中国西南陆面月尺度实际蒸散发(ETa)进行逐像元反演,结合袋外误差均方值(MSEOOB)、解释方差百分数(PVE)和均方根误差(RMSE)评价模型以及与其他典型数据集对比的方法对模型和反演结果进行精度评价,在对中国西南ETa的空间格局及时空演变特征进行分析的基础上,利用因子置换重要性评价模型(PIM)对特征因子进行重要性评价。结果表明:(1)MSEOOB均值为4.14,标准偏差仅为3.73,PVE均值为99.36%,标准偏差仅为0.33,模型基于2000年至2016年月尺度拟合结果的RMSE均值仅为1.04 mm/月,标准偏差为0.52,反演结果与GLDAS 2.1、2.0及MOD16数据的R2分别为0.99、0.89、0.95,总体而言模型及拟合结果可信度和精度较高;(2)西南地区ETa整体上表现出随着纬度的降低而增加的特征,从西北高原地区向东南沿海区域逐步增加,不同季节上西南的ETa空间分布差异较为明显,从春季到夏季先呈现出由东南向西北逐步增加的态势,夏季到冬季则呈现出从西北向东南减弱的特征,在每年的7、8月份左右各区域的ETa达到最大值,在1、2月份左右为最低值,并呈现起伏的周期特征;(3)以横断山脉为分界,横断山脉以南的丰水区的ETa主要受云覆盖百分数、月均气温日较差与月均日最高温共同驱动,而横断山脉以北的少水区域主要受云覆盖百分数、月霜日频率与月均水汽压共同驱动,而无论是在丰水区还是少水区,云覆盖百分数都是所有因素中最主要的驱动因子。
关键词: 实际蒸散发     反演     随机森林     时空演变     特征因子    
Inversion and spatiotemporal evolution of actual evapotranspiration in southwest China for the past 50 years
LI Huiwen 1,3,4, WANG Shijie 1,2, BAI Xiaoyong 1,2, TANG Hong 3, CAO Yue 1,4, WANG Mingming 1,3,4, WU Luhua 1     
1. State Key Laboratory of Environmental Geochemistry, Institute of Geochemistry, Chinese Academy of Sciences, Guiyang 550081, China;
2. Puding Comprehensive Karst Research and Experimental Station, Institute of Geochemistry, Chinese Academy of Sciences and Science and Technology Department of Guizhou Province, Anshun 562100, China;
3. Center for Lunar and Planetary Sciences, Institute of Geochemistry, Chinese Academy of Sciences, Guiyang 550081, China;
4. Univeersity of Chinese Academy of Sciences, Beijing 100049, China
Abstract: In this study, we combined the CRU 4.0 and GLDAS Noah 2.1 datasets and used the random forest algorithm (RF) to calculate the monthly actual evapotranspiration (ETa) in southwestern China between 1966 and 2016. The accuracy of our model and the inversion results were evaluated by the MSEOOB, PVE, and RMSE indices and a comparison with other typical datasets. Based on this, the spatiotemporal distribution and evolution pattern of ETa over multiple time scales were fully discussed. We extended our research by using the PIM model to evaluate the importance of the feature factors in each pixel. We obtained the following results:1) the mean value and the standard deviation (Std Dev) of MSEOOB were 4.14 and 3.73, and those of PVE were 99.36% and 0.33; in addition, the mean RMSE of the monthly inversion results for 2000 to 2016 was 1.04 mm per month and the corresponding Std Dev was 0.52; moreover, the R2 values for the inversion results of ETa from GLDAS 2.1, GLDAS 2, and MOD16 were 0.99, 0.89, and 0.95, respectively. All these evaluation indices illustrated the credibility and precision of the model and the inversion results were sufficiently high. 2) Our inversion results indicated that ETa increased with a decrease in latitude and gradually increased from the northwest plateau to the southeast coastal area; in addition, the spatial distribution patterns in southwestern China in different seasons were quite different:from spring to summer, high ETa expanded from the southeast to the northwest; in contrast, in winter, ETa diminished remarkably from northwest to southeast. The maximum value was reached around July every year and had the lowest value was reached around February, which showed periodic characteristics with fluctuations. 3) We found that the Hengduan Mountains were the boundary of the driving factors for ETa in arid and humid regions. The ETa in the humid regions south of the mountains was jointly driven by the cloud cover percentage (CCP), diurnal temperature range (DTR), and monthly average daily maximum temperature (TMX). In contrast, ETa in arid regions north of the mountains was mainly affected by CCP, frost day frequency (FDF), and vapor pressure (VAP). Notability, CCP is the most important driving factor weather in both humid and arid areas.
Key Words: actual evapotranspiration     inversion     random forests     spatiotemporal evolution     feature factors    

中国西南地形地貌极其复杂, 其生态系统多样而复杂, 该区域是研究地表复杂过程与生态系统演变规律的关键区域, 同时, 也是我国非常重要的生态屏障区。因此, 开展西南生态安全格局形成机制、演变机理等方面的研究具有典型性、代表性与紧迫性[1]。实际蒸散发(Actual Evapotranspiration, ETa)包含水面、土壤蒸发和植被蒸腾3个方面, 是全球水循环系统中影响着水体、土壤、植被和大气之间水交换过程的重要媒介, 是影响全球生态系统及全球变化的重要因素[2], 在区域乃至全球生态安全格局的形成及演变过程中扮演着不可或缺的角色。同时, 蒸散发是衡量生态系统活力的重要指标之一, 也是探讨区域生态系统空间格局组成的主要参数之一, 是区域及全球生态环境循环系统中不可忽视的组成部分, 在多类研究中都扮演着极其重要的角色[3-4]。因此, 要想全面了解西南地区水循环过程、生态系统能量转换以及生态安全格局演变机制, 就必须对长时间尺度的ETa有充分的掌握[5], 特别是对其时空演变特征和驱动因素进行系统的分析。

然而实际蒸散发不能在大尺度上(国家、区域、洲等)被直接观测[6]。因此到目前为止, 针对中国西南月尺度下长时间跨度, 特别是过去近50年的ETa数据非常欠缺, 这就使得全面解析西南地区生态安全格局及其演变的相关研究面临着数据缺失的困境。因此, 本研究拟从反演中国西南近50年来的ETa数据出发, 对其空间格局、时空演变规律以及影响因子空间特征进行探讨, 以期为西南生态安全格局的形成机制及演变机理相关研究提供基础数据和理论支撑。本文充分利用东英吉利大学气候研究中心(Climatic Research Unit, University of East Anglia, CRU)提供的1966年至2016年月尺度下的9类气象水文数据及GLDAS Noah 2.1同化数据集提供的2000年至2016年月尺度下的ETa数据, 结合机器学习中拥有较高可靠性和鲁棒性的随机森林算法(RF)对西南地区陆面月尺度ETa进行逐像元反演, 结合解释方差百分数(PVE)、袋外误差均方值(MSEOOB)和均方根误差(RMSE)评价方法以及与其他典型数据集对比的方法对模型和反演结果进行精度评价, 并对中国西南ETa的空间格局及时空演变特征进行分析。在此基础上, 利用因子置换重要性评价模型(PIM)对特征因子进行重要性评价, 得到基于像元的影响因子空间分布特征。本研究不仅为研究西南生态安全格局及其形成机理提供了长时间尺度的基础数据支撑, 也为生态安全演变特征等方面的研究提供了对气候变化的响应等方面的参考。

1 数据和方法 1.1 研究区

为全面解析中国西南ETa时空演变特征, 本研究以西南地区(21°06′—39°20′N, 78°22′—112°02′E)为研究对象, 包括西藏、青海、四川、重庆、贵州、广西、云南。研究区内涵盖丘陵、高原、盆地、山地等复杂地貌, 是我国山地最为集中的区域之一, 这也是西南地区生态安全格局复杂而重要的原因之一。研究区内涵盖了13种气候类型, 又以南部的南亚热带季风雨林气候、东部的中亚热带季风性常绿阔叶林气候以及西部的亚寒带寒带草原气候为主(图 1), 研究区内分布有8个主要的水文分区[7](图 1), 年均温在0—24℃之间, 年降雨量在600—2300 mm范围内[8]。本研究行政区划来源于国家科技基础条件平台-国家地球系统科学数据共享服务平台(http://www.geodata.cn)。

图 1 研究区内气候带和水文区划 Fig. 1 Spatial distribution of climate types and hydrologic regionalization in study area 图中代码对应的具体气候及水文区划类型名称请查看表 1

表 1 研究区气候带与水文区划代码对应类型 Table 1 Climate types and hydrologic regionalization of the corresponding code in Fig. 1
气候类型代码
Climate code
气候带类型
Climate type
水文区划代码
Hydrologic regionalization code
水文区划
Hydrologic regionalization
9 暖温带大陆性草原气候 秦巴、大别山北亚热带多水区
10 暖温带大陆性荒漠气候 东南亚热带、热带丰水地区
11 北亚热带季风性落叶阔叶、常绿阔叶林气候 西南亚热带、热带多水区
12 中亚热带季风性常绿阔叶林气候 滇西、藏东南亚热带、热带丰水区
13 南亚热带季风雨林气候 西北山地中温带、亚寒带、寒带平水、少水地区
14 热带雨林、季风雨林气候 西北盆地温带、暖温带干涸地区
15 温带落叶阔叶林气候 青藏高原东部和西南部温带平水地区
16 温带森林草原气候 羌塘高原亚寒带、寒带少水地区
17 温带草原气候
18 温带荒漠气候
19 亚寒带草原气候
20 寒带草原气候
21 寒带荒漠气候
1.2 数据来源与处理

本研究使用到的1966年至2016年月尺度下的9类气象水文特征因子数据来自于东英吉利大学气候研究中心CRU发布的TS 4.0数据集[9]。包括月降雨量(MPR)、云覆盖百分数(CCP)、月均气温日较差(DTR)、月霜日频率(FDF)、月均温(MAT)、月均日最低温(TMN)、月均日最高温(TMX)、月降雨频率(WDF)、月均水汽压(VAP), 数据集的空间分辨率为0.5°(表 2)。CRU数据在全球多类研究中得到了广泛的应用, 其精度和适用性也在多项研究中得到了验证[10-12]

表 2 使用到的数据的属性信息 Table 2 Attribute information of the data used
数据名称
Name
数据简称
Abbreviation
数据单位
Units
空间分辨率
Spatial resolution
时间跨度
Time span
GLDAS Noah 2.1蒸散发ETa from GLDAS Noah 2.1 ETa mm/月 0.25°×0.25° 2000.1—2016.12
月降雨量Monthly precipitation MPR mm/月
云覆盖百分数Monthly average cloud cover percentage CCP %/月
气温日较差Monthly average diurnal temperature range DTR ℃/月
月霜日频率Monthly frost day frequency FDF d/月
月均温Monthly average temperature MAT
月均日最低温Monthly average daily minimal temperature TMN 0.5°×0.5° 1966.1—2016.12
月均日最高温Monthly average daily maximal temperature TMX
月降雨频率Monthly wet day frequency WDF d/月
月均水汽压Monthly average vapor pressure VAP H-Pa/月

为计算西南地区近50年月尺度下的实际蒸散发数据, 本研究以2000年1月至2016年12月GLDAS Noah 2.1西南地区实际蒸散发数据为基准[13]。GLDAS结合卫星和地面观测资料, 采用陆面模型(LSM)和数据同化技术为相关研究提供了多类型的高精度数据集, 已成为全球变化与水循环等研究的重要数据源之一[13]。GLDAS Noah 2.1版本数据在1.0和2.0版本基础上, 充分利用全球资料同化系统(Global Data Assimilation System, GDAS)发布的大气分析资料、分类的全球降水气候学项目(the disaggregated Global Precipitation Climatology Project, GPCP)发布的降雨数据和美国空军气象局农业气象模拟系统(the Air Force Weather Agency′s AGRicultural METeorological modeling system, AGRMET)提供的辐射数据集以及修正的地表辐射收支数据集(Surface Radiation Budget, SRB), 对前期版本中主要存在的北半球短波下行通量异常和降雨数据异常问题导致的结果偏差较大等问题进行了修正[14]。该数据集较高的精度使得其在相关研究中得到了肯定和广泛的应用[15-16]

1.3 方法 1.3.1 基于随机森林算法的ETa回归模型构建

随机森林算法(RF)[17]利用bootsrap重抽样方法[18]从原始样本中抽取多个样本, 对每个bootsrap样本进行决策树建模, 然后组合多棵决策树的预测, 通过投票或者计算每棵决策树输出结果的均值得到最终的分类或回归预测结果。大量理论和实践研究都证明了RF算法具有很高的精度和稳定性, 对异常值和噪声具有很好的容忍度, 且不容易出现过拟合[19], 在机器学习、数据挖掘等领域被给予了很高的评价[20-21]。Cernadas等对179种主流的分类算法进行了测试评价, 其结果显示随机森林算法是这些算法中表现最优异的[22]

首先以独立同分布的随机向量Θt组成树形结构分类器{h(X, Θt), t=1, …, N}, 从本研究中影响ETa的9类特征数据集中随机有放回地抽取定量的训练样本, 对于每个样本集进行决策树建模, 然后组合多棵决策树的预测, 通过计算t棵决策树{h(X, Θt)}的预测均值得到最终的ETa预测值。对于待测样本集合X(ETa特征因子集), 其对应样本最终的结果(Yk)按下式计算:

(1)

式中, Yk表示随机森林算法对第k个待测样本特征因子集的计算结果, arg ave(·)表示取均值, ht(·)表示第t棵决策回归树的计算结果, xk|(·)表示第k个待测ETa对应的输入样本(该待测ETa对应的全部特征因子集), 第t棵决策树基于xk的各影响因子特征参数, 结合该棵决策树对应的回归拟合函数ft计算得到对应的ETa。

在自变量X(即各特征因子)和因变量Y(即ETa)组成的集合中, 随机抽样出的训练集(Xi, Yi)相互独立, 随机森林算法整体的预测结果精度用回归预测向量h(X)的泛化误差均方值来描述。其定义为:

(2)

式中, Y为训练样本集X对应的实际值, h(X)表示模型对样本集X的预测值。EX, Y(·)2算子表示:

(3)

为评估计算结果的精度, 模型首先采用可靠的精度估计方法, 即基于式(2), 模型利用未参与决策森林构建的袋外样本数据(out of bag, OOB), 对袋外数据的预测均方误差(Mean Squared Error of OOB, MSEOOB)进行计算:

(4)

式中, yk表示袋外数据中第k个ETa, 表示模型针对该ETa的预测值(所有回归决策树对该观测对象的预测均值)。MSEOOB越小表示模型精度越高。其次, 模型还采用了另外一种参数进行精度的评价, 即解释方差百分数(Percentage of Variance Explained, PVE), 该参数表示模型能显著解释样本数据的能力, 该值越接近100%表示模型的效果越好, PVE由以下公式计算得到:

(5)

式中, 表示模型针对袋外数据进行预测的预测结果的方差。通过MSEOOB与PVE这2个参数, 可以对构建的随机森林回归模型进行精度评价。

1.3.2 特征因子重要性评价模型

在本研究中, 针对逐像元构建的随机森林回归模型, 采用基于置换法的因子重要性评价模型(the permutation importance measure, PIM)对每个像元的9类因子的重要性进行评价, 最终筛选出每个像元上最重要的影响因子。

对各独立像元的9类因子中的因子Xj进行重要性评价的方法为将该因子的所有值进行随机置换, 该因子的重要性指标则通过由于置换该因子导致的预测结果精度的改变来确定。原因在于如果该因子纯粹是由随机噪声组成, 那么预测结果的精度将不会因为随机置换该因子的值受到影响[23]。因此, 某棵回归树t上因子Xj的重要性VI(t)(Xj)由以下公式定义:

(6)

式中, Bt为未参与第t棵回归树训练的所有OOB样本, L(Tt(xi), yi)为在第i个袋外样本集的因子Xj未置换前第t棵回归树对第i个袋外样本的预测精度, L(Tt(xi, πj, yi)为在因子Xj被随机置换为xi, πj之后的预测精度。回归模型预测精度L(·)由均方根误差决定:

(7)

式中, yi为第i个袋外样本的实际蒸散发, 为第t棵回归数对第i个袋外样本的预测结果, |Bt|为未参与第t棵回归树训练的袋外样本总量。因此, 因子Xj最终的重要性则由森林中所有回归树计算出的重要性的均值表示:

(8)

在本研究中对于每一个像元对应的9类因子进行了逐像元的重要性计算, 最终选择PIM值最大的因子作为该像元上影响ETa的最重要的决定性因子。

1.3.3 RMSE误差分析

为了对模型精度进行进一步的检验, 研究采用了目前使用较广的均方根误差指数(root mean square error, RMSE)对模型回归拟合结果进行逐像元地分析, RMSE值越低说明模型的拟合结果精度越高。RMSE计算公式如下:

(9)

式中, ETaact为实际蒸散发量, ETasim为模型拟合结果, n为样本容量。

2 结果与分析 2.1 实际蒸散发数据集精度评价 2.1.1 模型训练精度评价

通过对研究区内2000年1月至2016年12月的所有GLDAS Noah 2.1 ETa数据和对应像元上的9类气象水文特征因子进行逐像元的RF模型训练, 针对每一个像元都构建了独立的RF模型, 即每个像元都基于204组样本(17年×12月, 每组数据由对应月份的ETa数据和该像元对应的9类气象水文特征因子组成), 每个像元RF模型均由500棵回归决策树组成, 每棵回归决策树由随机选取的80%个样本数据集进行训练得到, 剩下的20%的袋外数据集用于对模型训练过程的精度进行评价。为评价每个像元对应的RF回归模型精度, 首先对各像元对应的随机袋外数据集进行了预测结果泛化误差分析, 模型MSEOOB空间分布如图 2。其次, 为进一步评价各个模型在对应像元上的回归拟合能力, 对每个模型的解释方差百分数PVE进行了计算, 其空间分布如图 2所示。

图 2 基于像元的RF模型袋外样本预测均方误差MSEOOB及解释方差百分数PVE空间分布 Fig. 2 Spatial distribution of MSEOOB and PVE for each pixel of the established RF model

模型整体的MSEOOB均值为4.14, 标准偏差仅为3.73, 说明模型对于袋外数据的拟合预测结果偏离程度非常小, 即模型在所有像元上的训练情况较好, 模型稳定性及可靠性较高, 多类领域的研究也充分证明了RF模型的稳定性和可靠性[23]。此外, 根据MSEOOB的等级百分比分布(图 3)可知, 98.7%的像元对应RF模型的MSEOOB指数都在15.7以内, 即袋外数据预测误差不超过3.9 mm/月, 这说明基于像元的RF模型训练精度比较高。

图 3 基于模型对2000年1月至2016年12月实际蒸散发(ETa)数据进行拟合的RMSE空间分布及MSEOOB、RMS、PVE各等级面积占比柱状图 Fig. 3 Spatial distribution of RMSE for each pixel predicted by the established RF model from January 2000 to December 2016 and the area percentage bar graphs of each classification of MSEOOB, RMSE and PVE 图中ETa:实际蒸散发, actual evapotranspiration;MSEOOB:袋外数据的预测均方误差, Mean Squared Error of data Out Of Bag;RMSE:均方根误差Root Mean Square Error;PVE:解释方差百分数, Percentage of Variance Explained

对于模型PVE指数, 其总体均值为99.36%, 标准偏差仅为0.33, 这同样说明模型整体的拟合性能比较稳定, 模型较强的拟合能力和鲁棒性在多种领域也得到了广泛的验证[24]。从PVE的等级百分比分布(图 3)可知, 有超过98.3%的像元对应RF模型的PVE参数大于98.5%, 这同样说明了模型的泛化能力较强。因此, MSEOOB和PVE 2个参数初步证明了RF模型可以适用于西南地区1966年至2016年的实际蒸散发拟合计算。

2.1.2 模型预测精度评价

为评价模型的预测精度, 首先利用构建的基于像元的RF模型对2000年1月至2016年12月所有原始GLDAS Noah 2.1 ETa数据集的每个像元进行了拟合计算, 采用RMSE指数对每个像元的回归拟合结果进行精度验证(图 3)。

在像元尺度上, RF模型对于2000年至2016年月尺度拟合结果的RMSE指数均值为1.04 mm/月, 标准偏差为0.52, 总体而言模型的拟合精度较高, RMSE的空间分布与MSEOOB的空间特征基本一致。从RMSE的等级百分比分布(图 3)可知, 98.38%的像元对应RF模型的拟合误差低于2.5 mm/月, 这充分说明了构建的RF模型其精度和稳定性的可靠。因此, 基于该构建的模型对研究区1966年至2016年月尺度上的实际蒸散发进行拟合计算是可靠的。利用上述构建的基于像元的RF模型, 对研究区1966年至2016年月尺度下的实际蒸散发数据进行了逐像元的拟合计算, 为充分验证计算结果的准确性, 通过与2000年1月至2016年12月的GLDAS 2.1版本原始ETa数据及2001年1月至2014年12月的MODIS MOD16数据和1966年1月至2010年12月的GLDAS 2.0版本ETa数据[13]进行了月均值的统计对比验证(图 4)。

图 4 1966年至2016年ETa月均值拟合计算结果与GLDAS 2.1版本、GLDAS 2.0版本及MOD16 ETa数据的对比 Fig. 4 Comparison of fitted monthly ETa with data from GLDAS 2.1, GLDAS 2.0 and MOD16 in study area from 1966 to 2016

拟合结果及对比数据集都表现出了一致的季节变化特征, 与生长季物候特征基本一致, 夏秋季节达到峰值, 春冬季节为低谷。在2000年至2016年的数据对比中, 拟合结果与计算基准数据集GLDAS 2.1的重合度非常高, 拟合结果与GLDAS 2.1数据集的R2系数达到了0.99(图 5), 说明模型的可靠性和精度是非常高的, 这也与MSEOOB、RMS、PVE三个参数表征的模型特征一致。在2000年至2014年与MOD16数据集的对比中可见, MOD16数据集在1月左右的低谷区域与拟合结果、GLDAS 2.1、GLDAS 2.0数据之间都存在一定的偏差, 相关研究表明, 在蒸散发过程较弱的11月到次年3月, MOD16 ETa产品在中国流域存在普遍性高估的问题[25], 这可能是研究区内MOD16在1月左右的低谷区域表现出显著的偏高原因, 拟合结果与MOD16数据集的R2系数虽然达到了0.95(图 5), 但从斜率和截距可知, 2套数据在趋势上保持了高度的一致性, 但是在量级上还是存在一定差别。从1966年至2010年区间内拟合结果与GLDAS 2.0数据集的对比情况可知, 拟合结果与GLDAS 2.0数据集基本吻合, 1999年之前在低谷区域拟合结果较GLDAS 2.0数据集略高, 在2000年之后, GLDAS 2.0数据集较拟合结果和2.1版本数据集在峰值区域相对偏低, 这主要是因为GLDAS Noah 2.1版本数据在1.0和2.0版本基础上, 对前期版本中主要存在的北半球短波下行通量异常和降雨数据异常问题导致的结果偏差较大等问题进行了修正[14]。从拟合结果与GLDAS 2.0数据集的散点图可知, 其R2系数超过了0.89(图 5), 说明拟合结果在2000年之前的可信度和精度是较高的。

图 5 1966年至2016年ETa月均值拟合计算结果与GLDAS 2.1版本、GLDAS 2.0版本及MOD16 ETa数据的散点图对比 Fig. 5 Scatter diagram of fitted monthly ETa with data from GLDAS 2.1, GLDAS 2.0 and MOD16 in study area from 1966 to 2016
2.2 实际蒸散发时空演变特征 2.2.1 ETa空间分布格局

ETa对陆地水储量具有巨大的影响[26], 也是影响陆地干旱状况的重要因子[27], 因此明确ETa空间分布格局对于解析区域乃至全球水循环系统及干湿状况在空间上的分布特征具有重要的意义。而对ETa在各季节内的空间分布特征进行了解可以更透彻地明确ETa对气候的调节作用以及对农业等的影响[28]。鉴于此, 研究基于计算得到的研究区内1966年至2016年月尺度ETa数据, 通过计算各个季度的ETa均值以全面了解研究区ETa在空间上的演变特征。由于受到太阳辐射、植被、温度、降水、以及地形地貌等因素的综合影响, 在各个季节上, 研究区ETa整体上表现出随着纬度的降低而增加的特征(图 6)。

图 6 研究区1966年至2016年各季节平均实际蒸散发(ETa)空间分布图 Fig. 6 Spatial distribution of the mean actual evapotranspiration (ETa) in study area at each season from 1966 to 2016 各个空间分布图右侧曲线为各个季节纬度带均值分布

不同季节中国西南的ETa空间分布差异较为明显, 各季节平均ETa的顺序为:夏季(6—8月, 66.3 mm/月)>秋季(9—11月, 37.9 mm/月)>春季(3—5月, 34.3 mm/月)>冬季(12—2月, 21.2 mm/月), 各季节内ETa均值分别占到年内ETa的41.52%, 23.73%, 21.48%, 13.27%。在空间上, ETa从西北高原地区向东南沿海区域逐步增加, 研究期间春季最大平均ETa约为81.3 mm/月, 夏季最大平均ETa在所有季节最高, 达到了123.8 mm/月, 而秋季的最大平均ETa为101.7 mm/月, 冬季的最大平均ETa在所有季节最低, 为75.4 mm/月。在秋季和冬季, ETa高值区域主要分布在云南和广西南部, 在夏季, ETa高值区域几乎覆盖了四川、重庆、广西、云南和广西以及西藏和青海的东南区域。夏季的ETa分布相对其他季节而言其最级都相对更大。在研究期间, ETa随着季节的变化从春季到夏季先呈现出由东南向西北逐步增加的态势, 夏季到冬季则呈现出从西北向东南减弱的特征。然而, 在全年尺度上, 西藏和青海西北部, ETa始终保持较低的值域范围, 这和该区域干旱特性是密切相关的[29]。值得注意的是, ETa对于区域乃至全球的气候变化尤其是干旱极其重要[30], 因此对西南干旱区域, 特别是西北山地中温带、亚寒带、寒带平水、少水地区, 西北盆地温带、暖温带干涸地区, 青藏高原东部和西南部温带平水地区以及羌塘高原亚寒带、寒带少水地区, 充分掌握其ETa的时空特征对评估生态系统和水资源的现状与演变规律, 以及制定适应政策是至关重要的。

2.2.2 ETa时空演变特征

以纬度为基准, 通过对研究区每个月的ETa拟合结果进行纬度带的均值计算, 计算结果如图 7所示, 横坐标代表研究期间各个月份对应的ETa在各个纬度上的均值, 纵坐标为对应的纬度区间。在横向上, 即以时间演变为基准, ETa在每年的7、8月份左右达到最大值, 在1、2月份左右为最低值, 并呈现起伏的周期特征。在纵向上, 即以空间坐标为基准, 同样可以得出与图 5相同的结论, 即ETa随着纬度的降低而增加。从1966年至2016年, ETa最大值呈现出逐步增加的趋势, 而最小值则呈现出减小的趋势。在2000年以前预测结果的纬度均值相对于2000年之后的纬度均值而言, 在春季和冬季像元的纯度相对较低, 而在夏季则表现出了较高的一致性, 这可能是由于CRU部分历史数据的不确定性造成的[31]

图 7 研究区1966年1月至2016年12月实际蒸散发ETa纬度均值分布 Fig. 7 Mean ETa in latitude of each month in study area from January 1966 to December 2016
2.3 特征因子空间格局分析

利用2000年1月—2016年12月每个像元的ETa及该像元上对应的9类因子, 采用PIM方法对各个像元上的9类因子的重要性进行评价, 最终筛选出每个像元上最重要的影响因子, 其空间分布如图 8

图 8 基于像元的最重要影响因子空间分布及其百分比 Fig. 8 Spatial distribution of the most important factor in each pixel and the percentage of each factor 图中影响因子缩写中英文含义请查看表 2

根据研究区最重要因子面积百分数占比情况(图 8), 在整个研究区, CCP是占比最高的影响因子(54.78%), 主要分布在西藏、青海、广西、四川西北地区、贵州东南地区, 其次为DTR(17.67%), 主要分布在四川、重庆、贵州与云南的交接地带, FDF(8.27%)与VAP(7.46%)相差不大, FDF主要集中分布在西藏西部, VAP除了研究区少量的零散分布外, 在西藏、四川和云南交界区域有较为集中的分布;TMX(5.75%)几乎全部集中分布于云南境内, 是影响云南ETa计算最主要的因子, 其他的因子在研究区分布则较为零散。在西藏、青海及广西, 云覆盖百分数是占比最大的因子, 在川渝及贵州, 月均气温日较差是占比最大的影响因子, 同时云覆盖百分数也具有较高的占比, 而在云南, 月均日最高温则是占比最大的因子。在所有区域, 云覆盖百分数都有显著的占比情况, 说明在各个区域利用RF模型对ETa进行计算时, 云覆盖百分数的重要性都是非常高的, 因为云覆盖百分数直接影响到地表及植被冠层所接收到的太阳辐射[32], 从而影响蒸散发的量级, 因此该数据的准确性和可靠性对于模型计算结果的精度具有较高的影响。大量研究表明, 气温日较差与实际蒸散发呈现出正相关关系[33], 是影响蒸散发的主要影响因子之一[34-35]。对于西藏、青海及川渝地区, 特别是川藏交界等干旱区域[36], 月均水汽压在这些区域较为集中地分布, 其原因在于, 对于干旱区域, 月均水汽压亏缺和土壤水分一样重要, 它直接决定了蒸散发的量级[37]。对于较为湿润的云南区域, 月均最大气温直接决定了蒸散发的量级。

为了对研究区内的ETa特征因子分布有更全面而系统地了解, 对研究区内各气候类型与水文区划下的最重要因子面积百分数占比进行了统计, 如表 3表 4

表 3 各气候类型像元尺度最重要因子面积百分数占比 Table 3 Area percentages of the most important factors of different climate types in study area
气候类型
Climate
最重要因子面积百分数占比Area percentages of the most important factors/%
MPR CCP DTR FDF MAT TMN TMX WDF VAP
9 0 100.00 0 0 0 0 0 0 0
10 0 100.00 0 0 0 0 0 0 0
11 0 27.91 60.47 2.33 0 0 0 4.65 4.65
12 2.18 19.64 51.27 1.09 0.36 1.82 11.64 4.36 7.64
13 0 56.57 2.02 1.01 1.01 2.02 28.28 5.05 4.04
14 0 20.00 20.00 0 0 6.67 53.33 0 0
15 0 40.91 38.64 0 0 0 0 6.82 13.64
16 0 60.00 11.00 3.00 0 3.00 0 5.00 18.00
17 0.43 70.39 4.72 15.45 0 0 0.43 1.29 7.30
18 0 95.77 0 1.41 0 0 0 0 2.82
19 0 83.51 1.60 11.17 0 0 0 2.66 1.06
20 6.93 53.47 0 26.73 0 0.99 1.98 0.99 8.91
21 7.02 40.35 7.02 14.04 3.51 1.75 0 7.02 19.30
表中气候类型代码所对应气候类型请查看表 1, 因子类型简写含义请查看表 2

表 4 各水文区划像元尺度最重要因子面积百分数占比 Table 4 Area percentages of the most important factors of different hydrologic regionalization in study area
水文区划
Hydrologic regionalization
最重要因子面积百分数占比Area percentages of the most important factors/%
MPR CCP DTR FDF MAT TMN TMX WDF VAP
0 11.11 88.89 0 0 0 0 0 0
0 100.00 0 0 0 0 0 0 0
2.11 25.61 50.88 0 0 1.05 12.98 2.11 5.26
0 15.12 6.98 2.33 2.33 5.81 36.05 15.12 16.28
0 92.86 0 3.06 0 0 0 0 4.08
0 96.00 0 0 0 0 0 0 4.00
0.47 66.20 12.59 9.32 0 0.70 0 2.56 8.16
4.20 53.78 2.10 23.53 0.84 0.84 1.26 4.20 9.24
表中水文区划代码所对应水文区划类型请查看表 1, 因子类型简写含义请查看表 2

研究区各个气候带内最主要的影响因子具有不同的空间分布特征, 除了北亚热带与中亚热带最重要的特征因子为月均气温日较差(DTR)以及热带雨林、季风雨林气候最重要的特征因子为月均日最高温(TMX)外, 其他气候带最主要的特征因子均为云覆盖百分数(CCP)。值得注意的是, 若根据分布面积考虑前3类最重要的因子, 可以发现以温带森林草原气候与温带草原气候为分界线, 该界线以北的寒冷干旱区域ETa主要受云覆盖百分数(CCP)、月霜日频率(FDF)与月均水汽压(VAP)共同驱动, 而该界线以南温暖湿润区域的ETa则主要受云覆盖百分数(CCP)、月均气温日较差(DTR)与月均日最高温(TMX)共同驱动。

对于水文区划, 平水、少水及干涸区域来说最主要的影响因子为云覆盖百分数(CCP), 秦巴、大别山北亚热带多水区与西南亚热带、热带多水区最主要的影响因子为月均气温日较差(DTR), 东南亚热带、热带丰水地区主要以云覆盖百分数(CCP)为主, 滇西、藏东南亚热带、热带丰水区则主要以月均日最高温(TMX)为主。在丰水、多水区与少水、平水区的分界线(即Ⅹ与Ⅴ的交界线), 亦即以横断山脉为分界, 可以看到横断山脉以南的丰水区的ETa主要受云覆盖百分数(CCP)、月均气温日较差(DTR)与月均日最高温(TMX)共同驱动, 而横断山脉以北的少水区域主要受云覆盖百分数(CCP)、月霜日频率(FDF)与月均水汽压(VAP)共同驱动, 该结论与按照气候带分区探讨基本一致。不可忽视的是, 无论是在丰水区还是少水区, 云覆盖百分数(CCP)都是所有因素中最主要的驱动因子。

3 结论与讨论 3.1 讨论

本文基于GLDAS 2.1及CRU数据集, 对中国西南地区1966—2016年月尺度的实际蒸散发数据进行了逐像元反演计算, 结合MSEOOB、PVE及RMSE评价方法以及与其他典型数据集对比的方法对模型和反演结果进行了充分的精度评价。在此基础上, 对中国西南地区ETa的空间格局及时空演变特征进行了分析。实际蒸散发对生态安全格局的形成及演变过程有着极其重要的影响, 因此, 在长时间尺度上对ETa数据的空间格局及其时空演变特征有充分的了解对充分明晰生态安全格局的形成机制及演变机理有着极其重要的参考价值。

通过分析各个因子对于实际蒸散发数据计算的重要程度, 可以在空间上辨析影响ETa计算的特征因子, 这些因子并没有包含所有影响ETa的因素, 如NDVI等, 但研究已经在长时间尺度限制的基础上充分考虑了所有可能影响ETa的因素。在不同的地理、地质条件下, 影响ETa的因素不同, 如云覆盖百分数在研究区内各个区域都有显著的占比情况, 这是因为云覆盖百分数直接影响到地表及植被冠层所接收到的太阳辐射, 从而影响ETa的量级;此外, 对于西藏、青海及川渝地区, 特别是川藏交界等干旱区域, 月均水汽压较为集中地分布, 其原因在于, 对于干旱区域, 月均水汽压亏缺和土壤水分一样重要, 它直接决定了蒸散发的量级。明晰中国西南不同气候类型、水文区划下的影响因子, 能够系统地了解区域安全格局形成机制以及其演变过程机理。此外, 对于开展生态区划、生态功能区划以及生态红线划定等相关工作也提供了指标和参数选择的科学支撑。因此, 辨析最关键的影响因素对全面了解ETa的演变和区域生态安全格局的形成和演变机理是至关重要的。

3.2 结论

(1) 对中国西南地区近50年陆面月尺度ETa数据进行逐像元反演, 构建了1966年至2016年各月尺度的实际蒸散发数据集, 模型整体的MSEOOB均值为4.14, 标准偏差仅为3.73, 模型整体PVE参数均值为99.36%, 标准偏差仅为0.33, 充分说明了模型整体具有较高的稳定性、可靠性和较高的泛化能力。此外模型对于2000年至2016年月尺度拟合结果的RMSE指数均值仅为1.04 mm/月, 标准偏差为0.52, 总体而言模型的拟合精度较高。

(2) 通过与GLDAS 2.1、2.0版本ETa数据集及MOD16 ETa数据对比, 拟合结果与GLDAS 2.1数据集的R2系数达到了0.99, 与MOD16数据集的R2系数达到了0.95, 与GLDAS 2.0数据集的R2系数超过了0.89, 充分说明了拟合结果在研究期间的可信度和精度是较高的。

(3) 研究区ETa整体上表现出随着纬度的降低而增加的特征, 不同季节上中国西南的ETa空间分布差异较为明显。在空间上, ETa从西北高原地区向东南沿海区域逐步增加。ETa随着季节的变化从春季到夏季先呈现出由东南向西北逐步增加的态势, 夏季到冬季则呈现出从西北向东南减弱的特征, 然而, 在全年尺度上, 西藏和青海西北部, ETa始终保持较低的值域范围。从1966年至2016年, ETa最大值呈现出逐步增加的趋势, 而最小值则呈现出减小的趋势。ETa在每年的7、8月份左右达到最大值, 在1、2月份左右为最低值, 并呈现起伏的周期特征。

(4) 研究区内, 以横断山脉为界, 横断山脉以南的丰水区的ETa主要受云覆盖百分数(CCP)、月均气温日较差(DTR)与月均日最高温(TMX)共同驱动, 而横断山脉以北的少水区域主要受云覆盖百分数(CCP)、月霜日频率(FDF)与月均水汽压(VAP)共同驱动, 而无论是在丰水区还是少水区, 云覆盖百分数(CCP)都是所有因素中最主要的驱动因子。

参考文献
[1] 刘国华. 西南生态安全格局形成机制及演变机理. 生态学报, 2016, 36(22): 7088–7091.
[2] Xiong Y J, Zhao S H, Tian F, Qiu G Y. An evapotranspiration product for arid regions based on the three-temperature model and thermal remote sensing. Journal of Hydrology, 2015, 530: 392–404.
[3] 刘世梁, 田韫钰, 尹艺洁, 安南南, 董世魁. 云南省植被NDVI时间变化特征及其对干旱的响应. 生态学报, 2016, 36(15): 4699–4707.
[4] Li H W, Wang S J, Bai X Y, Luo W J, Tang H, Cao Y, Wu L H, Chen F, Li Q, Zeng C, Wang M M. Spatiotemporal distribution and national measurement of the global carbonate carbon sink. Science of the Total Environment, 2018, 643: 157–170. DOI:10.1016/j.scitotenv.2018.06.196
[5] Dolman A J, de Jeu R A M. Evaporation in focus. Nature Geoscience, 2010, 3(5): 296. DOI:10.1038/ngeo849
[6] Ghilain N, Arboleda A, Gellens-Meulenberghs F. Evapotranspiration modelling at large scale using near-real time MSG SEVIRI derived data. Hydrology and Earth System Sciences, 2011, 15(3): 771–786. DOI:10.5194/hess-15-771-2011
[7] 王静爱, 左伟. 中国地理图集. 北京: 中国地图出版社, 2010: 34–41.
[8] 郑朝菊, 曾源, 赵玉金, 高文文, 赵旦, 吴炳方. 20世纪90年代以来中国西南地区土地覆被变化. 生态学报, 2016, 36(23): 7858–7869.
[9] University of East Anglia Climatic Research Unit, Harris I C, Jones P D. CRU TS4.01: Climatic Research Unit (CRU) Time-Series (TS) version 4.01 of high-resolution gridded data of month-by-month variation in climate (Jan. 1901-Dec. 2016). Centre for Environmental Data Analysis. (2017-12-04). http://dx.doi.org/10.5285/58a8802721c94c66ae45c3baa4d814d0.
[10] Harris I, Jones P D, Osborn T J, Lister D H. Updated high-resolution grids of monthly climatic observations-the CRU TS3.10 dataset. International Journal of Climatology, 2014, 34(3): 623–642. DOI:10.1002/joc.3711
[11] Briffa K R, Shishov V V, Melvin T M, Vaganov E A, Grudd H, Hantemirov R M, Eronen M, Naurzbaev M M. Trends in recent temperature and radial tree growth spanning 2000 years across Northwest Eurasia. Philosophical Transactions of the Royal Society B:Biological Sciences, 2008, 363(1501): 2271–2284.
[12] 张东, 宋献方, 张应华, 杨丽虎, 杨胜天. 基于CRU格点数据集的近百年渭河流域降水变化. 干旱区资源与环境, 2018, 32(2): 142–148.
[13] Rodell M, Beaudoing H K, NASA/GSFC/HSL. GLDAS Noah land surface model L4 monthly 0.25 x 0.25 degree V2.1, . Greenbelt, Maryland, USA: Goddard Earth Sciences Data and Information Services Center (GES DISC), 2007.
[14] Houser P R. Land data assimilation systems//Swinbank R, Shutyaev V, Lahoz W A eds. Data Assimilation for the Earth System. Dordrecht: Springer, 2003: 345-360. http://dx.doi.org/10.1007/978-94-010-0029-1_30.
[15] Park J, Choi M. Estimation of evapotranspiration from ground-based meteorological data and global land data assimilation system (GLDAS). Stochastic Environmental Research and Risk Assessment, 2015, 29(8): 1963–1992. DOI:10.1007/s00477-014-1004-2
[16] Lv M X, Ma Z G, Yuan X, Lv M Z, Li M X, Zheng Z Y. Water budget closure based on GRACE measurements and reconstructed evapotranspiration using GLDAS and water use data for two large densely-populated mid-latitude basins. Journal of Hydrology, 2017, 547: 585–599.
[17] Breiman L. Random forests. Machine Learning, 2001, 45(1): 5–32.
[18] Efron B. Bootstrap methods:another look at the jackknife. The Annals of Statistics, 1979, 7(1): 1–26.
[19] Verikas A, Gelzinis A, Bacauskiene M. Mining data with random forests:a survey and results of new tests. Pattern Recognition, 2011, 44(2): 330–349. DOI:10.1016/j.patcog.2010.08.011
[20] Iverson L R, Prasad A M, Matthews S, Peters M. Estimating potential habitat for 134 eastern US tree species under six climate scenarios. Forest Ecology and Management, 2008, 254(3): 390–406.
[21] Prasad A M, Iverson L R, Liaw A. Newer classification and regression tree techniques:bagging and random forests for ecological prediction. Ecosystems, 2006, 9(2): 181–199.
[22] Fernández-Delgado M, Cernadas E, Barro S, Amorim D. Do we need hundreds of classifiers to solve real world classification problems?. Journal of Machine Learning Research, 2014, 15(1): 3133–3181.
[23] Adelabu S, Mutanga O, Adam E. Testing the reliability and stability of the internal accuracy assessment of random forest for classifying tree defoliation levels using different validation methods. Geocarto International, 2015, 30(7): 810–821. DOI:10.1080/10106049.2014.997303
[24] Kamusoko C, Gamba J. Simulating urban growth using a random forest-cellular automata (RF-CA) model. ISPRS International Journal of Geo-Information, 2015, 4(2): 447–470.
[25] 姜艳阳, 王文, 周正昊. MODIS MOD16蒸散发产品在中国流域的质量评估. 自然资源学报, 2017, 32(3): 517–528.
[26] Syed T H, Famiglietti J S, Rodell M, Chen J L, Wilson C R. Analysis of terrestrial water storage changes from GRACE and GLDAS. Water Resources Research, 2008, 44(2): W02433.
[27] Vicente-Serrano S M, van der Schrier G, Beguería S, Azorin-Molina C, Lopez-Moreno J I. Contribution of precipitation and reference evapotranspiration to drought indices under different climates. Journal of Hydrology, 2015, 526: 42–54.
[28] Stoy P C, Katul G G, Siqueira M B S, Juang J Y, Novick K A, McCarthy H R, Oishi A C, Uebelherr J M, Kim H S, Oren R. Separating the effects of climate and vegetation on evapotranspiration along a successional chronosequence in the southeastern US. Global Change Biology, 2006, 12(11): 2115–2135. DOI:10.1111/gcb.2006.12.issue-11
[29] 邓兴耀, 刘洋, 刘志辉, 姚俊强. 中国西北干旱区蒸散发时空动态特征. 生态学报, 2017, 37(9): 2994–3008.
[30] Douville H, Ribes A, Decharme B, Alkama R, Sheffield J. Anthropogenic influence on multidecadal changes in reconstructed global evapotranspiration. Nature Climate Change, 2013, 3(1): 59–62.
[31] New M, Hulme M, Jones P. Representing twentieth-century space-time climate variability. Part Ⅱ:development of 1901-96 monthly grids of terrestrial surface climate. Journal of Climate, 2000, 13(13): 2217–2238.
[32] Qian Y, Wang W G, Leung L R, Kaiser D P. Variability of solar radiation under cloud-free skies in China:the role of aerosols. Geophysical Research Letters, 2007, 34(12): L12804.
[33] 王任超, 张利平, 徐霞. 南水北调中线工程水源区蒸散发计算方法比较及影响因素分析. 长江流域资源与环境, 2012, 21(S1): 127–133.
[34] 董煜.艾比湖流域气候与土地利用覆被变化的径流响应研究[D].乌鲁木齐: 新疆大学, 2016. http://cdmd.cnki.com.cn/Article/CDMD-10755-1016767956.htm
[35] 李修仓.中国典型流域实际蒸散发的时空变异研究[D].南京: 南京信息工程大学, 2013. http://cdmd.cnki.com.cn/Article/CDMD-10300-1013340778.htm
[36] Yang B, Kang S Y, Ljungqvist F C, He M H, Zhao Y, Qin C. Drought variability at the northern fringe of the Asian summer monsoon region over the past millennia. Climate Dynamics, 2014, 43(3/4): 845–859.
[37] Ray J D, Gesch R W, Sinclair T R, Allen L H. The effect of vapor pressure deficit on maize transpiration response to a drying soil. Plant and Soil, 2002, 239(1): 113–121. DOI:10.1023/A:1014947422468