基于多源数据与多模型集成的流域人为蒸散发变异评估

韦 林，段 凯，2，刘效东，林玉茹，陈晓宏，王小辣

（1.中山大学 土木工程学院，广东 广州 510275；2.南方海洋科学与工程广东省实验室（珠海），广东 珠海 519082；3.华南农业大学 林学与风景园林学院，广东 广州 510642；4.长江水利委员会 长江科学院，湖北 武汉 430010）

1 研究背景

蒸散发包括地表植被蒸腾、冠层截留、水面和土壤的蒸发，是水文循环和能量平衡的关键要素［1］。蒸散发不仅受气候变化的影响，还受水库调蓄、消耗性用水、跨流域调水、土地利用/覆被变化等人类活动的影响［2-4］。在人类活动日益加剧的背景下，评估流域或区域尺度上的人为蒸散发变异对水资源可持续管理至关重要［5-6］。

然而，蒸散发的实测数据稀少且难以获取，流域蒸散发的估算通常基于对站点/小流域尺度上观测数据的升尺度分析、流域水文模型模拟、或大尺度遥感数据产品的应用，面临着空间异质性与尺度差异的难题，仍存在较大不确定性［7-8］。目前研究人为蒸散发变异的思路主要有两种，一种是在传统水文模型中加入取水、耗水、灌溉等模块，用于模拟人类活动对蒸散发及其它主要水文过程的影响［9-12］。由于缺少人类活动的相关数据以及对其影响机制的认知有限，模型相关参数的区域性显著，一般很难在全球陆面水文模型中实现参数化［13］，因此大多数陆面水文模型侧重于模拟自然条件下的蒸散发。第二种思路是比较模型模拟蒸散发与水量平衡法得到的蒸散发之间的差异，从而衡量人类活动的影响程度［14-15］。重力恢复与气候实验（GRACE）卫星观测的陆地水储量变化被广泛用于水量平衡计算［16-17］。在较大的空间尺度上，陆面水文模型模拟的天然蒸散发与基于水量平衡法估算的实际蒸散发之间的差异可被认为是人类活动所引起的蒸散发变化［18-19］。例如，Castle 等［14］利用水量平衡与北美陆面数据同化系统（NLDAS）数据，识别了美国科罗拉多河流域生长季节期间人类活动引起的蒸散发变化；Pan 等［15］比较了水量平衡法与全球陆面数据同化系统（GLDAS）模拟的蒸散发，发现人类活动导致海河流域年蒸散发增加了12%；Chen 等［20］用类似的方法，发现人类活动导致松花江流域年蒸散发减少了12%；Liu 等［21］对比水量平衡法与VIP-RS 模拟结果，发现近年来人类活动引起子牙-大清河流域蒸散发量持续增长，增幅为4.88 km3/a。该类研究的评估结果高度依赖于模型的选择，尽管有研究利用人类活动影响之前的时期作为参照，校正陆面模式的天然径流数据［22］，但模型结构的不确定性仍是限制蒸散发评估精度的重要因素［23］。

多模型集成技术是弥补单个数据产品或单个模型系统偏差的有效工具，可从单值预测［24-25］或概率密度预测［26-27］的角度集成不同模型的优势。其中，以贝叶斯模型加权平均（Bayesian model averaging，BMA）为代表的多模型概率集成方法可综合多个模型模拟值的后验分布，提供更为可靠的概率分布预测，已在气象与水文预报等领域得到广泛应用，但在蒸散发研究中尚不多见［28-30］。为充分利用现有多源数据［31］，合理评估人类活动导致的流域蒸散发变异并量化其不确定性，本文构建了一种基于多源数据与BMA 多模型集成的人为蒸散发变异评估框架，通过综合利用地面观测、水资源统计、卫星遥感、与陆面模型等数据，并以概率方式集成多个模型的模拟结果，以提高评估结果的可靠性，为强烈人类活动影响下的流域水资源管理提供科学依据。

2 研究方法

本研究技术路线如图1 所示。通过BMA 集成多种陆面水文模型对天然条件下流域蒸散发的模拟结果（即“天然蒸散发”），利用流域实测降水、径流序列与GRACE 重力卫星提供的陆地水储量数据推算实际发生的流域蒸散发（即“实际蒸散发”），进而使用“天然蒸散发”与“实际蒸散发”的差值，评估人类活动所导致的流域蒸散发变异及其不确定性区间。

图1 人为蒸散发变异评估框架

2.1 数据预处理搜集整理研究区域的降水、径流、潜热通量等地面观测数据，相关地区的水资源统计数据，GRACE 卫星陆地水储量数据，陆面水文模型的降水、蒸散发、径流再分析数据产品等基础数据。结合地面实测数据，对统计数据、遥感数据与再分析数据进行尺度同化与校正处理。统计分析流域内的水资源供用耗排、水库调蓄与渗漏损失、跨流域调水等因素，应用逐项还原法，将实际径流还原为天然径流。

2.2 天然蒸散发计算

2.2.1 陆面水文模型优选与集成 由于可用于模型校验的蒸散发数据较为缺乏，本文根据备选模型对于天然径流的模拟精度进行模型优选，在结合降水输入数据进行径流校正处理的前提下，假设模型的径流模拟精度可有效反映流域蒸散发的模拟效果。模型优选按以下步骤进行：首先，分析单个模型的模拟效果，确定备选模型；其次，使用BMA 方法进行多模型集成，分析不同模型组合在集成后的模拟精度与不确定性；最后，综合分析单个模型与模型集成在研究区域的适用性。

BMA 是一种基于贝叶斯理论的统计后处理方法，可以结合多个模型进行联合推断和预测，其核心是以后验概率作为权重，对备选模型进行概率的加权平均［32-33］：

式中：y为多模型集成结果；f=f1，…，fk为k个模型的模拟值；yT为通过实测值；p(fk|yT)为第k个模型模拟值的后验概率，即模型权重；pk(y|(fk，yT))表示与单个模型的模拟结果fk相联系的条件概率密度函数。BMA 多模型集成的期望值和方差分别为：

式中：g(y|fk，σk2)表示均值为fk、方差为σk2的正态分布；wk为第k个模型的权重；方差包含两项，分别为模式间方差与模式内方差。模型集成中单个模型的权重wk和方差σk2可通过期望最大化（EM）算法求解，算法流程参见文献［26，34］。

2.2.2 模型评价指标 本文采用纳什效率系数（NSE）［36］、相关系数（r）、均方根误差（RMSE）［37］作为模拟精度的评价指标，使用覆盖率（CR）、平均相对带宽（RB）来衡量BMA 集成模拟的不确定性［33］，通过综合考虑模拟精度与不确定性范围来进行模型优选。表1所示为各个指标的计算公式。

表1 模型评价指标及计算公式

2.2.3 天然蒸散发数据的概率集成 通过将BMA 权重赋予备选模型模拟的蒸散发序列，最终得到蒸散发集成模拟结果与不确定性区间。步骤如下：（1）考虑到天然情况下的径流和蒸散发变量通常是非正态的，而BMA 方法假定条件分布需服从正态分布，首先使用Box-Cox 变换对各个模型的蒸散发序列进行正态转换；（2）使用蒙特卡罗组合抽样方法随机抽选出某个备选模型，重复抽选100 次以生成逐时段的样本分布；（3）将样本从小到大排序，使用5%和95%分位数分别作为上下边界，即可得到90%置信区间内的不确定性范围；（4）使用Box-Cox 逆变换，得到多模型集成蒸散发序列的期望值与不确定性范围。具体算法可参见文献［26，34-35］。

2.3 实际蒸散发与人类活动影响评估根据流域水量平衡原理，利用降水、实测径流、陆地水储量变化［16］计算逐月的流域实际蒸散发：

式中：P为降水量；Robs为实测径流量；TWSC为时段内的陆地水储量变化量。

TWSC计算公式为：

式中TWSAt为第t个月的陆地水储量距平值。

实际蒸散发与天然蒸散发的差值可视为人类活动引起的流域蒸散发变化（ETH）：

式中：ETobs为基于水量平衡法估算的实际蒸散发；ETnat为模型模拟的天然蒸散发。

3 案例研究

3.1 研究区域与数据本文以珠江流域为例。珠江流域（东经102°14′—115°53′，北纬21°31—26°49′）地处华南热带亚热带气候地区（图2），研究区域面积为44.21万km2，包括西江、北江和东江三大子流域。当前珠江流域的流域尺度蒸散发研究主要集中在时空变化规律分析［38-40］，针对人为蒸散发变异的系统研究较为缺乏。所使用的数据如表2所示。

图2 珠江流域水系、站点与土地利用类型分布

表2 研究数据

3.1.1 地面观测数据 降水数据来自中国区域地面气象要素数据集（CMFD）［41-42］，该数据集通过地面气象站点的实测数据插值生成，时间分辨率为3 h，水平空间分辨率为0.1°×0.1°。水文数据包括高要、石角、博罗三个主要控制站的逐月实测径流数据。此外，潜热通量观测数据可通过汽化系数（λ=2.45×106J/kg）换算为点尺度上的实测蒸散发，本文使用ChinaFLUX 鼎湖山站通量系统的数据代表珠江流域森林生态系统蒸散发［43］，作为模型模拟精度的参照标准之一。

3.1.2 GRACE 卫星数据 GRACE 卫星通过监测时变重力场变化得到地表质量变化，并转化为厘米级的等效水高变化［44］，数据序列最早始于2002年4月。本文使用CSR 和JPL 2 家研究机构公布的2003—2016年RL06 Mascon 数据产品，为了减少由不同GRACE 后处理策略带来的不确定性，采用这2 种数据计算结果的算数平均值作为研究区域的陆地水储量距平［45］。由于GRACE 卫星在其运行过程中因自身轨道调整等因素导致部分数据缺失，采用K 近邻算法（K Neareast Neighbor，KNN）插补缺失值［46］。

3.1.3 陆面水文模型数据 本文采用ERA-Interim［47］、GLDAS2［48］的天然蒸散发、天然径流、降水数据。ERA-Interim 数据为逐日时间步长，空间分辨率为0.25°。GLDAS2 采用目前发布的最新版本GLDAS2.1，包含Noah、VIC、CLSM 模式，三个模式由于缺少人类活动的参数仅模拟自然条件下的蒸散发［14］。其中，Noah 数据空间分辨率为0.25°，VIC 和CLSM 数据空间分辨率为1.0°，为保持一致性，用双线性插值法将VIC 和CLSM 数据从1.0°插值到0.25°。

3.2 水文模型适用性分析

3.2.1 降水 图3 为2003—2016年4 个陆面水文模型模拟的降水、径流与实测降水、Rnat的年内变化特征。受东南季风和西南季风的影响，珠江流域降水具有明显的季节性，4 种模拟数据均呈现出夏季高、冬季低的特点，其中，CLSM、Noah 的降水强迫数据是相同的，CLSM、VIC 的季节变化特征与CMFD 降水最为接近。在珠江流域，4 个模型模拟的降水与CMFD 降水的r均在0.9 以上，Noah 的精度最高，与CMFD 的r为0.98，ERA-Interim 与CMFD 的r最低，为0.91。从西江、北江、东江三个子流域来看，Noah、CLSM、VIC的RMSE（17.63 ～31.30 mm/月）小于ERA-Interim（47.05 ～65.37 mm/月）。在珠江流域多年平均降水量上，CMFD 为1418.38 mm/a，Noah 为1439.82 mm/a，CLSM 和VIC 为1440.52 mm/a，ERA-Interim 为1707.74 mm/a，ERA-Interim 高估了实际降水量。降水作为水文模型最重要的输入，其准确性对模型径流和蒸散发的模拟至关重要。为减小误差，在分析前使用CMFD 降水与模型降水数据的比值校正模型模拟的径流、蒸散发数据。

图3 四种陆面水文模型模拟的多年平均月降水与月径流

3.2.2 径流 在基于实测降水进行校正后，Noah、CLSM、VIC、ERA-Interim 模型模拟径流（图3）与还原计算后的天然径流的r分别由0.92、0.92、0.94 和0.77 提高至0.94、0.93、0.95 和0.86；RMSE则由31.64、17.12、29.50 和32.97 mm/月减小至27.31、15.58、27.93 和23.13 mm/月，说明校正后的模型模拟精度得到有效提高。校正后的Noah、CLSM、VIC、ERA-Interim模型模拟的珠江流域多年平均天然径流量分别为656.36、464.92、963.60 和486.73 mm/a，相较于实测结果（678.84 mm/a），Noah、CLSM、ERA-Interim 模型在不同程度上低估了天然径流量，VIC 模型则呈现高估，其中Noah 模型对于均值的模拟精度最高，相对误差为3%。各模型模拟的月径流在年内分布特征上均呈先增大后减小的单峰型分布，与实测数据（Rnat）相符。在西江、东江流域，Noah 与Rnat最为接近，而VIC 明显高估了Rnat。在径流最大的6月，ERA-Interim 在北江、东江流域存在明显低估。

3.2.3 蒸散发 图4 为2003—2009年模型模拟蒸散发与鼎湖山（DHF）站点实测数据的比较，其中，横轴和纵轴黑线为标准差，纵轴蓝色虚线为r，绿色虚线圆弧为RMSE。可以看出，4 种模型均与实测数据显著相关，r均达到0.6 以上。其中，VIC 模型总体精度最高，r为0.78，RMSE为15.24 mm/月，标准差与实测数据相近，为20.99 mm/月；Noah 模型精度次之，r为0.77，RMSE为17.74 mm/月，标准差为27.17 mm/月；ERA-Interim与CLSM误差略大，r分别为0.69和0.63，RMSE分别为34.17和28.04 mm/月，标准差均处于30 至50 mm/月之间。较高的相关性与相近的标准差表明，4 个陆面水文模型均能较好地捕捉到蒸散发的季节变化特征。

图4 四种陆面水文模型在鼎湖山站的蒸散发模拟精度

3.3 模型优选与概率集成图5 为4 种单一模型与11 种不同模型集成的径流模拟精度及不确定性。横坐标中的下标“N”“C”“V”“I”分别表示Noah、CLSM、VIC、ERA-Interim 等4 种模型，例如，BMANC表示利用BMA 方法集成Noah 与CLSM 模型，以此类推。从图5（a）可以看出，单一模型的性能无法得到保障，Noah 模型在东江流域、VIC 模型在西江流域均出现NSE为负值的情况，而BMA 方法可集成不同模型的优势，避免单一模型中的不确定性所带来的误差，由此提高模拟精度。例如，在西江流域，单一模型NSE最大为0.71（CLSM），集成Noah、CLSM、ERA-Interim 模型后NSE达到0.80（BMANCI），提高13%；在东江流域，单一模型NSE最大仅为0.44（CLSM），集成CLSM 与VIC、ERA-Interim 模型后NSE达到0.72（BMACVI），显著提高64%。

所有集成结果均具有较高的覆盖率（图5（c）），达到85%以上，但平均相对带宽（RB）的结果表明，不同模型集成产生的不确定性范围存在较大差异（图5（b）），因此，选择集成模型时需要权衡模拟精度与不确定性范围。在西江流域，平均模拟精度最高的模型集成为BMANCI（NSE=0.80）与BMACVI（NSE=0.78），但BMANCI的平均相对带宽为6.23，远大于BMACVI（RB=1.95），故在西江流域选择BMACVI更为合适。其中，CLSM、VIC、ERA-Interim 模型的权重分别为0.25、0.48 和0.27。在北江流域与东江流域，则分别优选出NSE最大而RB较小的BMACV与BMACVI模型集成。在北江流域，CLSM、VIC 模型的权重分别为0.26 和0.76；在东江流域，CLSM、VIC、ERA-Interim 模型的权重分别为0.26、0.46 和0.28。

图5 陆面水文模型的优选指标结果

3.4 人为蒸散发变异评估图6 所示为珠江流域实际蒸散发（ETobs）、天然蒸散发（ETnat）以及人为蒸散发变异（ETH）的逐月变化过程。利用优选出的模型集成与相应的权重，可得到珠江流域ETnat的多年平均值为677.57 mm/a。基于水量平衡原理得到的珠江流域2003—2016年ETobs的多年平均值为814.39 mm/a，与Zhong 等［49］的研究结果相似。受水汽条件与温度影响，珠江流域ETnat的月变化呈单峰型分布，主要集中在5—9月，占ETnat总量的60%，最大值出现在7月，为99.07 mm/月，最小值出现在2月，为24.97 mm/月。由于受到非平稳的人类活动影响，ETobs显示出比ETnat更大的年内振荡，尤其是在人类活动较为强烈的东江流域。西江、东江流域ETobs的月变化呈现双峰型分布，在北江流域为三峰型分布，峰值首次出现时间均比ETnat峰值出现时间早1～2 个月。

图6 珠江流域2003—2016年实际蒸散发（ETobs）、天然蒸散发（ETnat）与人为蒸散发变异（ETH）

就珠江流域全流域而言（图7（a）），2003—2016年人类活动导致蒸散发平均增大136.82 mm/a（21%），在西江、北江与东江子流域，人类活动分别导致蒸散发增大125.48 mm（19%）、232.19 mm/a（37%）、149.80 mm/a（18%）。珠江流域ETH整体呈增加趋势（图7（b）），平均增加速率为7.02mm/a，但趋势并不显著（R2=0.18，p＞0.05）。从年内分布来看，珠江流域ETH最大值出现在5月，为37.52 mm/月，最小值出现在7月，为-26.69 mm/月。西江和东江ETH最大值均出现在5月，分别为37.14和53.01 mm/月，北江流域ETH最大值出现在1月，为46.57 mm/月；而三个子流域的ETH最小值均出现在7月，分别为-27.00、-30.66、-16.44 mm/月。ETH最大值出现在5月可能与珠江流域以水稻种植为主的农业结构有关［50］。在5月份，水稻初插秧灌溉量增大，导致蒸散发显著增大。而在7月份，人类活动导致蒸散发显著减小，这可能与该月水库开闸泄洪、降水较6月和8月偏少（图3）以及数据误差等因素有关。值得注意的是，在经过模型优选与集成后，人为蒸散发变异的评估结果仍存在较大的不确定性。在90%置信区间内，珠江流域多年平均逐月ETH的不确定性范围在38.09～71.48 mm/月之间波动。在三个子流域之中，东江流域水资源开发利用程度最高，人类活动对水循环的干扰最为显著［51］，导致了该流域ETH存在较大的不确定性（53.73～133.31 mm/月），尤其是在4—7月。

图7 珠江流域水文变量多年平均值与蒸散发年际变化趋势

4 结论

本文构建了一种基于多源数据与BMA 多模型集成的人为蒸散发变异评估框架，并将其应用于珠江流域，主要结论如下。

（1）Noah、CLSM、VIC、ERA-Interim 等4 种陆面水文模型的模拟结果与实测结果在季节变化特征上较为一致，但仍存在一定的系统性误差，在结合地面实测降水数据进行校正后，模拟效果得到有效提升，其中Noah 模型与VIC 模型分别在珠江流域径流与天然蒸散发的模拟上体现出更强的适用性。

（2）4 种单一模型与11 种不同模型集成的径流模拟结果表明，BMA 方法可有效集成不同模型的优势，在一定程度上提高模拟精度。但是，由于大尺度陆面模型在应用于区域尺度时往往有较大的不确定性，优选集成模型时应权衡考虑平均模拟值的精度（如NSE）与不确定性区间（如RB）。在本例中，BMACVI、BMACV与BMACVI模型集成分别在西江、北江与东江流域体现出最佳的适用性。

（3）人类活动导致珠江流域蒸散发增大136.82 mm/a（21%），其中，西江、北江、东江流域蒸散发分别增大125.48 mm/a（19%）、232.19mm/a（37%）和149.80 mm/a（18%）。

（4）珠江流域人为蒸散发变异的年内最大值与最小值分别出现在5月（+37.52 mm/月）和7月（-26.69 mm/月），月际振荡特征可能主要受灌溉耗水、水库调蓄等因素控制。

本文侧重于在子流域尺度上解析人类活动对天然蒸散发过程与流域水量平衡的干扰及其不确定性，相较于MODIS-ET 等遥感ET 产品，GLDAS 等陆面水文模型能够提供更为完整的流域水循环过程数据，便于结合降水与径流的实测数据校验蒸散发，然而对于人类活动影响的识别仍受限于实测站点的数量。此外，本文中的分析结果主要反映了土地利用类型改变（如城市化、造林）［52］、水库调蓄与蒸发渗漏损失［53］、消耗性用水等人类活动对流域蒸散发的综合影响及其不确定性，将陆面水文模型所模拟的蒸散发视为气候演变作用下的天然蒸散发，并未考虑人类活动对陆面模型大气驱动要素及蒸散发的间接影响（如农业灌溉通过影响地表气温和湿度而干扰蒸散发过程［54］），未来需要针对人为影响空间异质性的精细化表征、陆-气双向反馈关系等问题展开进一步研究。