徐鹏程，刘楠楠，李 帆，仇建春，张昌盛，戴 笠
（1. 扬州大学水利科学与工程学院，江苏 扬州 225009；
2. 江苏省水利建设工程有限公司，江苏 扬州 225009；
3. 泰州市港航事业发展中心，江苏 泰州 225300）
气候变化和人类活动加剧将改变全球水文循环的现状，对降水、蒸发、径流、土壤湿度等关键水文气象要素造成直接影响，导致流域性大洪水的发生风险增加，而这些正成为制约国民经济的主要瓶颈［1-3］。现行的洪水频率分析是建立在平稳性假设前提下的，这意味着影响洪水特征变量的环境控制因素诸如气候要素、土地利用率在过去、现在和将来都遵循一成不变的物理机制［4］。而在气候变化、城市化等因素的多重影响下，基于平稳性假设的洪水频率分析方法的适用性受到了挑战［5，6］。为了充分保证洪水设计值的可靠性，探索非平稳性条件下的洪水频率分析方法以及工程设计值的推求显得至关重要。
目前，非平稳水文频率分析方法分为单变量非平稳和多变量非平稳分析方法。其中，单变量非平稳性已经逐步纳入到洪水极值事件［7］、干旱极值事件［8］和极端降雨事件［9，10］频率分析中。Rigby 和Stasinopoulos［11］提出了考虑位置、尺度、形状参数的广义可加模型（Generalised Additive Models for Location Scale and Shape， GAMLSS）， 其基于参数回归模型定量解析了极值序列的统计参数与物理解释变量间线性/非线性关系。GAMLSS模型涵盖了多种随机变量分布函数类型，克服了传统的分布函数对于特定分布类型（超峰度、高度正偏/负偏和平顶峰度类型的分布）拟合过度问题［12］。Chen 等［13］基于台湾地区降雨极值序列对比分析了单变量情景下3 种非平稳频率（非参数离散小波变换（Discrete Wavelet Transform， DWT）、经验模态分解法（Ensemble Empirical Mode Decomposition， EEMD）以及参数加权最小二乘法）分析方法的适用性。
传统的洪水频率分析主要聚焦于洪峰量级的分析，由于洪水事件（过程）并不能凭借单个特征变量（洪峰）得以完全概化，一次洪水过程往往包含洪峰、洪量、历时等要素信息，是一个包含频域、时域和空间域的复杂水文过程，仅从单变量角度无法有效刻画峰高量大的流域型洪水风险。为此，近年来，一些水文研究也逐步关注到多变量水文序列的非平稳性问题［14］:Sarhadi 等［15］提出了一种贝叶斯动态Copula 模型有效模拟了多维水文气象极端事件的混合连续与离散的多变量时变相依结构。Jiang 等［16］采用动态VineCopula 函数有效建立了西江流域的年最大日径流-年最大3 日洪量-最大7 日-最大15 日洪量的四维联合分布模型，为非平稳条件下多元水文极值设计提供了理论依据。Xu 等［17］基于动态Copula 函数量化分析了不同物理解释变量对洪水设计值的影响，并提出了变化环境下多变量洪水风险的短期预估的一般框架。
淮河流域人口密集，人水争地矛盾突出。因人类活动（城市化）、气候变化引发的洪水极端事件［7］，极端降雨事件［18］都呈现一定的非平稳性， 但是目前对于淮河流域洪水频率的非平稳研究主要聚焦于单变量角度开展的。过去几十年以来，城市化对流域洪水特性的影响逐步成为气候变化领域的研究热点［19，20］。由于多种因素的相互作用，高度城市化地区的洪峰流量增幅较为明显［21-23］。城市化地区不透水面积的增加很大程度上弱化了自然渗透量，导致暴雨径流量显著增加；
同时，城市规划设计层面考虑高效的人工排水沟替代天然河道，减少了径流响应的滞后时间。为此，如何从非平稳性入手，定量解析城市化和气候变化相关的驱动因子与多变量洪水风险的响应关系，是目前非平稳频率分析领域亟待解决的关键科学问题。为此本文基于淮河流域蚌埠（吴家渡）水文站的洪峰（P）洪量（V）序列为研究对象，构建非平稳条件下P、V相依性结构的二维动态Copula 模型，定量分析气候变化因子（极端降雨因子、长尺度气候涛动因子）、城市化因子（不透水面积率（Percent of Impervious Surface Area， PISA））对单变量、多变量洪水风险的影响程度。
由于不受限于边缘分布的形式，Copula 函数是一种备受推崇的水文多变量分析模型，主要用于拟合随机（极值）变量之间的相依性结构。针对气候变化和人类活动导致的单变量、多变量水文极值序列的非平稳性，本文拟构建如下的动态Copula 函数模型:假设［X，Y］为洪水-洪量的极值序列对， 基于时变Copula和边缘分布的多变量联合概率分布模型可以如下公式求得。
以上所有这些参数和时间变量或者物理解释变量呈现线性或非线性函数关系；
物理协变量因子集合主要包括气候涛动因子（NAO、SOI、NINO3.4）、极值降雨量和城市化率；
u和v是边缘分布的累积概率值。在本研究中，非平稳多变量风险评估模型分为2 个个阶段:①采用边缘分布和Copula 为主的时变矩多变量模型，并结合对数似然比法（LR）和最小化AICc（修正化的赤池信息）准则，从分布参数角度定量识别洪水极值序列中潜在非平稳性；
②基于不同非平稳分布模式下特定分位数的风险值的计算定量剥离城市化、气候变化因素对多变量洪水风险的影响。
1.1 非平稳边缘分布模型构建
本节讨论5种常见的极值分布函数为原型如何构建以物理相关因子作为分布参数的解释变量非平稳模型。表1 显示了5种常用极值分布的概率密度函数（PDF）形式。由于篇幅限制，本处着重以GEV分布的位置参数（μ）为例展示非平稳条件下的时变分布模型的形式:
表1 5种潜在的边缘分布的概率密度函数Tab.1 The probability density function （PDF） for the five potential marginal distribution
式中:Covk(t)为表示第k个潜在的时变物理变量序列｛大尺度气候涛动因子（NAO， SOI， NINO3.4），气象因子［极值面降雨量（Maxpre）］，城市化因子（PISA）｝。
式中采用了线性和非线性趋势型非平稳特征。当然多种因子的复合影响，可以增加公式（2）中的协变量因子。考虑到不同分布对应的参数数量存在差异，为此对于三参数分布（GEV 分布和PIII 分布），潜在非平稳边缘分布模型种类明显多于两参数分布的潜在非平稳模型种类。
1.2 非平稳Copula模型构建
为了模拟非平稳条件下洪峰和洪量之间的时变相依性结构，本节讨论物理解释变量作为参数协变量的动态Copula 的建立过程。在多变量水文频率分析中，可以考虑如下多种潜在的二元Copula 函数作为动态Copula 模型的母函数:5 种单参数Copula（Joe、Frank、Gumbel、Clayton 和Gaussian）和5 种双参数Copula［Clayton-Gumbel（BB1）、t、Joe-Gumbel（BB6）、Joe-Clayton（BB7）和Joe-Frank（BB8）］。以上Copula 模型的计算拟合可以通过R语言环境下CDVine包［24］中实现。
类似于非平稳边缘分布模型的构建过程，非平稳Copula 模型同样是通过构建参数与物理解释变量间的线性函数关系，可以建立如下公式:
式中:Cov1(t)，Cov2(t)，…，Covk(t)为k个潜在的时变物理变量序列。
术前常规使用0.5%左氧氟沙星滴眼液每日4次，连用3 d，术后使用0.3%妥布霉素地塞米松滴眼液和0.1%玻璃酸钠滴眼液每日4次，1周后停用妥布霉素地塞米松滴眼液，0.1%玻璃酸钠滴眼液连续使用1个月。
1.3 基于LR-AICc准则的最优非平稳模型的优选
一旦构建了多种潜在的非平稳边缘分布和Copula 函数模型，接下来需要综合运用极大似然比准则（LR）和AICc准则获得最优的非平稳/平稳模型。
假设ST0代表平稳模型，ST1代表特定的非平稳时变模型；
同时LL0和LL1为平稳和非平稳模型的对数似然值，其对应的对数似然比值可以按如下公式求得:
考虑到可能存在多个非平稳性模型能够通过以上的LR准则检验，需要计算各个非平稳模型的修正化赤池信息量［25］（AICc）:
式中:k代表非平稳模型对应的参数个数；
n代表研究序列的样本长度。AICc准则运用原则为:越小的AICc值代表了模型的拟合效果更好。
为了定量剥离城市化因子（PISA），气候变化因子（Maxpre，NAO，SOI，NINO3.4）对洪水风险的影响，因此可以通过重现期得出单变量和多变量风险。
对于单变量情形，每个洪水极值（P或V）的单变量重现期为:
因而单变量风险:
式中:FX(x0)为边缘分布的累积概率值（CDF）；
RP为单变量重现期；
UR为单变量风险（累积超越概率）。
AND情景［7］下的双变量联合重现期（JRP）可计算如下:
式中:FX(x0)，FY(y0)是洪水变量对应的边缘分布累积概率值；
JRP代表联合重现期；
BR是多变量风险（累积联合超越概率）。
淮河流域位于东经111°55"～121°25"、北纬30°55"～36°36"，东西长约700 km，南北宽约400 km，面积27 万km2。流域地跨湖北、河南、安徽、江苏、山东五省，覆盖40 个市、163 个县（市、区）的全部或部分区域。淮河流域地貌具有类型复杂多样、层次分明、平原地貌类型丰富的特点。流域西部、南部及东北部为山区和丘陵区，其余为平原、湖泊和洼地。淮河流域的植被具有明显的过渡性特点。流域北部的植被属暖温带落叶阔叶林与针叶松林混交林；
中部低山丘陵区属亚热带落叶阔叶林与常绿阔叶林混交林；
南部山区为常绿阔叶林、落叶阔叶林与针叶松林混交林，并夹有竹林，山区腹部有部分原始森林。淮河流域地处我国南北气候过渡地带，气候变化复杂，降雨时空分布不均。流域内众多支流多为扇形网状水系结构，洪水集流迅速。由于淮河流域位于多重过渡地带的孕灾区，特定的地理环境和下垫面条件、复杂多变的气候，加之黄河夺淮和人类活动的不利影响，导致流域水、旱、风暴潮等自然灾害的频繁发生。淮河流域的洪水呈现高洪峰-高水位-长历时的特征，易于发生全流域性大洪水。本研究选取淮河流域中游蚌埠（吴家渡）水文站（图1）作为研究对象，基于站点1959/1/1-2016/12/31 的日径流序列提取洪峰-洪量序列并构建动态Copula 函数风险分析模型。由于本文采用不透水面积率（PISA）作为物理协变量纳入到非平稳模型的参数中，所使用的PISA（城市化率）序列可以从如下的网址链接中下载:http://data.ess.tsinghua.edu.cn/。PISA的逐年演变图见图2。
图1 研究区域水文气象站点分布图Fig.1 Map of location of hydrometeorological stations in the study area
图2 淮河中游子流域（蚌埠水文站所控制的子流域）的城市化范围演变图Fig.2 The evolution map of urbanization scope of Huaihe River middle reaches sub-river basin （controlled by Bengbu hydrological station）
多个长尺度气候涛动因子序列［南方涛动因子（SOI），北大西洋涛动因子（NAO）和厄尔尼诺因子（NINO3.4）］由如下链接下载:https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries。同时基于蚌埠水文站控制的子流域边界内的气象站1959-2016年的日降雨序列提取子流域内的极值面雨量（Maxpre）序列作为非平稳模型潜在的气候变化因子。
4.1 非平稳优化结果
首先依据LR-AICc准则对潜在平稳/非平稳模型进行优化选择，具体模型优化结果见表2。对于蚌埠水文站提取的洪峰序列，非平稳模型GEV3 其对应的LR检验的p值为0.003 2 显著小于0.05 的阈值，表明洪峰序列的趋势型非平稳性特征在5%显著水平上能够被检测到；
尽管GEV1-GEV2 模型的p值可以通过LR检验，但是考虑到GEV3模型的AICc值相对其他两个非平稳模型更小，所以GEV3 是洪峰序列对应的最优边缘分布模型，其中位置参数和面极值雨量（Maxpre）呈现二次型统计学关系（非线性），尺度参数与NAO 因子呈现线性关系，形状参数和城市化率（PISA）呈现线性关系。按照同样的方法，模型GA3为洪量序列的最优非平稳模型，而由于AIC值最小化准则选定平稳Frank Copula（FR0）模型为最优的Copula模型。
表2 最优的非平稳边缘分布和Copula模型优选结果Tab.2 Results of best-fitted marginal distribution and Copula models
4.2 LR-AICc 准则对于趋势型非平稳特征检验的有效性
分别对蚌埠水文站的洪峰、洪量进行非参数的单变量、多变量Mann-Kendall（M-K）趋势型检验［26］，其检验结果为:ZQ=1.156，ZV= 0.689。由M-K 检验结果可知洪峰，洪量序列在5%显著水平上趋势不显著。由于M-K 检验侧重于变量均值的潜在的单调性趋势［27］，对于某一分布（如GEV 分布），将存在这样的特定情形:位置和尺度参数的反向趋势会使得M-K检验结果为不显著（|Z|＜1.96）， 这在一定程度上体现了本文采用的对数似然比（LR）检测方法的必要性。为了探究引发上述非参数检测法和LR检测法结果的不一致性，图3 绘制洪峰、洪量分布参数的随时间的演变关系。由于GEV 分布的尺度、形状参数（斜率分别为-0.196 7和-0.001 4）呈现下降趋势性而其位置参数呈现上升趋势性（斜率为6.474 1），由于分布参数间存在的反向趋势性导致了洪峰序列在M-K 检验时趋势不显著而LR检验时趋势性显著；
同样的结果适用于洪量序列，其尺度和形状参数也存在反向趋势性（斜率分别为0.004 3，-0.003 9）。
图3 分布参数趋势的协同性对极值序列非平稳性检测结果的影响Fig.3 The influence of the trend synergy of distribution parameters on the detection of nonstationarity of extreme value series
4.3 非平稳条件下多变量洪水设计值推求
由于洪峰、洪量序列存在一定的趋势型非平稳特征，为此特定联合重现期下的洪峰-洪量的设计分位数对存在一定的时变特性（图4）。虽然为Q-V之间的相依性结构的最优拟合模型为FR0 模型，表明其没有显著非平稳性，但在50 年一遇的联合重现期（50-JRP）水平（或图5 中的BR=0.02）下，Q-V的设计分位数等高线是逐年变化的。由于样本量为58（1959-2016），每年可以生成这两个极值变量在50-JRP水平上的等值线，这将绘制58 条JRP等值线。为了直观呈现设计分位数的时变性，1960 年、1970 年、1990 年、2000 年、2010 的JRP等值线如图5 所示。图5 中综合考虑气候变化、城市化因素的情形对应的50-JRP等值线表现出3 种时变特征:①1960-1970 年间， 洪峰、洪量的设计值呈增大趋势；
②1970-1990 年间呈现下降趋势；
③1990-2010年呈现逐渐增大趋势。单一考虑气候变化或者城市化因素都会出现一定程度上的高估或低估洪峰-洪量的设计分位数。
图4 考虑不同因素影响下的时变联合重现期等值线图（JRP=50）Fig.4 Contour map of time-varying joint return period considering the influence of different factors（JRP=50）
图5 城市化、气候变化对洪水单变量、多变量风险的影响Fig.5 Impacts of urbanization and climate change on flood univariate and multivariate risks
4.4 不同协变量影响下的单变量、多变量风险值对比分析
为了定量剥离城市化（PISA）、气候变化（Maxpre+NAO）对于单变量、多变量风险的影响，依据公式（7）和（9），考虑经验频率值分别为［70%，90%］（步长为0.01）对应洪峰、洪量分位数集合，获得不同非平稳模型（GEV1-GEV3， GA1-GA3）下对应的单变量（UR）和多变量风险值（BR）。为了量化城市化对于单变量和多变量风险的影响，可以分别依据GEV3（综合考虑了城市化+气候变化因素）模型和GEV2模型（仅仅考虑气候变化因素）中分别计算70%～90%分位数区间｛（Q70%～Q90%)， （V70%～V90%）｝对应的单变量和多变量风险值［22，28］。例如对于经验分位数Q70%，对于GEV3和GEV1 可以分别计算两个不同的单变量风险值（URGEV1-Q70%，URGEV3-Q70%）（如图5中左上角回归线图）。接着可以对风险样本序列（URGEV1，URGEV3）进行线性拟合，其中城市化对单变量风险的影响可以通过回归曲线斜率和对角线函数斜率之间的相对差值百分率来定量表示（对角线斜率等于1）。由图5可知，对于蚌埠水文站控制的子流域，城市化对于洪水单变量、多变量风险的影响明显弱于气候变化因素对于洪水风险值的影响。
由于频率分布模型地建立是基于非平稳性假设基础之上，为此采用有效的非平稳性诊断/识别方法是进行非平稳频率分析的必要前提。非平稳性诊断是从数理统计角度判断水文极值序列是否存在趋势的特征。目前，常用的水文序列趋势诊断方法是基于非参数法的Mann-Kendall（M-K）和Spearman。同时对于水文极值序列的趋势性检测在水文气象极端事件领域开展了大量的研究，如洪水趋势分析［29，30］、极端降雨趋势分析［31，32］，干旱趋势分析［3，33］。本文采用的LR-AICc趋势型非平稳检测方法相比于非参数的M-K检测法，有效解决了序列的非单调性趋势 （二次型的趋势:例如先上升后下降型）和参数的反向趋势性问题，为后续非平稳洪水频率分析模型的构建奠定良好的基础。
由于多变量水文序列的概率分布可以分解为单变量的边缘概率分布和变量间的相依结构描述两方面，非平稳性频率模型构建同样分为基于边缘分布的非平稳模型和基于联合分布的非平稳模型。尽管以物理因子作为洪水序列统计参数的解释变量较时间变量能够更好地描述与预测洪水序列中的非平稳性，但是洪水频率分析最终还是以洪水风险的计算作为最终落脚点，同时定量量化各个驱动因子对于洪水设计值的影响程度有利于解析洪水过程的驱动演进机制，本研究采用不同物理协变量组合下不同非平稳模型对比分析的思路有利于定量剥离气候变化因子、人类活动因子（城市化）对于洪水多变量风险的影响。
通过构建动态Copula 函数为核心的联合分布模型，采用LR-AICc 结合的趋势性检测手段对淮河流域蚌埠水文站的洪峰-洪量序列进行非平稳条件下的多变量频率分析，主要获得如下结论。
（1）由于城市化、气候变化因素的综合影响下，蚌埠水文站的洪峰、洪量呈现了一定的趋势型非平稳特征。其中以PISA因子、Maxpre以及NAO因子为分布参数的协变量非平稳广义极值分布模型（GEV3）是洪峰序列的最优边缘分布模型，而以PISA因子、Maxpre以及NAO因子为分布参数的协变量非平稳Gamma分布模型（GA3）是洪量序列的最优边缘分布模型；
平稳Frank Copula模型为Q-V相依性结构的最优联合分布模型。
（2）由于边缘分布参数的时变特性，50 年一遇联合重现期水平下的洪峰-洪量分位数对呈现一定的时变特性。蚌埠水文站控制的子流域内，气候变化因素是影响洪水单变量、多变量风险的主导因素。
猜你喜欢洪量洪峰平稳性赣江流域洪水峰量演变规律及联合分布研究江西水利科技(2022年4期)2022-08-04基于非平稳性度量的数字印章信息匹配数学物理学报(2021年3期)2021-07-19基于递归量化分析的振动信号非平稳性评价工程与建设(2019年5期)2020-01-19辽河干流主要控制站近75年最大洪峰及洪量变化特征分析研究水利规划与设计(2018年8期)2018-09-06淡定！新城乡(2017年8期)2017-08-26解禁洪峰纺织科学研究(2017年1期)2017-05-17高重合度齿轮传动的平稳性分析及试验厦门理工学院学报(2016年1期)2016-12-01信贷资源配置与我国经济发展平稳性湘潭大学学报（哲学社会科学版）(2015年5期)2015-11-25晨地火(2014年4期)2014-03-01适用于电算的设计洪水过程线放缩方法黑龙江水利科技(2014年7期)2014-01-21
