葛 路,詹良通,江衍铭
(1.浙江大学建筑工程学院水文与水资源工程研究所,杭州310030;2.浙江大学建筑工程学院岩土工程研究所,杭州310030)
IPCC第五次气候变化评估报告指出,全球气候系统的变暖是无可置疑的[1],并且我国科技部《第三次气候变化国家评估报告》中显示,我国气候变暖速率高于全球平均值[2]。气候变暖对于水文过程的影响也是显著的[3],国内外学者积极投入到气候变化对水文过程的影响研究中:周一飞等[4]研究了未来青海湖流域蒸散量、径流深和深层渗透量的变化趋势,曹丽娟等[5]预估气候变化对气温、降水、蒸发和径流的影响,Mantua N 等[6]分析夏季溪流温度和季节性流量以便评估太平洋鲑鱼淡水栖息地对气候变化的敏感性,王国庆[7]等研究发现我国淮河及其以南江河径流变化的主要原因是气候变化。
全球气候模式(Global Climate Model,GCM)是基于旋转球体的纳维-斯托克斯方程(Navier-Stokes equations)来描述地球大气、海洋和陆地物理过程的数学模型,被广泛应用于气候变化、气象预报等研究[8]。然而GCM 空间分辨率大,无法模拟区域气候特征且难以直接作为水文模型的输入,因此要对GCM进行降尺度处理。降尺度方法分为动力降尺度和统计降尺度,张徐杰等[9]比较了动力降尺度和统计降尺度方法对钱塘江流域设计暴雨的影响;Trigo 等[10]利用ANNs 对HadCM2 的降水数据进行降尺度。动力降尺度优点是物理意义明确,缺点是计算量大、费时,且拟合效果不佳;统计降尺度(Statistic DownScaling Model,SDSM)将GCM 输出中物理意义较好、模拟较准确的气候信息应用于统计模式来减少GCM的系统误差,且不用考虑边界条件对预测结果的影响,简单易行、效果好[11]。然而研究发现,SDSM 存在必须基于未来参数不变的假设、预报变量对预报因子过于敏感等问题,而BCSD 方法既保留了观测场的时间统计特征,又保留了观测场的小尺度空间相关结构[12],不依赖其他气象因子,被国内外广泛应用:王林等[13]得到M-BCSD 方法能够很好地降低GCMs 在中国区域降水模拟误差,且该方法不具有模式依赖性;Shrestha R R 等[14]比较GCM-M-BCSD-VIC 模拟的降水、温度和径流量和加拿大区域气候模式(CRCM)模拟结果与实测值的吻合程度,结果发现GCM-M-BCSD-VIC 吻合更好。因此本文选择M-BCSD(Bias Correction Spatial Disaggregation)方法对GCM气象数据进行降尺度和校正。
李昱等[15]利用SWAT 模型模拟澜沧江-湄公河流域历史和未来水文过程;冉思红等[16]利用多模式输出气象数据驱动SWAT-RSG 模型来研究天山地区的冰川覆盖率河流的径流对气候变化的影响;邓鹏等[17]构建应用CMIP5 多模式气候数据结果驱动VIC 水文模型评估未来潘阳湖流域径流量变化;王尹萍等[18]采用CMADS 驱动SWAT 模型模拟张家山水文站月径流。可以看出,前人在研究气候变化下的流域径流影响中大都采用分布式水文模型,并取得了丰厚成果,然而分布式水文模型需要大量的输入数据,获取数据是一大难题,且时空分辨率转换而导致的模型不确定性问题也很突出。而近年来许多研究也证明了数据驱动型方法在水文领域的优势,其在数学方程映射关系未知的情况下,系统软件存储大量的映射关系预测模式,因此在复杂的非线性问题处理时具有明显的优势,且该方法凭借出色的模拟结果赢得很多研究者的喜爱:范宏翔[19]等利用LSTM 构建鄱阳湖气象-径流模型,划分的子流域训练期NSE 均高于0.94。其中BPNN 可挖掘降雨径流之间的非线性关系,具有较强的泛化能力,适用于提取径流的复杂成分,同时该模型的使用条件不受地区限制,被广泛应用于月径流模拟中:王栋[20]等基于BPNN 及其组合模型预测青海湟水河流域的月径流,R2达0.85 以上;UR Berndtsson[21]对人工神经网络和概念径流模型在两个挪威河流域的月径流模拟进行比较,R2分别为0.82 和0.71。元江流域中的元江县一带属于干热河谷,气候较特殊,因此本文采用BPNN 模型对气候变化下元江流域的月径流进行预测,解决分布式水文模型建模的难题,使空间资料稀缺的地区和气候条件特殊的流域也可做相关的流域水文响应研究,可以被其他研究者所借鉴。
本文的研究区域元江流域上游,属于唯一发源于云南境内的一条重要国际河流,流至越南,研究该流域具有重要指导意义;同时,历史上元江流域发生过重大洪水灾害,据资料统计,1986年10月上旬元江流域的大暴雨导致山洪暴发、山体滑坡以及泥石流等洪水灾害,受灾损失严重;另外,元江流域水资源丰富,随着云南省经济发展和人口增加对水资源压力变大。相比全国其他流域,目前针对元江流域气候变化下的水文响应研究很少,探究气候变化对元江流域的气象变化特征与水文过程的影响,不仅对流域水文理论、分析方法和水资源合理配置具有重要意义,且可为当地应对极端天气及灾害提供科学预警。本文选取M-BCSD 方法作为降尺度工具,基于CMIP6 模型计算蒸发和相对湿度和分析未来降雨、气温,并将此气象数据作为输入构建BPNN 模型,分析气候变化下未来元江流域径流演变规律。
本文采用M-BCSD 方法对GCM 模式降尺度处理,采用算术平均的集合方法生成GCMs 集合数据,预测元江流域未来降雨和气温,然后利用基于交叉验证下的BPNN 模型模拟元江流域月径流,并分析气候变化对元江流域水文过程的影响,如图1所示为研究技术路线图。
图1 技术路线图Fig.1 The flow chart of this study
政府间气候变化专门委员会(IPCC)第六次评估报告中提出了一套与温室气体排放相关的社会经济发展情景——共享社会经济路径(SSPs),第六次国际耦合模式比较计划(CMIP6)的气候变化情景将共享社会经济路径与典型浓度路径(RCPs)科学组合,提供更可靠的气候变化可能结果。本文挑选了SSP126、SSP245、SSP370、SSP585 四种路径,四种路径的CO2排放、CO2浓度、人为辐射强迫和地表温度是依次增加的,同时考虑了SSP1(Sustainability)可持续路径、SSP2(Middle of the Road)中间路径、SSP3(Regional Rivalry)区域竞争路径以及SSP5(Fossil-fueled Development)以化石燃料为主的发展路径。
本文结合亚洲季风区域的气候特点,参考其他学者对GCMs 适用性分析[22~24],选取5 个GCM 模式,信息如表1所示。根据世界气象组织(WMO)推荐:选取基准期为1971-2000年、未来时期为2041-2070年(未来第一时期)和2071-2100年(未来第二时期)。数据包括月尺度的降水、气温、相对湿度。数据覆盖经纬度范围为:99.5°E~102.5°E 和23°N~26°N。BPNN 模型输入所需的潜蒸发量,使用Hargreaves公式进行估算得到[25]:
表1 GCMs信息Tab.1 Information of the GCMs
式中:a为参数值取0.002 3;TD为日最大气温和最小气温的温差,℃;Ta为日平均温度,℃;Ra为用相对应的蒸发单位表达的地球辐射值,mm/d。
式中:Gsc为太阳常数为0.082 MJ/m2/min;dr为地球与太阳的相对距离,rad;ws为日落角度,rad;ψ为纬度,rad;δ为太阳高度角。
式中:J为儒略日;
式(2)需经单位转换代再入式(1)中,转换式为:
BCSD方法是Wood等于2004年提出的统计降尺度方法[26],可以解决GCMs 分辨率低和模拟偏差的问题,本文对5 个GCMs模拟的基准期和未来时期的降雨和温度数据按月进行偏差纠正,因此称为M-BCSD。该方法分为两个部分:偏差纠正和空间降尺度。具体做法如下:
(1)分位数映射法(Quantile Mapping,QM)[27]进行偏差纠正。QM 是基于频率分布的误差订正方法,利用非参数转换使GCMs 模拟变量的频率分布与实测变量一致。传统的QM 是使用累计分布函数(cdf),它的局限性是只考虑湿日的降雨概率,本文使用经验累计分布函数(ecdf)来同时考虑湿日和干日的降雨概率,经验累计分布函数包括零值,因此该方法不仅可以纠正降雨、气温数据的量值,也考虑了降雨事件发生与否。先用反距离权重插值法(IDW)将所有GCMs 的空间分辨率统一到1.5°×1.5°,对高分辨率的实测数据进行升尺度处理,即将流域内及周围10 个气象站点降水和气温数据也插值到与GCMs 相同位置;然后对模拟降水、气温和实测降水、气温的ecdf建立对应关系,引入偏差纠正因子,按月使模拟降水(气温)和实测降水(气温)的经验累计分布函数(ecdf)相同来订正模拟降水(气温)时序;最后利用ecdf的反函数结果可推出未来时期GCMs 的输出变量。QM方法的非参数转换函数如下式:
式中:X是GCMs 模拟的基准期的降雨或气温结果;Y是对应变量在基准期偏差纠正后的结果;Xfuture是GCMs 模拟的未来时期的降雨或气温结果;Yfuture是对应变量在未来时期偏差纠正后的结果;ecdfsim是变量X的经验累积分布函数;F是偏差纠正因子。
(2)空间统计降尺度处理。先将1.5°×1.5°的实测值和偏差纠正的GCMs 分别反距离权重插值到10 个站点,然后引入校正因子,根据插值前后的实测值和插值后的GCMs 来得到MBCSD方法降尺度后的GCMs模拟值。具体如下式:
式中:Si,j,t和分别为1.5°×1.5°的实测值和其经过IDW 插值到10 个站点的值;CFi,j,t为校正因子;为GCMs 集合数据经过IDW 插值到10 个站点的值;为经过M-BCSD 处理后的GCMs集合数据。
由于本文只有10年元江水文站流量数据,数据量少,容易造成BPNN 模型数据集划分不合理而导致训练过拟合或欠学习的问题。交叉验证是一种评估模型的泛化性能的方法,可以从有限的数据中获取尽可能多的有效信息,有效解决这类问题。交叉验证法分K 折交叉验证法(K-fold Cross-Validation)和留一法(Leave-One-Out Cross Validation),K 折交叉验证法是将原始数据均分为K份,固定其中1组为测试集,将其余每个子集数据分别作为一次验证集,其余K-2 组子集作为训练集,这样会得到K-1 个模型,将这K-1 个模型的误差结果平均值作为模型性能指标,留一法实际上就是将原始样本总数N作为K折交叉验证法中的K。由于留一法计算成本高,需要建立的模型数量与原始样本数量相同,因此本文使用K折交叉验证法。
本文采用四折交叉验证,将1991-2000年的气象水文数据按时序分为5 组,固定1993-1994年的数据作为测试集,其余4组分别作为一次验证集,剩下的数据作为训练集。这里M1、M2、M3 模型的参数在下文说明。用确定好的参数反复轮流迭代,取4 个试验的测试相对误差平均值和均方根误差平均值来评估模型的泛化性能。
倒传递神经网络(Back-Propagation Neural Network,BPNN)是由Rumelhart 和McClelland 为首的科学家在1986年提出的概念,是一种按照误差反向传播训练的多层前馈式神经网络,目前得到广泛的应用[28-30]。BP 算法包括信号的前向传播和误差的反向传播两个过程,误差反向传播会使权重不断更新。
本研究BPNN 模型利用MATLAB R2017b 版本中的feedforwardnet 函数构建,基本参数通过试错法得到[31],模型的输入变量组合通过相关分析得到,如表2所示,其中Q(t)表示t时刻站点的流量,P(t)表示t时刻站点的平均面雨量(流域内4 个气象站点数据通过泰森多边形法加权得到平均雨量),T(t)、Tmax(t)、Tmin(t)分别表示t时刻站点的平均、最高、最低温度,H(t)表示t时刻站点的相对湿度。
表2 各模型参数表Tab.2 Parameters of each model
采用实测气象数据、GCMs 集合数据、GCMs 集合与实测气象融合数据分别作为BPNN 模型的输入因子,建立以下3 种BPNN 模型,分别记为M1、M2、M3。其中,变量下标为S代表实测数据、变量下标为G代表GCMs模型集合结果。
选取的时段是1991-2000年。为了确保模型的有效性,本文对元江流域的数据进行了划分,训练、验证和测试阶段的数据比例为3∶1∶1,根据峰形选择1992、1994-1997年为训练数据,1991、1998年为验证数据,1993、2000年为测试数据。考虑到历史数据作为输入变量的变化范围可能未涵盖未来气候模式输出的变量的范围,因此将训练数据归一化为[0.2,0.8],保证模型的泛化能力和预测精度。
本研究采用均方根误差(RMSE)、纳什系数(Nash Sutcliffe Efficiency,NSE)和相对偏差(PBIAS)作为BPNN 模型评价指标,具体公式如下:
式中:Qobs(t)、Qobs(t)分别为t时刻的预测径流量和实测径流量;为率定期内实测径流量的平均值。RMSE值越小,说明模拟结果越好;NSE的范围从负无穷到1,值越接近1,说明模型效果越好;PBIAS最优值为0,绝对值15%~25%、10%~15%、小于10%分别表示模拟结果较好、好、非常好。
本研究采用相关系数(CC)、均方根误差(RMSE)、平均误差(ME)、和相对偏差(PBIAS)来评估降尺度处理后GCMs 与实测降雨相比的精度和可靠性[32],具体公式如下:
式中:Si为实测值;Gi为降尺度处理后的GCMs 模型值;-S和-G分别为平均实测值和GCMs 模型平均值。CC范围为0~1,值越接近1,说明处理效果越好;RMSE、ME和BIAS的绝对值越接近0,说明处理效果越好。
本文选取云南省元江流域的上游段为研究区域,北邻金沙江流域,西与澜沧江以无量山为分水岭,东接南盘江流域,南面与越南接壤。流域面积2.59 万km2,属中部高原温和区,海拔1 600~2 200 m,年平均气温15~18 ℃,常年不结冰,其中气温在3-10月较高,年平均降水量800~1 000 mm,降水一般集中在5-10月份。其中元江站属于干热河谷地区,海拔为400.9 m。图2为元江流域位置及其水文气象站点。本文数据主要包括从国家气象局获得流域内的10 个气象站的1971-2000年逐日降雨、气温(最高、最低、平均气温)和相对湿度数据,从云南省水文局获得元江水文站的1991-2000年逐日流量数据,逐月数据通过统计计算得到。调研可知[33],元江水文站建于1953年,1954年下迁15 m,迁站前后断面水量控制情况无系统性变化,资料系列具有一致性;1960-2012年元江流域内1999年降水量最多,且1976-1995年为平偏枯水期,1996-2002年平偏丰水期,选取的水文资料具有代表性。综上,1991-2000年流量数据可作为BPNN模型的历史模拟资料。
图2 元江流域位置及其水文气象站点Fig.2 Location of Yuanjiang River Basin and the hydrometeorological stations
通过M-BCSD 方法,得到1971-2000年各个GCM 降尺度处理后的结果,图3所示为元江流域(元江站)在基准期实测、BCC-CEM2-MR 模拟与偏差纠正后的降雨经验累积分布函数对比情况。可以看到,12个月除10月和11月,其余月份降雨模拟值都偏低,经过偏差纠正后,BCC-CSM2-MR 模型的经验累积分布函数跟实测降雨的经验累积分布函数几乎重合,而其余GCMs 模型的偏差纠正后结果与BCC-CEM2-MR 模型相似,便不一一赘述。通过偏差纠正后,接着以算术平均的集合方法生成GCMs 集合结果,比较基准期GCMs 与实测降雨和温度数据,结果如图4所示。可以看到,比较降雨数据,CC值都在0.7 以上,RMSE都小于50 mm,ME绝对值都小于2.5 mm,BIAS绝对值都小于3.5%;比较温度数据,CC值都接近1,RMSE都小于1.2 ℃,ME绝对值都小于0.4 ℃,BIAS绝对值都小于2%。从偏差纠正前后降雨经验累积分布函数图和评价指标表明:MBCSD 方法在元江流域的降水和温度上的效果都比较明显,其中对GCMs 温度数据处理的相对更好,验证了M-BCSD 方法在元江流域的适用性,采用此方法可对元江流域2041-2100年降雨和气温序列进行偏差纠正和降尺度。
图3 QM方法对12个月降雨的偏差纠正结果(BCC-CEM2-MR)Fig.3 Results of bias correction for rainfall in 12 months by QM(BCC-CEM2-MR)
图4 降雨、温度偏差纠正后评价指标结果图Fig.4 Results of the rainfall and temperature after bias correction
以元江站作为流域出口断面,取该站1991-2000年月径流数据。采用四折交叉验证法,根据前文数据的输入因子的划分,分别对3 种模型进行模拟,每个模型进行4 次试验,得到的测试集相对误差和均方根误差及其平均值如表3所示。由表可知,M3 模型4 次试验的RMSE为86、PBIAS为1.7%,远小于M1、M2 模型,而M1、M2 模型的PBIAS均小于20%,因此3 种模型的预测性能和泛化能力均良好。
表3 3种模型的交叉验证评价结果Tab.3 Results obtained from three models calibrated with cross validation
分别以3 种不同预报因子作为模型(M1、M2、M3)输入,对元江流域月径流量进行模拟。BPNN 模型在进行训练时,为了避免局部最优,采用100组随机初值进行训练,择优得到各模型1991-2000年的月流量过程曲线及测试阶段的评价指标如图5所示。从图5(a)和图5(b)的流量过程曲线可以看出,M1 整体拟合度都比M2 高,M1 在丰水期是径流量模拟值比实测值低,M2 在丰水期和枯水期模拟值和实测值都有较大偏差,且洪峰没有模拟出来;而从图5(c)和图5(d)的流量过程曲线可以看出,M3 枯水期拟合程度高,丰水期除了个别模拟值偏高,总体结果比前两者好,且洪峰模拟值偏大对于防灾减灾影响不大;为了验证上述分析结果,使用评价指标从不同角度定量分析,从图5中可以得到,不同指标体现的各模型从优到次排名,其中NSE为M3、M1、M2,PBIAS为M3、M2、M1,RMSE为M3、M1、M2,从排名中可以看到M3模型3个指标都较前两者好,NSE分别比M1、M2 提高10.3%和25%,PBIAS分别比M1、M2 减少了33.1%和17.8%,RMSE分别比M1、M2 减少了3.9%和14.1%。综上所述,说明基于GCMs 集合与实测气象融合数据可以有效模拟元江流域的月径流,后面就用M3 模型来预测元江流域未来两个时期的月流量。
图5 各模型的流量过程曲线及其评价指标Fig.5 The streamflow hydrograph obtained from different models and the corresponding evaluation index
为了分析元江流域未来降雨和气温相对基准期的变化,根据历史实测气象数据和GCMs 集合结果,计算不同路径下多年平均降水和多年平均温度相对基准期的变化值。如表4所示,2041-2070年除SSP370路径下的年均降雨量减少1.8%,其余路径下的年均降雨量都增加。2071-2100年所有路径下的年均降雨量都增加,且最大增幅达到24.6%,同一路径下2071-2100年比2041-2070年年均降雨增幅大;未来两个时期所有路径下的年均温度都增加,从低碳排放量到高碳排放量路径下年均温度变化幅度依次增大,最大增幅为SSP585 下的2071-2100年达到20.4%,同一路径下2071-2100年比2041-2070年年均温度增幅大。总体来说2041-2070年降雨量可能有6%的增加,温度可能有7%~16%的升高;2071-2100年降雨量可能有10%~25%的增加,温度可能有7%~20%的升高。
为了分析元江流域未来径流量相对基准期的变化,利用率定好的M3 结合2041-2100年GCMs 集合结果,预估未来两个时期径流。如表4所示,2041-2070年除SSP370 的年均径流量减少0.6%,其余路径下的年均径流量都增加,2071-2100年所有路径下的年均径流量都增加,且同一路径下2071-2100年比2041-2070年径流量增幅大,最大增幅超过10%。总体来说2041-2070年径流量可能有1%~3%的增加,2071-2100年径流量可能有2%~10%的增加,可以发现径流的变化与降水变化相似。
表4 未来不同排放情景下水文过程变化Tab.4 The changes of hydrological processes under different emission scenarios in the future
分析图6可知,从季节变化来看,总体上不同路径下季节径流预测结果一致,2041-2100年春季径流量减少而秋季径流量增加,夏季无明显变化,2071-2100年冬季除SSP126 外,其余按从低碳排放量到高碳排放量路径下径流量增幅依次增大。进一步分析元江流域未来月径流的变化,如图7所示是未来4 个情景下的月平均径流相对基准期的变化值图,总体上不同路径下月径流预测结果一致,未来1月和4月径流有较明显的减少趋势,2月和11月有较明显的增加趋势,其中SSP585 径流增幅明显大。
图6 未来不同排放情景下的季径流变化Fig.6 The changes of seasonal runoff under different emission scenarios in the future
图7 未来不同排放情景下的月径流变化Fig.7 The changes of monthly runoff under different emission scenarios in the future
元江流域未来径流变化与降水变化呈正相关,与气温变化呈负相关,气象、水文变化关系与赵红玲[34]研究红河流域径流对降雨、气温响应的结论类似。未来径流量增加原因分析如下:元江流域属于干热河谷,如今森林覆盖率不足5%,全是红色沙页岩地层,不透水,生态系统脆弱[35],而本文研究发现气候变化下未来降水量会显著增加,在森林植被少的地方降水对地表径流的响应更明显,因此未来径流量增加是由于降水的变化引起的。未来降雨量增加,而元江流域不透水的沙页岩使得水资源利用率不高,如今元江流域水利化程度为26.2%,是全国水利化程度比较低的地区之一,建议大力修建水利工程用于灌溉耕地和林地。
为了探究未来元江流域水文过程的演变趋势和提供防灾减灾的参考依据,本文基于不同路径情景下5 个GCMs 模型的结果和BPNN 模型,结合M-BCSD 降尺度方法和交叉验证的方法,预测元江流域未来降雨、气温和径流量并分析其演变趋势。得到结论如下:
(1)M-BCSD 方法对GCMs 在元江流域的降水和温度的结果进行降尺度处理效果显著。本文用4个评价指标对历史实测值与GCMs 降尺度后的结果进行比较,CC值分别在0.7 和0.9 以上,RMSE分别小于50 mm 和1.2 ℃,ME绝对值分别小于2.5 mm和0.4 ℃,BIAS值都在3.5%以内,结果表明基于M-BCSD 方法的GCMs降雨和气温结果可靠;
(2)BPNN 可以用于预估气候变化对元江流域径流量的影响。引入交叉验证法通过划分数据集对3 种不同预报因子的BPNN 模型进行训练,结果证明了建立的3 种BPNN 的泛化性能,可以用于模拟元江流域的月径流;
(3)GCMs集合与实测气象融合数据作为输入因子的BPNN模型模拟元江流域的月径流效果最好。本文建立基于三种不同预报因子的BPNN 模型进行模拟,通过比较流量过程曲线和评价指标,融合数据在丰水期和枯水期的流量曲线拟合程度比基于实测气象数据或GCMs集合数据的拟合程度高,RMSE分别减少3.9%和14.1%,NSE分别提高10.3%和25%,PBIAS分别减少33.1%和17.8%;
(4)研究得到元江流域未来水文过程演变趋势:年尺度下,元江流域未来降雨量、温度和径流都增加,分别为24.6%、20.4%、10.2%;季节尺度下:春季径流减少而秋季径流增加;月尺度下,4月径流减少而11月径流增加。不同路径下未来降雨、气温和径流变化趋势大致相同,仅增减幅度有所差别。
本文对元江流域在气候变化下未来径流量预测可为应对极端水文事件提供参考依据,但是气候变化预测往往具有较大的不确定性,使得决策者无法制定有效的防洪减灾措施。本文主要包括了5 个GCMs 模型、4个共享社会经济路径、BCSD 降尺度以及BPNN 模型的不确定性,今后在不确定性分析方面还亟待进一步的研究。□