任梅芳,徐宗学,庞 博

(1. 北京师范大学水科学研究院,北京 100875；2. 城市水循环与海绵城市技术北京市重点实验室,北京 100875)

近年来,变化环境下河川径流发生了明显的变化,在气候变化和人类活动双重影响下,干旱地区流域径流显著减少[1- 2],湿润地区流域径流增加[3],而城市流域的洪涝事件增多增强[4]。变化环境下流域径流演变及归因是国际水文科学协会(International Association of Hydrological Sciences,IAHS)最新10年科学研究计划的核心研究主题之一[5],研究水文过程变化及其驱动机理对评价气候变化和人类活动对全球和区域水循环演变过程的影响十分重要[6]。通常情况下,降水是流域洪水变化最主要的驱动要素[7],在全球气候变暖及城市化高度发展的背景下,全球及区域降水发生了很大的改变,同时,流域洪水的变化也受到其他因素的影响,如气温、地下水变化及人类活动[8]。近年来,有关学者对全球及区域径流的驱动机制进行了分析。Zhang等[9]对中国降水和径流的关系研究发现,不同流域的人类活动伴随着降水的变化对径流表现出不同的影响；Wang和Hejazi[10]研究气候变化和人类活动对美国413个流域径流的影响发现,气候变化导致流域径流量增加,而由人类活动导致的径流变化在不同流域具有不同的特性；王国庆等[11]研究表明,人类活动对中国北方河川径流的影响大于气候变化带来的影响,而在中国淮河以南,由气候变化引起的径流变化占主导地位；Li等[12]研究表明,降水、气温、农业用地和水库的修建对黄河流域流量的变化有显著影响,由水库和农业用地相关因素变化引起流量的变化是错综复杂的,量化每个驱动要素对流量变化的影响仍然面临着较大的挑战。在城市流域中,城市化的发展对城市流域流量的变化起到至关重要的作用。城市化进程使得城市区域人口急剧上升,城市区域不透水面积显著增加,对城市流域的水循环过程造成一定的影响,而在高度城市化地区,多种人类活动共同作用使得城市洪水的发生机理更为复杂。

大多数传统的水文变量统计分析是基于一致性假设基础之上的,换言之,假设水文变量在时间序列上没有发生突变[13- 14]。然而,由于受到气候变化和人类活动的影响,无法保证水文资料在过去、现在和未来都服从同一分布,这就导致了水文变量在时间序列上的非一致性[15]。大量研究表明一致性假设在水资源规划及危险性评估中不再适用,非一致性洪水频率分析已成为前沿水文科学问题[16- 17]。目前,国内常采用还原法或还现法来考虑非一致性水文频率计算[18]。还原法是指将非一致性时间序列修正为满足一致性条件的时间序列,还现法是将水文资料修正为在现状条件下满足一致性要求的时间序列[16],而大量研究表明传统的还原法或还现法存在“还原失真”和“还原失效”的现象[18- 19]。对于非一致性水文资料的分析常采用3种方法,分别为基于混合分布、基于时变矩和基于条件概率分布的非一致性水文频率分析方法。由于时变矩法可以灵活地选取与洪水变化相关的解释变量来描述洪水的变化,并且能够对非一致性进行成因分析,因此,时变矩法成为近年来国内外研究非一致性洪水频率分析最常采用的方法之一。广义加和模型(Generalized Additive Models for Location,Scale and Shape,GAMLSS)是一种引入位置、尺度和形状的广义参数可加模型,该模型能够灵活描述统计参数和解释变量之间的线性或非线性关系,为时变矩法提供了强大和便捷的工具[20- 23]。

北京市作为中国的政治、文化和经济中心,是城镇化程度最高、经济最发达的地区之一[24]。温榆河流域位于北京市的核心区域,该流域流经北京市昌平、海淀、顺义、朝阳和通州5个行政区。近年来,北京市洪涝事件频繁发生,全面科学识别不同等级洪水发生的关键要素,研究气候变化和城市化如何影响城市流域流量的变化对城市防洪工作具有重要的科学意义。

本文以高度城市化的温榆河流域为研究对象,借助GAMLSS模型,对温榆河流域夏季全概率洪水演变的关键驱动要素进行识别,揭示变化环境下城市流域不同频率洪水演变的主要驱动机制,以期对中国海绵城市建设和城市洪水管理提供科学依据。

1 研究区概况

温榆河流域位于北京市中部,是北京市城区主要的防洪和排水河道,全长47.5 km,流域面积为2 478 km2,被誉为北京市的“母亲河”,流域属于温带大陆性季风气候,夏季炎热多雨,冬季寒冷干燥。多年平均降雨量为619.0 mm,大多数雨量集中在6—9月,约占全年的84%。流域位于北京市核心区域,自20世纪80年代以来随着北京市城市化的高度发展,温榆河流域的土地利用发生了显著变化,不透水面积显著增加,城市化率从1985年的约4%增长到2016年的约42%[25- 26]。

2 数据与方法

2.1 数据及来源

本研究采用温榆河流域出口北关闸1985—2016年夏季日尺度流量数据和流域内10个气象站日尺度降水数据,水文、气象站地理位置如图1所示。将流域内10个气象站日降雨量的平均值作为流域平均降雨量,并计算夏季累计日降雨量作为降雨量驱动因子,水文数据来自北京市水文总站。气温数据来源于中国区域高时空分辨率地面气象要素驱动数据集(China Meteorological Forcing Dataset,CMFD),该数据集是中国科学院青藏高原研究所研发的近地面气象与环境要素再分析数据集,CMFD数据集包括降水、近地表气温、地表气压、近地表空气比湿、近地面全风速、向下短波辐射和向下长波辐射等7个气象变量,时间分辨率为3 h,空间分辨率为0.1°×0.1°,数据下载地址为http:∥westdc.westgis.ac.cn/data/7a35329c- c53f- 4267- aa07- e0037d913a21[27]；将CMFD数据叠加转换为日尺度气温数据,计算研究区域季节气温平均值作为流域气温。本研究采用北京市逐年地下水埋深反映温榆河流域地下水埋深状况,该数据来源于北京市水文总站。采用全球逐年人工不透水面积地图来分析温榆河在过去几十年间人工不透水面积的变化过程[28],该数据集是采用Landsat遥感影像数据,借助于Google Earth Engine平台建立的1985—2018年全球不透水面积地图,数据空间分辨率为30 m×30 m,数据来源于http:∥data.ess.tsinghua.edu.cn/。图2为温榆河流域在1985—2016年不透水面积空间变化过程示意图(红色区域为不透水面)。

图1 温榆河流域位置及水文气象站点分布Fig.1 Location of the Wenyu River basin and the spatial distribution of hydrological and meteorological stations

图2 研究期温榆河不透水面积空间变化过程Fig.2 Spatial changes of impervious areas in the Wenyu River basin during the study periods

2.2 研究方法

为了对温榆河夏季不同概率洪水的非一致性进行诊断,本文选用Pettitt突变检验法对温榆河夏季不同概率流量数据进行突变检验,Pettitt突变检验能够较好地识别出数据序列的突变点,其优点为计算简便,可以明确突变的时间[29]；采用GAMLSS模型对研究区不同概率洪水建立模型；最后采用皮尔逊相关系数法评估GAMLSS模型模拟效果及检验优选模型的鲁棒性。

2.3 GAMLSS模型简介

GAMLSS是一种引入位置、尺度和形状的广义参数可加模型,该模型于2005年由Rigby和Stasinopoulos提出[30],GAMLSS模型是假设因变量服从一种参数分布的广义回归模型,这种分布的所有参数可以通过解释变量的函数来估计,因此,GAMLSS模型的主要特点是根据解释变量的值允许不同因变量的分布形式[30]。近年来,国内外学者将GAMLSS模型应用到水文分析中,该模型不仅能够对非一致性的水文时间序列进行模拟分析,并且能够对非一致性进行成因分析。

(1)

式中:θk和ηk为长度为n的向量；βk=(β1k,…,βJkk)T为长度为Jk的回归参数向量；Xk为一个已知的n×Jk的解释变量矩阵；hjk为解释变量的函数,可以较灵活地采用协变量通过所选的分布函数来描述。

2.4 GAMLSS模型建立与优选

表1 GAMLSS分布函数表达式及连接函数形式

表2 GAMLSS模型表达式及协变量因子

为了选取各概率流量的最优拟合模型,引入广义AIC(Akaike Information Criterion)准则GAIC(GAIC,Generalized Akaike Information Criterion)作为评判标准[34]:

GAIC=DG+ωdf

(2)

式中:DG为模型的全局偏差；df为模型自由度；ω为惩罚因子。当ω=2时称为AIC准则;当ω=lgk时称为SBC(Schwarz Bayesian Criterion)准则,k为解释变量的个数；AIC和SBC准则是GAIC准则的2种特例,取GAIC最小的模型作为最优模型,本文选用SBC准则作为评判准则。

由于SBC准则不能较直观地反映出优选模型的拟合效果,GAMLSS模型中worm图可反映出优选模型的残差分布,可作为模型的评判标准。为了评估优选模型的鲁棒性,本研究采用留一法对最优模型的预报流量值进行交叉验证,即为将各个概率各个年份的流量逐一作为预报项,将其余年份的驱动因子作为预报因子,采用最优模型来预报去除年份的流量值,逐一进行这一过程,直到得到所有的流量值,采用Pearson相关系数法,分别计算观测流量值与模型交叉验证预报值之间的相关系数,从而验证模型对不同概率流量值的拟合效果。

2.5 皮尔逊相关系数法

皮尔逊(Pearson)相关系数是描述2个随机变量线性相关的统计量[35],取值范围为(-1,1),假设2个变量x(x1,x2,…,xn)和y(y1,y2,…,yn),2个变量之间的Pearson相关系数(r)计算公式式(3)所示:

(3)

当r>0时,表明2个变量呈正相关,当r越接近1表明正相关越显著；当r<0时,表明2个变量呈负相关,当r越接近-1表明负相关越显著；当r=0时,表明2变量相互独立。

3 结果与讨论

3.1 Pettitt test径流突变点分析

采用Pettitt法对温榆河夏季不同概率流量数据进行突变检验,由于篇幅所限,本文仅列出95%、75%、50%、25%、10%和5%等6个概率流量的Pettitt突变点检验结果(如图3所示)。从图中可以看出,不同概率的流量均存在突变点,95%和75%概率的流量突变点发生在2011年,50%概率流量突变点发生在1996年,25%、10%和5%概率流量序列突变点均发生在1998年。

图3 温榆河夏季不同概率洪水Pettitt突变点检验结果Fig.3 Pettitt test results of different probability discharges in summer in the Wenyu River basin

3.2 GAMLSS模型优选

本研究对温榆河夏季不同概率(仅列出部分概率模拟图)的流量数据建立GAMLSS模型进行拟合分析,GAMLSS模型可对预测变量的概率分布进行拟合,换言之,GAMLSS模型拟合的不是1个单一值,而是全概率分布,图4为各概率流量最优模型拟合的全概率分布(主要计算95%、75%、50%、35%、20%和5%分位数)。从图4中可以看出,温榆河夏季高频流量呈现出逐渐上升的趋势,而中高频流量比较稳定,从图中还可以看出,绝大多数观测流量点据位于模拟值范围内,说明优选模型能够较好地捕捉到观测流量的变化特征。

图4 不同概率洪水最优模型拟合全概率分布Fig.4 Probabilistic distribution of selected best models for different probability of floods

根据SBC准则,将SBC值最小的模型作为最优模型,表3给出了各个概率优选模型残差分布的平均值、方差、偏态系数、峰度系数和Filliben相关系数,全局偏差和SBC值,残差分布的平均值越接近0,方差越接近1,偏度系数越接近0,峰度系数越接近3,Filliben相关系数都大于0.95,证明所选模型的残差分布符合正态分布[35]。由于SBC准则无法直观反映出优选模型的拟合效果,采用worm图来反映GAMLSS模型的拟合效果,图5给出了各个分位数优选模型的残差分布worm图,图中模型拟合残差值位于两条黑色曲线包围的区域,代表模型拟合效果较好,从worm图可以看出,各个概率的流量拟合残差满足要求。为了验证优选模型的鲁棒性,采用留一法对模型拟合结果进行验证,计算观测值与预报值序列的皮尔逊相关系数(如表3所示),经计算,预测值与观测值之间的r平均值为0.74,说明所选最优模型的模拟效果较好。

表3 优选模型拟合效果相关参数及相关系数

图5 不同概率洪水模拟残差worm图Fig.5 Worm plots for different probability of floods

3.3 协变量分析

由于GAMLSS模型位置参数(μ)为反映预测变量量级变化的参数,为了识别温榆河流域夏季不同概率流量非一致性的关键驱动要素,分析不同概率洪水非一致性的成因,探讨该流域全概率洪水演变的驱动机制,本文对μ表达式的系数进行讨论分析。图6为不同概率洪水优选模型μ表达式的系数变化图,从图中可以看出,μ的截距(如图6(a))随着洪水发生概率的降低而增高,这与前人研究结果一致[12,33]；从图6(b)可以看出,发生概率高于70%的夏季小洪水主要受到城市不透水面积比变化的影响,其不透水面积变化比的系数要远大于降水的系数；而发生概率为45%到70%的夏季洪水受到降水和不透水面积比变化的双重影响,但降水的系数要大于不透水面积比的系数,说明这一发生概率阶段的洪水,相比于降水,受到不透水面积比变化的影响较小；对于低于概率45%的夏季洪水,不透水面积比变化不再选为关键驱动要素,降雨的系数越来越大,说明夏季低于45%概率洪水演变的主要驱动因素为降水的变化。

图6 夏季不同概率洪水优选模型μ截距及协变量系数变化Fig.6 Intercept and the coefficients of covariates of location parameters μ of selected models for different probability floods in summer

综上所述,温榆河夏季中小洪水(中高频洪水)的变化主要受到城市化的影响,而温榆河在过去几十年间城市化发展十分迅速,这也就较好地诠释了温榆河夏季中小洪水呈现出逐渐上升的趋势(如图4)；而夏季大洪水(低频流量)变化的主要驱动因素为降水的变化。

4 结 论

由于受到气候变化和高强度人类活动的影响,城市洪水在长时间序列上表现出非一致性特性,本文采用GAMLSS模型对温榆河夏季全概率洪水变化的主要驱动机制进行分析探讨。主要结论如下：

(1) 温榆河流域夏季不同概率的洪水在研究期呈现出非一致性,并且不同等级的洪水出现突变点的年份不一致；GAMLSS模型能够较好地应用到城市洪水分析中,该模型可以较好地捕捉到城市洪水的变化特征,优选模型预测的中位数值与观测流量值之间的Pearson相关系数平均值为0.74,说明模型拟合效果较好。

(2) 降水和不透水面积比的变化是温榆河夏季洪水演变的主要驱动要素,不同等级洪水呈现出的非一致性特性具有不同的驱动成因。温榆河夏季高于概率70%的小洪水变化的主要驱动要素为不透水面积比变化；概率为65%到50%之间的中小洪水变化的主要驱动要素为降水和不透水面积比变化,其中受到降水的影响要显著大于城市不透水面积变化的影响；而对于低于45%概率的大洪水,城市不透水面积变化比不再作为主要驱动要素,而其主要受到降水变化的影响。

(3) 下垫面类型改变的影响主要体现在中小洪水上,而大洪水及极端洪水的演变主要受到降水变化的影响,这一结论可对中国海绵城市建设提供一定的科学借鉴价值。

(4) GAMLSS模型可以灵活地选取与洪水变化相关的解释变量来描述洪水序列的非一致性,并且能够对非一致性进行归因分析,但GAMLSS模型更擅长于评估水文资料在时间系列上的非一致性成因,而在对水文序列的预报方面可能存在一定的缺陷。