基于GAMLSS模型的三峡库区主汛期降雨非一致性分析

高 斌, 肖伟华, 鲁 帆, 崔 豪, 侯保灯, 张星娜, 杨 恒,3

(1.中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室, 北京 100038; 2.北京交通大学 计算机与信息技术学院, 北京 100044; 3.中国长江三峡集团有限公司 科学技术研究院, 北京 100038)

水文序列的一致性是指水文样本在过去、当前和未来均满足独立随机同分布假设。但是，很多地区由于人类活动及气候变化的影响使水文序列不再具有一致性特征[1-4]。Milly等[5]指出，变化环境下的“一致性”已经不再满足水文风险评估的假设，采用现有的一致性方法将会面临由变化环境带来的设计频率“失真”风险。

极端降水往往会引发自然灾害，例如洪水、干旱、山体滑坡和泥石流，对国家和人民群众的生命财产安全造成了严重威胁[6]。近年来，在全球范围内多地区已经监测到极端降水事件的频率和强度有所增加[7-9]。过去几十年来，人们广泛讨论了造成极端降水变化的原因。IPCC第五次会议指出，极端降水事件不可避免地可归因于自然和人为因素的综合影响[10]。具体而言，全球变暖和人为气溶胶被认为是影响极端降水变化的主要因素[11-13]。这些因素是如何干扰降水序列的一致性，还在不断讨论中。

长江流域是中国最容易受洪涝灾害的流域之一，主汛期降雨量占全年降雨总量的41%左右，高强度的降水是导致流域内主汛期洪涝灾害的主要原因[14]，且一些年份降雨的异常偏少也会引起旱灾[15]。三峡库区位于三峡大坝的上游、长江流域的腹地，对其主汛期降雨的非一致性研究，对变化环境下库区水利工程校核、水资源开发利用和长江中下游防洪安全等具有重要意义。

近年来，基于位置、尺度和形状参数的广义可加模型(GAMLSS)，在水文非一致性研究中应用较为广泛[16-19]，它可以基于R软件平台的程序包方便地对水文序列的多种分布参数(如均值、均方差和偏态系数等)进行非一致性建模。因此，本文采用GAMLSS模型，分别以时间t，4种气温形式(年平均气温Tave、年平均最高气温Tmax、年平均最低气温Tmin和年平均温差ΔT)为协变量对三峡库区的主汛期降雨序列进行非一致性研究，从时间上和空间上分别分析库区Pmfs的分布参数(均值、均方差)的非一致性特征，为三峡库区防洪减灾工作提供理论依据。

1 研究区概况及数据资料

三峡库区(28°10′—31°50′N，105°10′—110°50′E)，位于长江上游末端，库区面积约8.1万km2，是我国典型暴雨中心之一。地形特点是沟壑纵横，谷深水急，山体破碎，容易发生滑坡，生态环境十分脆弱。从地貌上看，三峡库区处于川鄂湘黔隆起褶皱带、川东平行岭谷三大构造单元和大巴山褶皱带交汇处，地形以山地(约占75%)和丘陵(约占22%)为主。库区雨量充沛，年降雨量为1 000～1 300 mm，主汛期(按6—8月)降雨量占全年降雨量的41%左右。多年平均气温17℃，气温垂直差异明显。库区由库首(包含兴山、秭归、巴东、奉节气象站)、库腹(包含镇坪、梁平、万州、利川气象站)和库尾(包含江津、合川、沙坪坝、长寿气象站)三部分组成。

本文所用数据资料来源于中国气象局国家气象信息中心提供的1961—2016年逐日降水、平均气温、最高气温和最低气温数据。根据研究需要，由原数据统计得到三峡库区1961—2016年12个站点的主汛期(6—8月)降雨量、年平均气温Tave、年平均最高气温Tmax、年平均最低气温Tmin和年平均温差ΔT作为备用数据。

2 GAMLSS模型

GAMLSS模型是由Rigby和Stasinopoulos(2005年)提出的(半)参数回归模型。与传统统计模型相比，GAMLSS模型建模更灵活[20]，主要体现在：(1) 分布的所有参数都可以建模为解释变量的函数；(2) 可以灵活地对响应变量使用多种函数分布类型，包括一系列高偏度和高峰度的离散分布和连续分布；(3) 在模型中它允许分布参数使用各种附加项。基于上述特征，使得GAMLSS模型可以有效地解决传统指数分布族无法拟合具有超峰度和平顶峰度、高度正偏或负偏的随机变量序列的难题。

2.1 模型的定义

在GAMLSS模型中，假设某一时刻t的观测值yt服从概率分布函数F(Yt|θt)，其中θ=(θt1,θt2,…,θtp)是时刻t对应的p个分布参数的向量。记任意时刻的第k个分布参数组成的向量θk=(θ1k,θ2k,…,θnk)T,(k=1,2, …,p)。若不考虑随机效应项，θk与解释变量Xk的函数关系：

gk(θk)=Xkβk

(1)

式中:Xk为解释变量矩阵,维度为n×Ik;βk=(β1k,β2k,…，βnk)是回归参数向量,长度为Ik。

在GAMLSS模型中，位置参数向量被定义为第一参数向量θ1(用均值mu表示)，尺度参数向量被定义为第二参数向量θ2(用均方差或变差系数sigma表示)，如果分布中存在其他参数，就被定义为形状参数(用nu表示)。

2.2 模型评价准则及残差评价

GAMLSS采用全局拟合偏差GD来筛选最优分布、分布参数与解释变量的最优函数类型，如下定义：

(2)

为了防止过度拟合，引入GAIC(Generalized Akaike Information Criterion)准则，即：

GAIC=GD+#df

(3)

式中:df为模型中的整体自由度；#是惩罚因子。GAIC的值越小,模型拟合效果越好。AIC准则(Akaike Information Criterion)和SBC准则(Schwarz Bayesian Criterion)是GAIC准则中最广泛的两种特例,罚惩因子#分别为2和ln(n),n为样本序列长度。

模型残差评价是衡量模型拟合好坏的一个重要标准，其可以通过残差正态QQ图、worm plot图、分位图等定性验证模型拟合效果，残差评价原理及标准参见文献[20]。

2.3 概率分布

根据1960—2016年三峡库区主汛期降雨总量系列频率分布特点，这里选用比较常见的Gumbel(GU)、Weibull(WEI)、Gamma(GA)、Logistic(LO)和Normal(NO)分布作为GAMLSS模型的备用分布函数。

3 结果与分析

3.1 Mann-Kendall趋势检验

采用Mann-Kendall趋势检验法对主汛期降雨量序列进行趋势性分析，库区12个站点有5个站点主汛期降雨量呈不显著减小趋势，7个站点主汛期降雨量呈增加趋势，见图1。从空间上看，三峡库区主汛期降雨量由库首秭归地区的显著增加(显著性水平α=0.05)过渡到库腹的不显著减少，再到库尾的不显著增加。

图1 库区多年平均主汛期降雨量空间变化趋势

3.2 基于协变量时间t的非一致性频率分析

以时间t为协变量，建立主汛期降雨量序列分布参数与时间t的关系。使用9种组合方式和5种分布函数类型(共计45种拟合方式)对分布参数进行拟合，组合方式如下：(1) 均值(mu)、均方差(sigma)均为常数；(2) 均值(mu)为常数，均方差(sigma)随时间t线性变化；(3) 均值(mu)为常数，均方差(sigma)随时间t抛物线型变化；(4) 均值(mu)随时间t线性变化，均方差(sigma)为常数；(5) 均值(mu)、均方差(sigma)均随时间t线性变化；(6) 均值(mu)随时间t线性变化，均方差(sigma)随时间t抛物线型变化；(7) 均值(mu)随时间t抛物线型变化，均方差(sigma)为常数；(8) 均值(mu)随时间t抛物线型变化，均方差(sigma)随时间t线性变化；(9) 均值(mu)、均方差(sigma)均随时间t抛物线型变化。将AIC和SBC值作为择优的标准，AIC和SBC值越小，表示模型的模拟效果越好。前文模型评价准则及残差评价中提到，AIC值为惩罚因子为2时的GAIC值，而SBC值为惩罚因子为ln(n)时的GAIC值[ln(56)=4.025，AIC值恒小于SBC值，差值为2.025·df]，因此择优标准AIC值和SBC值寻优效果一致，下文选用AIC值作为典型指标进行迭代寻优。

三峡库区12个站点主汛期降雨量GAMLSS模型的AIC值见图2。结果显示：除长寿站最优分布类型为Weibull外，其他11个站点的最优分布均为Gamma分布；江津、万州、奉节、秭归、镇坪和沙坪坝的AIC值拟合效果在考虑分布参数具有非一致性特征的情况下比认为分布参数为常数时要有所减小，6个站点AIC值依次减小了4.8，1.0，1.6，7.1，2.0，0.9，其中秭归的改善幅度最大；其他6个站点AIC值在分布参数为常数时较小，可认为这些序列仍然是一致的。

上述根据AIC值准则，选出了各站点以时间t为协变量的最优GAMLSS模型，需进一步验证模型拟合的好坏。根据分布参数的回归方程绘制各站点相应的分位图,统计结果见表1，结果表明：仅有长寿站、梁平站和利川站的实测点据频率与相应理论概率最大差值略微超过5%，其他的站点实测点据频率与相应理论概率最大差值均在5%以内，表明统计参数方程与实测点的拟合关系良好。从分位图中还可以得到，江津地区主汛期降雨量的方差不断增大。秭归、沙坪坝地区主汛期降雨量均值不断增大，其中秭归地区增长速度约为29.6 mm/10 a；万州地区主汛期降雨量方差、奉节和镇坪地区主汛期降雨量的均值在1990年以前越来越大，之后逐渐恢复。其他地区主汛期降雨量的均值和方差均为常数。

表1 以时间t为协变量的实测点据频率与理论分位曲线概率对比 %

注：横坐标1—9分别表示9种组合方式，其中1表示统计参数均为常数，为传统的一致性模型，而2—9是分布参数与时间t相关，为非一致性模型。

虽然主汛期降雨序列采用了多种分布函数以及以时间t为协变量的多种组合方式得到了较优的GAMLSS模型，但在气候变化的背景下，缺乏相应的物理意义。于是在这里进一步用与气候变化密切相关的气温要素作为协变量拟合模型，研究其对主汛期降雨量搭建的GAMLSS模型是否具有修正意义。

3.3 基于协变量气温的非一致性频率分析

在对年平均气温Tave、年平均最高气温Tmax、年平均最低气温Tmin和年平均温差ΔT作为协变量分析前，先计算了它们分别与主汛期降雨量的线性相关系数r，见图3。

图3 响应变量与各气温序列相关系数对比

结果显示：主汛期降雨量与4种气温形式均呈负相关，其中长寿、梁平、万州、镇坪、利川、合川和沙坪坝共7个站点的主汛期降雨量与年均日温差ΔT相关性很强，相关系数绝对范围为0.31～0.42；奉节、兴山和秭归的主汛期降雨量序列与年平均气温Tave相关性最强，相关系数分别为-0.22，-0.29和-0.53；江津、巴东站的类似相关系数绝对值均没有超过0.2，说明从线性趋势上来看，这两个站的主汛期降雨量序列受气温影响不明显。

与前文以时间t为协变量时的9种组合方式相同(各分布参数与气温的组合方式：常数、一次函数、抛物线关系)，分别以Tave，Tmax，Tmin和ΔT为协变量对各站点的Pmfs进行非一致性频率分析。

根据AIC值最小原则，在GAMLSS模型里择优得到最优分布函数回归方程，见表2。结果表示：除了长寿站以ΔT为协变量拟合时最优分布类型为Logistic分布外，其他组合方式的最优分布均为Gamma分布(11/12)；长寿、梁平、万州、镇坪、利川、合川和沙坪坝在以平均日温差为协变量拟合响应变量时，AIC值最小，模拟效果最优；奉节、巴东、兴山、秭归在以年平均气温为协变量拟合响应变量时，模拟效果最优，均与相关系数分析结果对应。

表2 以气温因子为协变量的GAMLSS模型最优分布参数

模型拟合的好坏需要通过残差分析进行评价。首先对模型的残差进行正态标准化处理，然后进行PPCC检验[20]。表3为残差分布统计结果，表明各模型残差序列分布情况接近标准正态分布，Fillben相关系数也非常接近1，这表明模型拟合效果十分理想。

表3 以气温因子为协变量的残差分布计算

以气温因子为协变量的回归模型分位图见图4，江津地区主汛期降雨量的均方差增大；秭归地区在1998年左右发生了均值突变，主汛期降雨量由从原来的480 mm突变至580 mm左右；其他站点的均值也出现了以气温为协变量的非一致性，但方差为常数(由于篇幅有限，其他站点未列出)。

3.4 模型拟合结果对比及不确定性分析

从GAMLSS模型运算结果中，得到位置参数μ、尺度参数σ为常数或在不同协变量下的AIC值，见图5。(1) 长寿、梁平、巴东、兴山、利川和合川的最优AIC值比较：AIC(时间t为协变量的最优分布组合)>AIC(分布参数为常数的最优分布)>AIC(气温为协变量的最优分布组合)，说明这6个站点在以时间t为协变量拟合时分布参数是固定不变的，但考虑气温因子后，其表现出分布参数随气温变化的非一致性。(2) 江津、万州、秭归、镇坪和沙坪坝最优AIC值比较：AIC(分布参数为常数的最优分布)>AIC(时间t为协变量的最优分布组合)>AIC(气温为协变量的最优分布组合)，这说明这5个站点的分布参数不仅对时间t变化，也随气温变化，但考虑以气温为协变量拟合模型不仅具有物理意义，拟合效果也最优。(3) 奉节站的最优AIC值规律：AIC(t为协变量的最优分布组合)>AIC(气温为协变量的最优分布组合)>AIC(分布参数为常数的最优分布)，奉节站相应变量的分布参数既没随时间变化也没随气温变化，可能人类活动或者其他气象因素对其的干扰掩盖了气温对主汛期降雨量的影响，有待进一步验证。

注：A代表以时间为协变量；B代表气温为协变量。

根据分布参数非一致特征的空间分布情况，可以得出：(1) 以t为协变量时，库尾江津地区主汛期降雨量的方差不断增大，主汛期降雨量的量级的可能增大；库首秭归地区、库尾沙坪坝地区主汛期降雨量均值不断增大，其中秭归地区增长速度约为29.6 mm/10 a，理论上这两地主汛期降雨量未来会越来多。(2) 以气温为协变量时，库首各站主汛期降雨与年均日温差ΔT相关系数最高(-0.42～-0.31)，并在以ΔT为协变量时模型拟合效果最好；库区中部和西南角各站点主汛期降雨与年平均日气温Tave相关系数最高(-0.53～0.22)，并在以Tave为协变量时模型拟合效果最好。库尾江津地区在以Tave为协变量时模型拟合效果最好。

本研究以气温要素为协变量对主汛期降雨量进行了非一致性分析，分位图中部分站点的经验频率与理论概率的差值仍在5%以上，说明还有其他因素或对主汛期降雨的非一致性影响较大，例如：土地利用变化因素[21]、温室气体排量因素[22]、自然气候变化因素[23](北太平洋涛动(NPO)、太平洋年代际涛动(PDO)、北极涛动(AO)、南方涛动(SO)和东部热带太平洋海表温度异常等)，这些因素是否能够改善模型，有待进一步验证。

注：1表示分布参数为常数；t表示以时间为协变量拟合分布参数。

4 结 论

(1) 从趋势检验上看，从1961—2016年，三峡库区12个气象站点中有5个站点的主汛期降雨量(Pmfs)呈不显著减小趋势，7个站点主汛期降雨量呈不显著增加趋势。从空间上看，三峡库区主汛期降雨量由库首的显著增加(显著性水平α=0.05)到库腹的不显著减少，再过渡到库尾的不显著增加。

(2) 以t为协变量拟合时，主汛期降雨量(Pmfs)的均值、均方差具有非一致性的站点个数分别为4，2。其中库尾江津地区主汛期降雨量的均方差不断增大，主汛期降雨量的量级的在未来可能增大；库首的秭归地区、库尾的沙坪坝地区主汛期降雨量均值不断增大(秭归地区增长速度约为29.6 mm/10 a)，理论上这两地主汛期降雨量逐渐增加。这些地区在未来可能引发洪涝、滑坡、干旱等自然灾害。

(3) 以气温为协变量拟合时，库区站点主汛期降雨量均呈现非一致性(均值、均方差具有非一致性特征的站点个数为10，2)，模型总体拟合效果优于前者且具有物理意义。库首以Tave为协变量时，模型拟合效果最好；库腹、库尾总体上以ΔT为协变量时，拟合效果最优。值得注意的是，库首秭归地区主汛期降雨量在1998年左右发生了均值突变，均值由从原来的480 mm突增至580 mm左右，可能是由于人类活动或气候变化导致年平均气温Tave变化引起的；库尾的江津地区主汛期降雨量的均方差呈增大趋势，这两地出现自然灾害的风险同以t为协变量的拟合结果。

(4) 3种模型(分布参数均为常数、以时间t为协变量和以气温为协变量)的比较表明，有近92%的站点最优分布函数类型为Gamma分布，且以气温为协变量的非一致模型更好地反映了观测数据的变化。这些结果证实，在主汛期降雨序列中非一致分析考虑物理协变量是必要且有效的，可以根据气温因素来减小三峡库区主汛期降雨量非一致性分析的不确定性。