GAMLSS模型在黄河流域蒸发皿蒸发量序列统计建模中的应用研究

张彩玲，宋昕熠，肖伟华，方 文，剌美如

（1.中国水利水电科学研究院，北京 100038；2.北京恒华伟业科技股份有限公司，北京 100011）

1 研究背景

在气候变化和人类活动影响下，水文要素序列越来越复杂，部分序列不再服从同一分布，出现突变趋势，使序列不再满足一致性要求。基于一致性假设理论下的水文频率计算结果受到质疑，为进一步研究水文频率计算方法和变化环境对水文序列的影响，需对水文序列进行非一致性分析，它对明确和减少水文频率分析计算结果中的不确定性具有重要意义［1-2］。目前，在水文序列非一致性研究方面应用较为广泛的是GAMLSS模型，国内外诸多学者对其进行了研究和应用。

GAMLSS （generalised additive models for location，scale and shape）模型是由Rigby等［3］和Stasinopoulos等［4］在2005年提出的应用于非一致性分析的（半）参数回归模型。随着应用需求的增加，出现了以R软件为平台的GAMLSS工具包［5-6］，为时变矩法在非一致性水文频率分析中的应用提供了强有力的工具，可根据水文序列特征选择对应的函数，极大方便了水文非一致性研究工作。Gabriele等［7-8］应用Gumbel、Gamma等5种两参数分布模型研究了罗马地区长期降水和气温序列的趋势和成因，并采用GAMLSS模型分析了高度城市化的Little Sugar Creek流域的最大洪峰流量的非平稳性，取得较好的效果。张冬冬等［9］采用Gumbel、Gamma、Lognormal、Weibull、和Logistics5种分布，以GAIC为GAMLSS模型拟合评价指标，研究了大渡河流域降水频率的非一致性。江聪等［10］采用GAMLSS模型对宜昌站的年平均流量和年最小月流量序列的趋势进行了研究，结果表明年最小月流量是非平稳序列。

本文利用GAMLSS模型，通过对比分析黄河流域蒸发序列一致性与非一致性的拟合结果，分析黄河流域三级区不同时间尺度上蒸发序列的变化特征，为今后黄河流域水文分析和水资源评价等工作提供理论依据。

2 研究方法

2.1 模型定义 GAMLSS模型提供了多种方式产生的不同分布，极大地扩展了分布的种类；引入了更为复杂的半参数或非参数部分及随机效应部分，不但可以建立均值与解释变量间的模型，也可以建立分布的其他参数与解释变量间的模型，从而可以研究因变量统计特征随解释变量变化而变化的非一致性特征。GAMLSS模型［11］定义gk(·)（k=1，2，…，p）作为参数向量θk与解释变量和随机效应项之间关系的单调连接函数（Monotonic link function）。本文使用本体对应连接（Identity link function）和对数对应连接（Logarithm link function）这两种连接函数，连接函数统一的表达式为［9］：

式（1）中的参数向量和随机效应参数可以通过最大化惩罚似然函数来估计。若假定随机变量服从三参数概率分布，则GAMLSS模型的表达式为：

若分布函数表现出时间的非一致性，则GAMLSS模型的表达式为：

2.2 模型拟合优度准则与残差评价 采用GAIC准则（Generalized Akaike Information Criterion）在对序列进行拟合时，自由度取值为1～3，忽略惩罚因子带来的影响。本文以全局拟合偏差GD（Global Deviance）作为GAMLSS模型评价指标，最小GD值对应的模型为最优模型。以AIC准则值和SBC准则值［5］为评价指标评价非一致性的拟合效果，AIC和SBC准则值越小，说明非一致性拟合效果越好。本文采用概率点据相关系数法［12］，检验正态标准化的残差序列是否服从标准正态分布，Filliben值越接近1，则待检验序列越接近标准正态分布。

2.3 概率分布函数类型 选择线性函数、三次样条函数、幂次分段多项式函数作为参数解释变量之间的联系函数，选择Gumbel（GU）、Gamma（GA）、Logistic（LO）、Generalized Gamma（GG）4种分布函数，通过计算每次拟合的全局拟合偏差（Global Deviance）选取最优分布。

3 数据来源与处理

本文选取90个气象站作为黄河流域的代表站点。金堤河和天然文岩渠、花园口以下干流区间2个三级区没有气象监测站，取金堤河和天然文岩渠附近3个气象站，花园口以下干流区间附近2个气象站，1963—2013年的蒸发皿蒸发数据作为2个三级区的日蒸发量。以95个气象站的20cm口径蒸发皿日蒸发量作为研究蒸发序列非一致性演变规律的基础数据。气象站点分布情况如图1所示。

本文中蒸发量的计算通过ArcGIS实现，按照不同时间尺度上三级区各站点的点蒸发数据×面积权重，取和作为相应三级区的面蒸发数据，面积权重通过泰森多边形法实现。通过Matlab将不同时间尺度逐日蒸发数据转为年、汛期和非汛期数据，站点缺测数据通过Matlab软件使用临近站点的数据进行插补。本文中关于GAMLSS模型的计算是基于R软件平台的GAMLSS程序包，在RStudio软件中完成计算。

4 黄河流域三级区蒸发序列拟合结果分析

图1 黄河流域气象站点分布及三级区蒸发趋势结果

表1 趋势检验结果

4.1 趋势性分析 黄河流域29个三级区年、汛期（6—10月）、非汛期（11—次年5月）蒸发序列变化趋势采用Kendall检验法进行检验。统计量U≥1.98，表明蒸发序列显著增加；统计量U≤-1.98，表明序列显著减少；统计量U介于两者之间，表明序列无显著变化。因篇幅有限，蒸发序列组数较多，这里Kendall检验结果如表1所示。年、汛期、非汛期蒸发序列均呈显著下降趋势的三级区有9个，主要集中在黄河流域下游；有4个三级区年、汛期、非汛期蒸发序列均呈显著增加趋势，主要集中在黄河流域中游区域；5个三级区年、汛期、非汛期蒸发序列无明显变化趋势，主要集中于黄河流域上游中部和中游下部区域。另外，龙羊峡-兰州干流区间、石嘴山-河口镇南岸、内流区3个三级区的汛期蒸发序列无明显变化趋势，但年蒸发和非汛期蒸发序列均呈明显的下降趋势。

4.2 一致性拟合结果分析 黄河流域29个水资源三级分区在不同时间尺度上的蒸发序列采用一致性模型进行拟合。将概率密度函数中的参数设为固定值，根据GU、GA、LO、GG 4种分布函数拟合黄河流域29个三级分区的蒸发序列，根据4种备用分布函数计算所得的全局拟合偏差GD值，以GD值最小为准则优选三级区年、汛期、非汛期蒸发序列的分布函数类型，结果见图2。

图2 各三级区蒸发最优空间分布

由图2可看出，不同时间尺度上各三级区的最优分布有明显的差别，但多数三级区的蒸发序列以GG分布最优。针对年蒸发序列，以GG为最优分布函数的三级区最多，LO次之，GU最少，无以GA为最优分布的三级区；汛期蒸发序列以LO为最优分布函数的三级区最多，GG次之，GU最少，无以GA为最优分布的三级区；非汛期蒸发序列以GG为最优分布函数的三级区最多，龙羊峡-兰州干流区间以LO分布为最优分布函数，大夏河与洮河蒸发序列以GU为最优分布函数，渭河宝鸡峡-咸阳、石嘴山-河口镇南岸非汛期蒸发序列的GG分布和GA分布拟合结果一致，本文选择GA为最优分布对这两个三级区的蒸发序列进行拟合。

4.3 非一致性拟合结果分析 对于黄河流域各三级区年蒸发、汛期蒸发和非汛期蒸发序列的拟合结果，选择适合序列的最优分布函数，以时间t作为分布函数位置参数的解释变量，不同时间尺度上各三级分区蒸发序列的拟合结果，如表2所示。

表2 非一致性拟合结果

由表2可看出，位置参数与时间的函数关系种类较多。其中，年蒸发序列最优表达式为线性函数的三级区最多，自由度为3的三次样条函数次之；汛期蒸发序列最优表达式为自由度为3的幂次分段多项式函数的三级区最多；非汛期蒸发序列最优表达式为自由度为3的三次样条函数的三级区最多，自由度为1的三次样条函数次之。

对残差序列进行标准化处理，计算黄河流域各三级区蒸发序列最优分布残差的前四阶中心矩值及Filliben系数值，由于篇幅原因，这里只给出Filliben计算结果，如表3所示。三级区概率点据相关系数在0.926～0.996之间，越接近于1，说明残差序列越服从标准正态分布。对于年蒸发序列非一致

性拟合结果来说，小浪底-花园口干流区间、大通河享堂以上、大夏河与洮河、兰州-下河沿、汾河、湟水等6个三级区的Filliben系数均小于0.97；汛期蒸发序列拟合结果显示，大通河享堂以上、湟水、北洛河状头以上等3个三级区的Filliben系数均小于0.97；非汛期蒸发序列的非一致性拟合结果显示各三级区的Filliben系数均大于0.97。相对来说，非汛期蒸发序列拟合结果最好，年蒸发序列较差。

表3 模型Filliben相关系数

表4 年蒸发序列非一致性模型和一致性模型结果对比

4.4 拟合结果对比分析 年、汛期和非汛期蒸发序列通过非一致性模型和一致性模型两个模型的计算，得到两组GD、AIC和SBC准则值，其拟合效果对比情况分别见表4、表5和表6。

由表4可以看出，除了渭河宝鸡峡-咸阳、清水河与苦水河、下河沿-石嘴山、石嘴山-河口镇北岸、河口镇-龙门左岸、大汶河6个三级区外，其余三级区的拟合效果均有所改善，改善幅度不同，但相对较大。除石嘴山-河口镇南岸的SBC值没有改善外，22个三级区的全局拟合偏差GD、AIC和SBC值拟合效果比不考虑统计参数非一致性特征的情况下有不同程度改善。SBC值改善幅度最大的三级区是龙门-三门峡干流区间，为58.965，改善幅度最少的三级区是龙羊峡-兰州干流区间，为0.679。以AIC值为例，AIC值减少最多的是龙门-三门峡干流区间，为66.693，渭河咸阳-潼关减少最少，为3.354。可见龙门-三门峡干流区间是改善最大的三级区，说明其年蒸发序列的非一致性特征最为明显。总体上来说，在考虑统计参数非一致性特征情况下拟合效果有一定改善的三级区占79.31%，说明黄河流域年蒸发序列存在较为明显的非一致性特征。

表5 汛期蒸发序列非一致性模型和一致性模型结果对比

表6 非汛期蒸发序列非一致性模型和一致性模型结果对比

汛期蒸发序列的拟合结果对比情况如表5所示。考虑统计参数非一致性特征情况下，河源-玛曲、玛曲-龙羊峡、泾河张家山以上、渭河咸阳-潼关、清水河与苦水河、下河沿-石嘴山、内流区、吴堡以下右岸、河口镇-龙门左岸、大汶河10个三级区的拟合效果没有改善，且有变差趋势。龙羊峡-兰州干流区间、渭河宝鸡峡以上、汾河3个三级区仅有GD值和AIC值有所改善，受自由度影响，函数曲线复杂度增加，SBC值也有所增加。兰州-下河沿等16个三级区的全局拟合偏差GD、AIC和SBC值拟合效果比不考虑统计参数非一致性特征的情况下有一定善。SBC值改善最为明显的三级区金堤河和天然文岩渠，为45.327，改善幅度最小的三级区是兰州-下河沿，为0.456。同样以AIC准则为例，金堤河和天然文岩渠的AIC值减少最多，为47.259，渭河宝鸡峡以上的AIC值减少最少，为0.667。金堤河和天然文岩渠的AIC值和SBC值改善幅度最大，表明其汛期蒸发序列的非一致性特征最为明显。综合来说，除少部分三级区外，非一致性模型的拟合效果均比一致性模型好，且改善幅度较大，其占比为65.5%，说明黄河流域大部分三级区汛期蒸发序列存在非一致性特征。

非汛期蒸发序列的拟合结果对比情况如表6所示，龙羊峡-兰州干流区间、渭河宝鸡峡-咸阳、渭河咸阳-潼关、清水河与苦水河4个三级区的拟合结果没有改善。在基于时变的非一致性模型拟合后，内流区仅有全局拟合偏差GD和AIC值有所改善，AIC值减少1.18，幅度较小，改善不明显。河源-玛曲、玛曲-龙羊峡等24个三级区的全局拟合偏差GD、AIC和SBC值拟合效果均比不考虑统计参数非一致性特征的情况下有所改善，不同三级区的改善值跨度较大。以SBC值为例，花园口以下干流区间减少最多，为43.684，渭河宝鸡峡以上减少最少，仅为0.109。以AIC值为例，AIC值减少最多的三级区是花园口以下干流区间，为47.548，减少最少的三级区是内流区，为1.18。花园口以下干流区间的AIC值和SBC值改善幅度最大，表明其非汛期蒸发序列的非一致性特征最为明显。综合来说，86.2%的三级区在考虑统计参数非一致性特征情况下拟合效果有一定改善，说明黄河流域非汛期蒸发序列存在较为明显的非一致性特征。

5 结论

（1）黄河流域各三级区蒸发序列趋势检验结果显示：大通河享堂以上、汾河、兰州-下河沿、清水河与苦水河、河口镇-龙门左岸5个水资源三级区的年、汛期、非汛期的蒸发序列无显著变化趋势；大夏河与洮河、渭河宝鸡峡-咸阳、北洛河状头以上、吴堡以下右岸4个水资源三级区的年、汛期、非汛期的蒸发序列呈显著增加趋势；湟水、龙门-三门峡干流区间、石嘴山-河口镇北岸、伊洛河、沁丹河、三门峡-小浪底区间、小浪底-花园口干流区间、金堤河和天然文岩渠、花园口以下干流区间9个水资源三级区的年、汛期、非汛期的蒸发序列呈显著减少趋势。

（2）黄河流域各三级区蒸发序列一致性拟合结果显示：年蒸发序列最优分布为3种分布函数，以GG为最优分布的三级区最多，LO分布次之，GU分布最少；汛期蒸发序列最优分布为3种分布函数，以LO为最优分布的三级区最多，GG次之，GU分布最少；非汛期蒸发序列最优分布为4种，以GG为最优分布的三级区最多，渭河宝鸡峡-咸阳、石嘴山-河口镇南岸两个三级区蒸发序列的GG分布和GA分布拟合结果一致，LO分布和GU分布最少。

（3）黄河流域各三级区蒸发序列非一致性拟合结果显示：年蒸发序列最优表达式为线性函数最多；汛期以自由度为3的幂次分段多项式函数最多；非汛期以自由度为3的三次样条函数最多。另外各三级区的Filliben系数均接近于1，说明不同时间尺度蒸发序列拟合效果较好，以非汛期蒸发序列拟合效果最好。

（4）黄河流域各三级区蒸发序列一致性与非一致性拟合结果对比分析：年蒸发序列拟合效果较好的三级区有23个，其中22个水资源三级区的三个指标值均有改善；汛期蒸发序列拟合效果较好三级区为19个，其中16个三级区的三个拟合指标值均有改善；非汛期蒸发序列有25个三级区拟合效果较好，其中24个三级区的三个指标值均有所改善。综合来说，非汛期蒸发序列拟合效果相对最好，年、汛期和非汛期的蒸发序列非一致性拟合结果表明大部分三级区具有明显的非一致性特征。