城市化背景下年最大日径流演变及影响因素研究——以长江下游秦淮河流域为例*

罗 爽,许有鹏**,王 强,于志慧,林芷欣,唐 仁,江如春

(1:南京大学地理与海洋科学学院,南京 210023)(2:江苏省水旱灾害防御调度指挥中心,南京 210029)

随着全球性的气候变化加剧以及区域性的人类活动增强,极端水文事件日趋频繁[1],城市化地区洪涝问题日益严重[2-3],这严重威胁到了社会生产、人类生活以及生态环境的安全与稳定,因此城市发展与水安全成为了研究关注的焦点。与此同时,变化环境下水文事件存在的非一致性已成为共识[4],检测水文事件的变化特征并量化不同驱动因素的作用对于区域洪涝灾害防御具有重要意义[5]。城市化地区的水文过程发生较大改变。一方面,自然条件下的降水、蒸散发和产汇流等水循环过程受到了显著影响[3],城市地区产生“热岛”与“雨岛”效应[6],径流出现明显变化[7-8];另一方面,土地利用类型变化使得下垫面特征发生转变,进而影响了自然条件下的产汇流,而城市河道渠化及排水系统管网化使洪峰提前,并降低了天然河道的调蓄能力,不断增长的不透水面成为加剧城市化地区洪涝问题的新的重要因素[9],因此亟待开展城市化背景下洪水变化及驱动机制的研究。

随着暴雨洪水问题的加剧,极端径流变化规律越来越受到关注[10],而在统计意义上具有代表年际极端洪水量级变化特征的年最大日径流是探究极端径流演变的重要途径[11]。相关研究发现了不同地区年最大日径流存在的增减趋势以及突变情况[12-14],但是变化机制较为复杂[15]。在城市化地区,虽然年最大日径流的演变对降水量变化的响应较强,但日益加剧的人类活动使得两者之间的关联性降低[16]。已有研究表明,年最大日径流对温度[17]、气候指数[18]和用地类型[19]等因素变化的响应加强,甚至人口可作为解释变量来反映人类活动强度变化对年最大日径流的潜在影响[20]。非一致条件下,对于极端径流演变的驱动机制分析仍然存在较大挑战,而对年最大日径流变化进行归因分析的研究同样欠缺。除此之外,在识别年最大日径流变化影响因子的方法上,基于明确物理过程的水文模型法应用较多[21],但由于模型结构、参数校正和尺度问题的影响,研究结果可能存在一定偏差。一般的统计方法如多元线性回归法、累积双曲线法等则由于原理较为简单而仅适合描述较长时间尺度的平均变化特征,导致在分析年最大日径流变化机制的应用上也有所不足。广义可加模型(generalized additive model for location scale shape,GAMLSS)是包含位置、尺度和形状的广义参数可加模型,它可以根据解释变量的值允许不同因变量的分布式,从而灵活描述统计量与解释变量之间的线性或非线性关系[22]。如L’opez等应用GAMLSS模型对西班牙内陆河流的年最大日径流进行了非一致分析,并选择4项气候因素和1项人类活动因素将它们与分布参数建立非线性关系[23];任梅芳等应用GAMLSS模型分析了北京温榆河流域夏季洪水的演变机制,发现降水和不透水面是关键驱动因素[24]。因此,通过运用该方法选择不同因变量构建水文事件统计模型,可以较好地揭示不同驱动因素对水文过程的影响。

长江下游地带是我国重要的经济活动地区,而随着近年来城市化快速发展,自然河道受到强烈的人为干扰,区域面临严峻的洪涝风险。以往关于城市化地区极端径流演变机制的研究通常以气候变化和人类活动两个角度进行概括性归因分析,缺乏对关键驱动因子的探讨,而极端径流对驱动因素的响应也存在较大不确定性。鉴于此,本文以我国长江下游典型城市化地区的秦淮河流域为例,探讨城市化背景下年最大日径流变化趋势和非一致性特征,选取不同的潜在驱动因素并借助多元线性回归和GAMLSS方法对影响因素建立年最大日径流演变模型,综合判别关键驱动因子并量化其贡献作用,以便为长三角区域防洪减灾和城市水资源环境管理提供一定参考。

1 研究方法与数据来源

1.1 研究区概况

图1为研究区的地理位置、地形地貌、水文站点分布及不透水面现状示意图。秦淮河流域地处江苏省西南部长江的东岸,主体为南京市,流域总面积约2631 km2。由于该流域上游位于丘陵区且支流众多,中下游流经城镇地区地势平坦,汇流速度较快,因此洪涝威胁较大。已有研究主要通过水文模型对秦淮河流域不同气候和下垫面条件下典型洪水事件的变化进行模拟[25-26],结果表明降水和不透水面的改变强烈影响了区域水文特性,但针对该区洪水非一致性特征的研究较为缺乏,尤其是年最大日径流的时序演变及归因分析。该流域在20世纪80年代后期大规模城市化而区域不透水面在2000年后扩张显著[9],因此本文将1987-2018年作为研究的时间段以揭示城市发展对洪水变化的潜在影响。

图1 研究区示意Fig.1 The location of study area

1.2 数据来源

本研究中的年最大日径流为全年日平均流量的最大值,使用的水文数据主要包括径流资料和降水资料,其中径流资料为秦淮河流域出口武定门站和秦淮新河站的逐日径流数据,降水数据为研究区内7个雨量站的逐日实测资料,均来源于长江流域水文年鉴。土壤湿度和气温数据来源于ERA Interim数据集,下载网址为https://apps.ecmwf.int/datasets/data/interim-full-daily,获取到1987-2018年的近地表(2 m)气温和土壤湿度数据,时间分辨率均为12 h,空间分辨率均为0.125°×0.125°,由于土壤湿度数据包含3个层位则取平均值代表该地区的土壤湿度。不透水面数据来源于全球不透水面(global artificial impervious area,GAIA)数据集[27](http://data.ess.tsinghua.edu.cn/gaia.html)。

2 研究方法

2.1 趋势及突变分析方法

本研究综合采用了线性回归趋势方法、Mann-Kendall趋势及突变分析(以下简称M-K方法)、Pettitt突变检验等方法,检测年最大日径流的长序列趋势性和突变性特征。M-K方法和Pettitt突变检验是广泛应用于气象学和水文学领域的趋势和突变分析方法,结合两种方法共同确定突变分析的结果,详细介绍及具体原理见参考文献[28-29]。

2.2 驱动因素选取

参考国内外洪水演变和归因分析的相关研究,针对常见驱动因素和表征形式进行归纳得到表1。首先,降水和气温是反映气候变化的两个重要指标[19]。对于研究区而言,全年的降水集中于汛期因而汛期降水量级对于流域的防洪排涝等实际工作有着直接的影响,因此选取汛期降水量可以充分反映流域的降水结构变化对年最大日径流的作用。与此同时,流域气温因素的年际变化可以通过汛期近地表平均气温进行表征。其次,土壤水分对于流域的产流存在一定作用,因而降雨过程发生前无雨日的平均土壤湿度(前期土壤湿度)选为一项背景因子。此外,土地利用方式和人口规模等在统计上存在较大的不确定性,而城市化地区的不透水地表变化是人类活动对流域下垫面影响的典型体现且其变化特征稳定[9,24]。综合考虑了气候变化和研究区近年来下垫面的改变等实际情况,同时为减少模型冗余带来的不确定性,本文选定了汛期降水量、汛期均温、前期土壤湿度和不透水面率4类指标作为年最大日径流变化分析的潜在驱动因素。

表1 洪水变化分析中常见驱动因素及表征形式Tab.1 The common driving factors and their representations in flood change analysis

2.3 多元线性回归方法

多元线性回归法是一种较为常用的多因子拟合方法,可使用该方法分析水文事件对不同因素存在的线性变化响应和各因素作用程度。首先将所有潜在驱动因素序列值采用归一化处理,使得所有因素的数据序列在相同范围,消除量纲差异以满足可比性,之后通过因素的组合建立如表2所示的多元线性回归方程。其中,一类模型是为了描述所有因子共同对径流量级变化的作用,二类模型是描述降水和其他任意2个因子组合的作用,三类模型是描述降水和其他任意1个因子组合的作用,四类模型是描述降水作为单一因子的作用。最终,根据R2与校正R2(adjustedR2)的值选取出最优模型。参考相关研究[30],通过计算最优模型中各因子回归系数的权重即可量化各因素对因变量作用的相对贡献率。

表2 模型表达式*Tab.2 The expression of different models

2.4 GAMLSS模型

多元线性回归方法虽然可以对年最大日径流序列的变化进行模型拟合,但是无法拟合变量的非线性变化响应且因素权重赋值不够客观,导致确定性有所欠缺[31],因此本研究引入GAMLSS方法,通过结果的比较验证,进一步明确不同潜在驱动因素对年最大日径流事件的作用。

GAMLSS模型假设随机观测值yi(i=1,2,…,n),服从分布函数F(y|θi),θi=(θ1,θ2,…,θp)表示具有p个参数(位置、尺度和形状参数)的相邻。fk(θk)表示解释变量θk与随机响应变量Xk之间的函数关系式:

(1)

式中,ηk和θk是长度为n的向量;Xk是一个n×Jk的解释变量矩阵;βk是长度为Jk的回归参数向量;hjk表示分布参数和解释变量xjk之间的联系函数,本研究为三次样条函数(cubic spline,CS)。

先将所有潜在驱动因素序列值采用归一化处理,之后在R软件里使用“gamlss”程序包对年最大日径流序列建立GAMLSS模型。根据前人研究结果[18],本研究选用5类分布函数,包括Gamma(GA)、Generalized Gamma(GG)、Gumbel(GU)、Lognormal(LOGNO)以及Weibull(WEI)来拟合因变量对不同因素的响应,5种分布函数的密度函数和连接方式见参考文献[22]。由于GAMLSS模型为广义的回归形式,而模型中的μ参数是表征预测指标量级变化的参数[19],因此该参数中的协变量系数可充分反映自变量对因变量的作用。本研究设置GAMLSS模型时首先将μ参数里各因素按表2进行组合,5类分布函数每类对应8个,总计得到40个模型。之后,基于拟合的AIC(akaike information criterion)值最小的原则选择最优模型。最后,提取最优模型μ参数中的协变量系数,计算其权重得到因素作用的相对贡献率。

此外,多因素的综合作用具有边缘效应[32],对于本研究而言, 年最大日径流对汛期降水量因素的响应可能也会随着不透水面率的变化而发生一定改变。为了更清楚得到不同阶段年最大日径流对降水的响应以及不透水面变化的具体作用,确立最优GAMLSS模型后,通过设立汛期降水量和不透水面率的相互作用项(x1·x4)来代替表1中单独设置的不透水面率因子(x4),之后可通过下式来计算边缘效应[32]。

(2)

式中,ΔD为边缘效应;αi(i=1,2,3,4)分别为各因素和相互作用项在μ参数中的协变量系数,j为最优模型中除相互作用项外的因素个数;t1和t2为年最大日径流序列起止或发生特征性变化的时间节点,根据实际情况确定。

3 结果分析

3.1 年最大日径流演变规律

秦淮河1987-2018年最大日径流序列如图2所示,使用M-K和Pettitt检验分析方法所得的结果如图3所示。综合两种突变检验方法的结果可知,秦淮河流域年最大日径流序列在2001年发生显著突变,其演变趋势上存在非一致性。秦淮河流域年最大日径流在突变前为减少趋势(-9.62 m3/(s·a)), 而后为增加趋势(23.41 m3/(s·a)), 但均不显著。突变前平均值为573.28 m3/s, 而突变后为820.47 m3/s,且对于不同阶段中的低值来说,突变后明显高于突变前,说明年最大日径流有明显增大趋势。序列整体线性变化速率为14.77 m3/(s·a)(P<0.05),因此可知秦淮河流域1987-2018年的年最大日径流呈显著增加趋势。

图2 秦淮河流域1987-2018年最大日径流序列Fig.2 Trend of annual maximum daily runoff in Qinhuai River Basin from 1987 to 2018

图3 年最大日径流序列突变检验结果Fig.3 Mutation test results of annual maximum daily runoff series

3.2 年最大日径流变化的潜在驱动因素分析

为了准确识别年最大日径流变化的关键驱动因素,探究不同环境变量对其演变的影响,本研究对潜在驱动因素(汛期降水量、前期土壤湿度、汛期均温和不透水面率)的变化特征进行分析。表3为各驱动因素的线性变化趋势斜率,以及使用皮尔逊相关性分析得到的各因素与年最大日径流之间的相关性系数。由表3可知,在与年最大日径流的相关性上,汛期降水量相关系数最大,其次为不透水面率,且两者均呈现显著的相关性,汛期均温的相关系数最小。汛期降水量呈现不显著的增加趋势,汛期均温和不透水面率均呈现显著的增加趋势,前期土壤湿度为不显著的减少趋势。根据不透水面的时空变化(图4)可知,2000年后流域地表硬化加剧,总体不透水面率呈增加趋势并从1987年的2.5%增到了2018年的24.1%。

表3 潜在驱动因素变化趋势及与年最大日径流相关性Tab.3 Trend of potential driving factors and its correlation with annual maximum daily runoff

图4 秦淮河流域不透水面时空变化Fig.4 Temporal and spatial variation of impervious surface in Qinhuai River basin

3.3 年最大日径流变化的驱动机制

根据建立的多元线性回归方程,提取能反映不同类别模型稳定性和可靠性的R2和adjustedR2(表4)。由表4可知,模型1的R2最接近1,7号模型的adjustedR2最接近1而其R2与1号模型相差无几。因此,综合来看汛期降水量和不透水面率两者的组合为最优模型,这也表明加入不透水面率变化的影响比仅考虑降水量变化的影响(模型8)对于模型的拟合效果有明显的提升。

表4 多元线性回归模型拟合结果Tab.4 Fitting parameters of the multiple linear regression model

根据GAMLSS模型拟合的实际情况,提取其AIC值绘制得到图5。遵照AIC最低为最优拟合原则,可知秦淮河流域年最大日径流变化的最优GAMLSS模型为使用LOGNO分布函数的模型7。该结果与多元线性回归模型结果相印证,表明包含不透水面率与汛期降水量变化的模型是较优的选择,由此可知秦淮河流域年最大日径流变化的关键驱动因素是汛期降水量和不透水面率。由于GAMLSS模型可对预测变量的概率分布进行拟合,即GAMLSS 模型拟合的不是一个单一值而是全概率分布[24],因此根据所选定的最优模型绘制拟合分位数图(图6)。由拟合分位数序列可知,绝大多数观测流量点据位于模拟值范围内,说明优选模型能够较好地拟合出不同概率下的年最大日径流,准确地反映其变化和分布特征。

图5 GAMLSS模型AIC值 Fig.5 AIC value of different GAMLSS models

图6 最优GAMLSS模型拟合序列分位数Fig.6 Quantile diagram of the optimal GAMLSS model fitting sequence

依照所建立的最优多元回归模型提取关键驱动因素和回归系数,而最优GAMLSS模型提取μ参数中的协变量系数并计算因素的相对贡献率,结果如表5所示。由表5可知,经两种方法最优模型拟合的汛期降水量和不透水面率的系数均通过了P<0.05水平的显著性检验,表明两者产生的变化对于年最大日径流演变的影响是显著的。对于秦淮河流域1987-2018年的年最大日径流演变过程,多元线性回归方法所得的汛期降水量贡献率(77.1%)大于GAMLSS所得结果(73.7%),相应的不透水面率的贡献率(22.9%)则低于后者(26.3%)。由此可知,选用不同的统计拟合模型对于结果产生了影响,但是两者相差并不大,说明结果得到了相互验证,同时这也指示了降水和不透水面的增加可能会导致年最大日径流增长,它们的共同影响起着重要的驱动作用。

表5 最优模型中驱动因素的系数及贡献率Tab.5 Coefficients and contribution of the drivers in the optimal model

此外,一般在分析径流变化响应时往往将各因素视为相互独立的变量,而本研究为了进一步探究不透水面率的变化在年最大日径流对汛期降水量变化响应中的影响,将降水和不透水面率的共同作用视为交互作用项来代替单独的不透水面率因素,之后提取最优GAMLSS模型协变量系数。同时,考虑到秦淮河流域年最大日径流在2001年发生显著突变,故将(2)式中的t1和t2分别设为研究起止时间和突变发生时间,得各阶段边缘效应。年最大日径流发生突变前(1987-2000年)、发生突变后(2001-2018年)以及整个阶段(1987-2018年)期间边缘效应ΔD分别为0.067、0.104和0.164,这指示着不透水面率的增加改变了年最大日径流和汛期降水量的一致性与相关性,而上述3个阶段中影响程度分别为6.7%、10.4%和16.4%。与此同时,突变发生后的边缘效应(10.4%)大于突变发生前的(6.7%),表明年最大日径流对降水因素的响应在秦淮河流域不透水面大规模扩张后出现较大改变,地区水文事件的非一致性有所增大。

4 讨论与结论

4.1 讨论

相关研究指出,近几十年来秦淮河流域径流显著改变,径流深在2001年左右突变且21世纪以来年均径流深相较于1987-2000年增长了50%以上[33],洪水洪峰流量呈现逐年增加的趋势[9]。本文在此基础上发现,1987-2018年期间秦淮河流域年最大日径流呈现显著上升趋势,而且在2001年左右发生显著改变,表明年最大日径流的演变存在明显的非一致性特征,充分证实区域面临的洪水威胁增大。

在环境变量影响的分析方法上,以往研究多侧重于驱动因素的选取,如Li等[19]、López等[23]、任梅芳等[24]、夏露等[34]通过构建多种因素的多元线性回归或广义可加模型进行径流变化归因分析。但是仅考虑相关关系建立起的多因素统计模型具有局限性,尤其是回归变量存在共线性以及系数难以精确的问题[35]。因此本文在前人研究基础上做出一定改进,包括采用了多元线性回归和GAMLSS模型两种方法综合分析以减少方法选择对结果造成的不确定性,也相应地对不同因素进行了组合并对模型进行了优选。结果表明,不同因子的组合对模型解释力产生了影响,建立适宜的多因子模型可较好地进行分析。

本文所得结果与目前有关城市化地区年最大日径流演化的研究结论一致,证实了降水和不透水面率是重要的影响因素。首先,年最大日径流演变受到降水变化的决定性作用,如林娴等通过HIMS模型计算得到气候变化对广东武江流域年最大日径流增加的贡献率为94%[21]。本研究当中,得到汛期降水量改变对秦淮河1987-2018年期间年最大日径流变化的贡献超过70%,由于较大空间尺度下降水的增多尤其是汛期降水增多导致了流域产流量增加,进而引起年最大日径流的量级显著变化。其次,大量研究指出不透水面率同样是城市化地区洪水变化分析中不容忽视的要素。Shao等[36]通过对武汉市设置多级流域径流监测模型发现,区域不透水面率达20%时的流量是4%时的两倍以上;孙延伟等[9]运用水文模型分析了秦淮河流域不透水面扩张对于洪水的影响,发现大小洪水的洪峰均随不透水面率的增加而显著增加。由于城市化地区剧烈的人类活动,大量天然下垫面被人造的不透水面取代,这阻碍了地表水下渗,切断了城市区域地表水与地下水之间的水文联系[3],致使降雨过后地面汇流速度变快,这对洪峰流量起到增强效应。本文通过最优拟合模型揭示了秦淮河流域洪水演变的重要特征,计算表明不透水面率贡献超过了20%,证实了该地区不透水面的增加与年最大日径流存在的显著关联,它对年最大日径流的增加趋势存在重要影响。此外,值得注意的是,研究区前期土壤湿度的变化并不显著,因此它对于年最大日径流增加的作用可能也会受到不透水面变化的间接影响,如下渗减少和土壤退化压实等现象都导致了秦淮河流域土壤水动态响应规律发生改变[37],进一步影响了该地区的产流过程。还有研究指出,21世纪以来秦淮河流域雨洪一致性关系显著改变[38],本文中前后阶段边缘效应(1987-2000年,6.7%;2001-2018年,10.4%)的差异则证实了该结论,表明不透水面的作用在增大,这也反映了气候变化和人类活动的耦合作用进一步深化。城市化改变了原有的降水-径流响应关系,城市建筑用地的增加可能导致流域的气温、降水以及蒸散发量等气象因素改变进而影响了水循环过程[3]。

4.2 结论与展望

城市化背景下,人类活动和下垫面条件的变化深刻影响了水文事件和水文过程。本文针对秦淮河流域1987-2018年的年最大日径流的演变过程,基于多元线性回归和GAMLSS模型综合识别了关键驱动因素并量化其贡献作用。主要结论如下:

1)秦淮河流域1987-2018年期间年最大日径流呈现显著的增长趋势,平均线性增长速率为14.77 m3/(s·a);年最大日径流在2001年发生显著突变,突变前呈现减少趋势,突变后呈现增加趋势,且突变后的流量平均值远大于突变前,表明近年来区域洪水威胁加剧。

2)以汛期降水量、前期土壤湿度、汛期均温与不透水面率为潜在驱动因素,使用多元线性回归和GAMLSS方法构建了多种年最大日径流演变分析模型,通过因素的组合和模型的优选最终确定了汛期降水量和不透水面率为关键驱动因素,两者与年最大日径流的相关系数分别为0.84和0.33(P<0.05)。

3)根据所构建的最优模型计算得汛期降水量贡献率超过70%,不透水面率贡献率低于30%,表明了降水的增加决定了年最大日径流的上升趋势而城镇建设用地的扩张对年最大日径流的增长有重要作用。2000年之后不透水面率因素对年最大日径流与汛期降水量因素的一致性关系的影响加大,表明了近年来快速城市化过程对区域的雨洪响应关系造成了显著改变。

本文揭示了长江下游典型快速城市地区秦淮河流域的年最大日径流的演变趋势,基于两种方法定性且定量分析了非一致条件下影响其变化的关键驱动因素,有助于城市地区洪水演变及驱动机制相关研究的开展。但是,由于径流响应的复杂性,如何更好地通过构建多因子变量模型厘清不同因素对水文事件的影响是需要进一步考虑的问题。此外,在城市化背景下如何更全面准确地分析人类活动的特征和作用还有待进一步深入的研究来进行揭示。