高 洁

(水电水利规划设计总院，北京100120)

1 研究背景

水文系列的平稳性是工程设计中频率计算、抽样统计的基础。在工程实践中尤其重视水文系列的三性检验，要求满足可靠性、一致性和代表性。当气候变化、人类活动等多种因素影响系列的一致性时，通常进行还原计算，以保证采用平稳序列进行频率适线。但是，随着人类活动加剧，众多江河水文情势发生显著变化，水文极端事件频发[1]；区域性的非平稳序列广泛存在，非一致性问题并非个案，流域水文律情的“一致性”背景已不复存在。现有的基于平稳序列的工程水文分析计算方法将面临变化环境带来的设计频率失真风险[2]。在普遍非平稳的环境中，需要采用针对非平稳序列的统计模型和计算流程[2]。

针对非平稳序列统计参数的时变特征，Strupczewski等[3- 5]通过对时变一阶矩、二阶矩采用极大似然法、加权最小二乘法进行参数估计；并基于AIC准则优选模型的原理和方法，以波兰河流洪峰流量系列为例，在Mann-Kendall趋势检验的基础上，分析序列平稳性及统计参数变化趋势。叶长青等[6]采用时变矩模型，通过5种分布、8种趋势选择最优拟合模式，分析武江坪石水文站1964年～2008年和东江龙川水文站1953年～2008年历年最大日流量系列。研究表明，两者均符合对数正态分布，分别满足均值均方差线性变化、均值均方差抛物线型变化模式。基于平稳序列坪石站、龙川站100年一遇洪峰流量分别为4 170、6 020 m3/s；采用时变矩模型后，坪石站4 170 m3/s的洪水重现期从1970年前大于200年一遇变为2000年以后小于60年一遇。受水库调蓄影响，龙川站6 020 m3/s洪水重现期从1962年前小于100年一遇变为2000年以后大于400年一遇，研究认为传统平稳序列的洪水重现期概念应修正。刘丙军等[7]通过时变矩模型分析西江和北江水系的马口水文站和三水水文站1960年～2009年历年最大日流量系列，发现马口站符合GLO分布，均值抛物线型变化，均方差与变差系数Cv相关；三水站符合GEV分布，均值抛物线型变化，均方差平稳不变。马口站100年一遇洪水洪峰洪量，1960年为56 000 m3/s，1975年为52 000 m3/s，2009年增至67 000 m3/s；三水站100年一遇洪水洪峰洪量1960年为15 000 m3/s，1990年为16 000 m3/s，2009年增至20 000 m3/s。同一频率设计洪水在20世纪70年代较小，21世纪后明显增大，差异化明显。

考虑到城市化进程对小流域洪水的影响，Villarini等[8]以美国北卡罗来纳州集水面积110 km2的Littele Sugar Creek为对象，通过GAMLSS(Generalized Additive Models for Location， Scale and Shape)研究发现，在83年时间序列中，传统方法计算平稳序列100年一遇洪水洪峰流量设计值为3.2 m3/s。非平稳序列分析发现，设计流量3.2 m3/s的重现期从1957年的5 000年一遇减小为2007年的8年一遇，与传统频率计算成果差异较大。

本文详细介绍基于开源的专业语言R及软件平台RStudio开发的GAMLSS[9]在非平稳序列中的应用，以美国Little Sugar Creek的经典案例为例，阐明GAMLSS在非平稳时间序列分析中的流程方法。

2 研究基础

2.1 广义可加模型GAM

广义可加模型(Generalized Additive Model，GAM)是在广义线性模型(Generalized Linear Model，GLM)和可加模型(Additive Model)的基础上，由Hastie和Tibshirani于1990年提出。模型的解释变量不限于正态分布，响应变量形式灵活，采用非参数的方法进行拟合，将存在复杂非线性关系的解释变量通过不同的函数形式，以加和的方式构成模型，可反映变量间非单调、非线性关系[10- 12]，包括了可加项(additive component)、随机项(random component)和联系两者的连接函数(link function)。即

广义线性模型

g(θ)=η=α+β1X1+…+βJkXJk

(1)

可加模型

E(Y/X1，…，XJk)=α+h(x1)+…+h(xJk)

(2)

广义可加模型

(3)

式中，η为预测值或观测值；X为解释变量；Jk为解释变量的个数；α为常数项；β为解释变量的线性参数；g(·)为连接函数;θ=E(Y/X1，…，XJk)为Y的数学期望值；hj(·)为平滑函数，包括光滑样条函数、核函数、局部回归平滑函数等各种不同形式。其灵活非参数形式，可反映解释变量的非线性效应[13]。

2.2 基于位置、尺度和形状参数的广义可加模型GAMLSS

GAMLSS模型是由Rigby和Stasinopoulos于2005年提出的(半)参数回归模型。GAMLSS模型在拟合位置、尺度和形状参数的基础上，通过可加的半参数项或非参数项、随机效应项，建立响应变量统计参数(位置、尺度、形状等)与解释变量的关系[14]。GAMLSS模型在医学、经济学[15- 16]等领域已得到广泛应用，对于超峰度、平顶峰度、高度正偏或负偏等高偏度和高峰度的随机变量系列具有较好的拟合效果[17]。江聪和熊立华[17]将GAMLSS模型应用于宜昌站1882年～2009年历年平均和年最小月流量系列趋势分析。研究结果表明，该站的年最小月流量系列非平稳，均值呈非线性变化，偏度系数呈线性变化。

2.2.1模型I——不考虑随机项

GAMLSS模型假设独立随机变量第i时刻观测值的概率密度函数为f(yi|θi)。其中，θi可通过位置参数(以均值、中位数等表征)、尺度参数(方差均方差等表征)和形状参数(偏度系数，峰度系数)反映。即

(4)

针对位置参数θ1和尺度参数θ2的两参数模型，以时间t为解释变量不考虑随机效应的全参数模型

g1(θ1)=η1=X1β1=β11+β21t+…+βI1tI1-1

(5)

g2(θ2)=η2=X2β2=β12+β22t+…+βI2tI2-1

(6)

式中，可设定解释变量矩阵Xk为时间线性相关或抛物线型相关矩阵，并通过RS法[3- 4]对β进行参数估计。

2.2.2模型II——半参数化可加项

GAMLSS还引入一些子模型，如半参数可加模型、非线性半参数模型、非线性参数模型等[9]。相同的分布函数有多种参数化方案，需要针对具体问题建模优选。设Zjk=In，In是m×m的单位矩阵，γjk=hjk=hjk(xjk)，则GAMLSS半参数可加模型

(7)

式中，Xkβk是解释变量的线性函数矩阵，hjk是解释变量xjk的未知函数，xjk(j=1，2，…，Jk)均是n维向量，hjk(xjk)体现了对解释变量模拟效果更为灵活的平滑函数。时间序列模拟如不考虑其他变量的线性可加项，仅保留时间变量的平滑函数，针对位置参数θ1和尺度参数θ2的半参数化模型可表示为

(8)

(9)

上述模型均可通过AIC(Akaike Information Criterion)或SBC(Schwartz Bayes Criterion)准则对拟合效果进行优选，选出AIC值或SBC值最小者为最优模型，最终求得非平稳序列的统计参数分布函数及模型形式

(10)

AIC=-2lnML+2k

(11)

SBC=-2lnML+ln(n)

(12)

式中，n为序列长度(观测数目)；ML为似然函数极大值；k为模型参数个数。

3 研究方法

3.1 研究步骤

通过GAMLSS模型对水文系列进行分析，可分为模型拟合、残差判别和成果分析三个步骤。

(1)在模型拟合中，一是根据模型I的特点，以水文系列时间变量或大气环流因子作为解释变量，针对水文系列的位置参数、尺度参数等选择不同的概率分布和时间变化趋势进行组合，通过AIC准则选择最优模型。二是根据模型II的特点，在平滑函数中考虑时间变量或大气环流因子，选择不同的概率分布，通过AIC准则优选位置参数、尺度参数等有效自由度，并确定模型。

(2)模型残差判别，一是通过plot函数进行总体检验，二是通过worm plot检验正态标准化残差序列是否符合标准正态分布。所有散点位于两条椭圆线之间的95%置信区间，模型选择合理，模拟结果可接受。

(3)在成果分析中，通过统计参数随时间的变化趋势，获取不同时期各种频率的流量值，以揭示同一流量不同时期测算的重现期之间的差异，反映水文频率特征值是否随时间变化及变化情况，为研究水文时间序列的非平稳性提供基础。

3.2 研究思路

(1)概率分布。根据河川年最大流量年际分布特点，本研究初选两参数的Gamma、Lognormal(对数正态)、Gumbel、Weibull和Normal(正态)分布函数作为模型优选的基础。位置参数和尺度参数分别采用均值和均方差表征。

(2)模型优选：①模型I，根据水文时间序列特点，综合考虑模型拟合效果和预测外延性，选择统计参数随时间的变化函数。本研究将均值和均方差随时间变化模式，简化为三种情形：无趋势(S)、线性变化(L)和抛物线变化(P)。根据AIC准则，对均值和均方差选择最优概率分布函数和时间变化模式。②模型II，根据水文时间序列特点，本研究选择三次样条插值函数作为平滑函数。通过AIC准则，分别选出均值、均方差拟合效果最优的有效自由度。综合考虑模型应用的预测外延性，在拟合效果基本不降低、模型残差满足要求的条件下，可对自由度进行微调和降维。

(3)残差判别。通过残差散点图、QQ图、worm plot图等验证模型拟合效果。

(4)成果分析。首先判断水文时间序列是否平稳，再基于统计参数随时间的变化特征分析频率设计值。

表1 概率分布

4 研究应用

本研究主要借鉴Villarini等[8]采用的美国Little Sugar Creek 1924年～2006年共83 a年最大洪峰流量数据，介绍GAMLSS模型进行非平稳序列分析的方法和流程。

4.1 软件基础

采用开源的专业软件平台RStudio，下载并加载gamlss，gamlss.data及gamlss.dist软件包和数据包，导入水文系列数据集。

4.2 概率分布

在gamlss.family中，Gamma、Lognormal、Gumbel、Weibull、Normal概率分布，分别对应于的GA()、LOGNO()、GU()、WEI()和NO()函数。

4.3 模型拟合

4.3.1模型I

(1)对位置参数(mu)和尺度参数(sigma)，采用不同的概率分布函数，分别按照平稳过程(S)、随时间线性变化(L)、随时间抛物线型变化(P)进行模型拟合。

(2)选择AIC值最小者为最佳模型，具体AIC值见表2。

由表2可知，采用Gamma分布和对数正态分布模拟效果均较好。Gamma分布中，位置参数抛物线型、尺度参数平稳模型最优，位置参数抛物线型、尺度参数抛物线型模型其次。经多模型比较，尺度参数与时间不相关、线性相关、抛物线型相关的模拟效果差异不大，随时间的变化规律不明显。

表2 不同概率分布模式和变化趋势模型AIC值

注：下划线标识不同分布函数的最优拟合值，斜体标识最优分布函数的最优拟合值。

4.3.2模型II

对位置参数和尺度参数，采用不同的概率分布函数，选取三次样条平滑函数，通过迭代算法，推求均值和均方差最优自由度，选择AIC值最小者为最佳模型。

模型II仍以Gamma分布的拟合效果最优(见表3)，位置参数和尺度参数的自由度分别为2.0和0.9。根据Villarini等[8]研究采用的正态分布模拟结果，其流量系列均值和均方差的自由度分别为1.71和5.02，与本文相应正态分布的结果(1.7，5.1)基本一致。Villarini等为简化模型复杂度，在不显著降低模拟效果的情况下，最终将尺度参数的自由度调整为1.7。根据AIC值初步判断，本文所选用的Gamma分布更优于Villarini采用的正态分布拟合方式。

表3 不同概率分布和自由度模型AIC值

注：下划线和斜体标识分别为最优分布函数和自由度的最优拟合值。

4.4 残差判别

通过plot函数统计残差的均值、方差、偏度系数、峰度系数、Filliben相关系数，并根据核密度估计图、QQ图等检验正态标准化残差是否服从标准正态分布。当模型残差服从标准正态分布时，均值、方差、偏度系数和峰度系数分别接近0、1、0和3。经检验，模型残差满足标准正态分布要求。

4.5 模型成果

(1)根据模拟效果较优的模型I(mu-P-Sigma-S)和模型II(Gamma(2.0，0.9))，该流域年最大流量分布均值(位置参数)随时间先减小后增大，在1950年左右出现转折。1950年之前，流量呈现缓慢减小趋势，1950年后流量分布均值增幅逐渐增大，尤其在1965年以后增幅显著增大(见图2)。结合文献[8]的背景介绍，在83年观测期间，该区域城市化进程主要分为两个时期。以1965年为界，后半期的城市化进程明显加剧。关于流量分布均方差(尺度参数)的时变特征，与时间无关、线性相关、抛物线型相关均适用，自由度变化范围较大，对模拟效果影响不显著，不是模型的敏感参数。

图1 模型I均值和均方差的时变特征

图2 模型II位置参数和尺度参数的时变特征

(2)结合水文系列实测点，分析模型模拟效果。实测水文系列的点据分散，存在一定带宽(见图3)，通过单一频率的拟合线仅能从趋势上逼近，很难反映散点的分布特征。对不同时期、不同量级的实测流量赋予时间和频率双重属性，图4和图5可反映了流量随时间的变化趋势和分布的频率特征，更全面和直观的解释实测流量数据的意义。

图3 模型拟合效果

图4 模型I流量分位数

图5 模型II流量分位数

(3)根据模型成果，分析水文系列在不同时期的频率特征值。通过centiles.pred函数可导出不同时间各种频率相应流量。根据文献[8]成果，经自由度优选模型，该流域100年一遇洪水洪峰流量，在计算之初的1924年为2.8 m3/s，20世纪50年代后期最小值为2.1 m3/s，截至2007年达到5.1 m3/s。在本研究中，根据模型I和II模拟结果，100年一遇洪水洪峰流量设计值在1924年分别为2.5 m3/s和2.7 m3/s，截至2006年分别达5.1 m3/s和4.6 m3/s。其中，模型I在1947年～1950年降至最小值2.1 m3/s，模型II在1951年～1956年最小值为1.9 m3/s。多种方法计算成果基本一致，均有效揭示了该区域水文时间序列的非平稳特征。

由于水文系列的非平稳特征，在不同时期，不同频率设计值波动较大；因此，频率值的时间属性尤为重要。

在模型I中，1940年～1958年期间，1000年一遇洪水洪峰流量约2.6 m3/s，该流量级相当于1976年的100年一遇、1987年的20年一遇、1992年的10年一遇；在模型II中，2006年10年一遇洪峰流量约3.3 m3/s，相当于1999年的20年一遇，1989年的100年一遇，1926年和1980年的1000年一遇。随着时间变化，由于下垫面条件不同，某个时期10年一遇洪水量级甚至可能超过其他时期1000年一遇洪水。GAMLSS模型揭示的相同流量重现期差异正反映了水文时间序列的非平稳现象。

5 结 语

本研究通过GAMLSS模型及RStudio开发平台研究水文时间序列的非平稳特征。基于GAM模型及GAMLSS模型原理，根据对水文时间序列模型参数化方案不同的处理方式，在前人研究的基础上，按照模型I和模型II两种方案，对水文系列进行模拟。借鉴美国Little Sugar Creek 1924年至2006年共83年历年最大洪峰流量数据的研究案例，介绍了GAMLSS模型研究水文系列非平稳性的方法，通过概率分布时变模式选择、模型优选、残差判别、成果分析的系列流程，多模型多途径对比分析模型的有效性。

研究发现，GAMLSS模型研究方法和流程可有效揭示水文系列的非平稳性特征。在气候变化和人类活动加剧的区域，水文系列的非平稳性已显著影响统计参数和频率设计值。以研究案例区域为例，受城市化影响，某个时期10年一遇洪水量级甚至可能超过其他时段1000年一遇洪水，由此造成了工程设计和环境评估的极大不确定性，应在推求频率设计值时，综合考虑设计值的频率要素和时间趋势。