1. 引言
气候变化和人类活动的耦合扰动下，全球很多地区的河川径流数据都表现出了显著的非平稳性(或称非一致性)，已经引起了国际水文、气象、环境等领域专家学者的广泛关注 [1]。随着近年来我国社会经济发展进程的持续推进，人类社会对水资源需求的日益增加，加之脆弱的生态环境和不合理的开发方式，我国正面临着严重的水资源短缺危机及生态环境安全问题，特别是非平稳水文变化的出现，使得流域或区域的水文情势变化更趋复杂，水安全问题凸显 [2]，亟待发展非平稳条件的水文分析方法，以提供变化环境下科学的水资源开发与生态保护。
所谓的河川径流非一致性，即水文序列的概率分布或期望、方差等统计参数时变，导致很难用一种模型实现对观测数据的拟合 [3]。为此，国内外学者对此进行了系统的研究，并提出了许多不同层面的研究成果，如极端降水的非一致性研究 [4]、降水与径流的非一致性研究、全球变暖背景下水文非一致性研究 [5]、洪水非一致性研究 [6] 等。而在方法上，也有很多具有说服力的成果，如Kharin和Zwiers [4] 使GEV分布的位置和尺度参数随时间发生变化，对非一致性条件下未来全球气温和降水极值的变化进行了研究。刘丙军等 [7] 运用时变矩法对西北江三角洲地区非一致性洪水频率进行了分析，当水文序列分布的统计参数随时间变化时，能较好反映变化环境下水文要素的非一致性变化特征。叶长青等 [8] 采用时变矩法对年最大日流量序列进行了非一致性分析，得到了较优的拟合效果。熊立华和郭生练 [9] 应用Gumbel-Hougaard Copula函数建立了长江流域水文站点的洪峰和洪量的联合分布，取得了更好的洪水频率分析效果。
基于位置、尺度和形状参数的广义加法模型(GAMLSS)是由Rigby和Stasinopoulos [10] 提出，该模型具有广泛的随机变量分布类型，包括了高偏度和高峰度的分布，并且可以对序列任何分布参数随协变量的线性变化或非线性变化进行描述。Villarini等 [6] [11] 于2015年最先将GAMLSS模型用于分析水文序列的非一致性，之后GAMLSS模型被广泛应用于水文领域非一致性研究中序列的趋势性变化。李婧等 [12] 采用GAMLSS模型对渭河流域和珠江流域5个站点的年最大洪水序列进行非一致性频率分析，有较好的拟合效果。鲁帆等 [13] 采用GAMLSS模型对黄河干流6个水文站的年径流系列进行拟合，得到各站年天然径流对应的设计径流量有随时间减少的趋势。高洁 [14] 以美国Little Sugar Creek为研究案例，利用GAMLSS模型对其83年的洪峰流量进行了分析，得到该流域2006年10年一遇的洪峰流量相当于1999年20年重现期的洪峰，并和1989年100年一遇洪峰量级接近。综上，大量的研究结果已表明，GAMLSS模型可以灵活地描述随机变量分布的任何参数与解释变量之间的线性或非线性关系，在非平稳水文频率计算方面具有有效性，且已被众多学者认可和接受。
因此，本文拟选取GAMLSS模型进行非一致性水文频率模型构建，并对马莲河流域不同水文站点的观测数据进行水文频率计算，得到不同频率下的径流设计值。通过与传统方法的对比，指出在变化环境下水文数据筛选和方法应用的建议。
2. 研究区概况与数据来源
马莲河流域地处泾河流域东北部，界于东经106˚37'~108˚34'，北纬35˚18'~37˚23'，流域主要包括了甘肃省庆阳市除镇远县外的其他七个县区，涉及了黄土高原沟壑区、黄土丘陵沟壑区和黄土丘陵区三种地貌类型，总面积19,086 km2。马莲河上游主源西川，自发源地麻黄山南流到洪德镇，与东北方向来的发源于白玉山西麓的东川汇合。西川全长82 km，流域面积4700 km2，多年平均年径流量0.378亿m3。东川全长52 km，流域面积2093 km2，多年平均年径流量0.079亿m3。两川汇合后称为环江。环江过洪德镇东南流，右岸依次纳入玄城沟、马坊川、城西川等主要支流后到达环县县城。环江继续向东南流，左岸纳入安山川，右岸纳入合道川，继而到达曲子镇，即出环县县境进入庆阳县境。干流继续东南向流动，分别纳入马岭东沟、柳树河沟、蔡家沟等小支流到达庆阳县城(庆城镇)，左岸纳入柔远河，自此以下环江改称马莲河。马莲河过庆阳县继续东南流，左纳合水川后呈正南流向，后左岸纳入固城川和九龙河，最后于甘陕边界处的政平村汇入泾河，整个马莲河水系图如图1所示。马莲河全长374.8 km，多年平均流量13.9 m 3/s，年输沙量1.34亿吨。
Figure 1. Malian River drainage and hydrological site distribution map
图1. 马莲河水系及水文站点分布图
本文选取马莲河干流洪德站、庆阳站、雨落坪站的实测径流数据(序列长度为1964~1997、2006~2012和2015~2019年)和支流柔远河贾桥站的实测径流数据(序列长度为1979~1997、2006~2012和2015~2019年)，所有径流数据均来源于黄河水利委员会印制的水文数据年鉴。图2显示了不同时期四个水文站实测径流的箱图，由图可知，马莲河径流量在过去的50余年中，已表现出了较为明显的减少趋势。
Figure 2. Map of observed runoff box at different stations of Malian River
图2. 马莲河不同站点实测径流箱图
3. 研究方法
GAMLSS模型是一种用于非一致性分析的(半)参数回归模型，可考量广泛的随机变量频率分布类型 [3]。GAMLSS模型是对GLM (Generalized Linear Models)广义线性模型和GAM (Generalized Additive Models)广义加法模型的进一步优化与修正。模型的建立不仅只针对于均值的分布，也考虑了标准差、峰度、偏度等统计量的分布情况，对非一致性序列具有良好的适用性。
本文将不同水文站点的月径流量作为研究对象，构建GAMLSS模型，使用poly多项式模型对各月径流量进行拟合，并藉由评价准则，确定最佳多项式次数和拟合分布类型，最后依据最优拟合模型，绘制分位数图，得到不同频率下的月径流量。
评价准则采用AIC和SBC准则。AIC和SBC值越小，表明模型拟合的越好，因此以最小AIC和SBC值对应的模型为最优模型 [15]。
4. 结果与讨论
应用GAMLSS模型对马莲河流域洪德、贾桥、庆阳、雨落坪站的实测水文数据进行概率分布拟合，并根据拟合分布类型的评价结果，均选择正态分布LOGNO，多项式次数选择3次，可得到不同站点的不同月份不同频率下的径流量设计值，如图3所示。
将其与P-III曲线所计算得到的结果进行对比，如图4所示。由图可知，不考虑径流序列本身的非一致性，直接应用P-III曲线进行拟合适线，所得到的不同频率的设计径流量与面向水文非一致性的GAMLSS方法差异较大。总体上表现为，特枯水年型(95%和90%)，GAMLSS模型的计算结果小于P-III曲线，而在其他年型，则表现为GAMLSS模型的计算结果大于P-III曲线。
Figure 3. Design values of annual runoff distribution measured at different hydrological stations
图3. 不同水文站实测径流年内分布设计值
Figure 4. Deviation of the calculated results between GAMLSS model and P-III curve
图4. GAMLSS模型与P-III曲线计算结果的偏差
将不同频率下两种方法的计算结果的最大偏差进行统计分析(图5)，可知特枯水期(95%, 90%)计算的最大偏差主要集中在−67%~−90%，即GAMLSS模型的计算结果小于P-III曲线，枯水期(75%)计算的最大偏差主要集中在70%~136%，GAMLSS模型的计算结果大于P-III曲线，平水期(50%)计算的最大偏差主要集中在86%~210%，丰水期(25%)计算的最大偏差主要集中在88%~223%，特丰水期(10%, 5%)计算的最大偏差主要集中在92%~296%。
Figure 5. Maximum deviation of the calculated results between GAMLSS model and P-III curve
图5. GAMLSS模型与P-III曲线计算结果的最大偏差
从年内时程分布上看，两种模型的计算偏差具有较强的聚集特征，主要集中于4~8月份，95%和90%年型，主要偏差集中于7~8月，而在75%、50%、25%、10%和5%年型下，这主要集中于4~6月份。
透过以上分析，不难发现对于非一致性径流序列，若不使用正确的频率计算方法，所得的径流设计值可能产生极大的偏差，继而影响工程规模的确定。从影响的程度看，丰水期 > 平水期 > 枯水期。
5. 结论
本文使用GAMLSS模型，对马莲河流域不同水文站点的观测数据进行水文频率计算，得到如下结论：
1) 通过对马莲河流域水文站点观测数据进行箱图分析，得知该流域50余年径流量呈现明显下降趋势。
2) 应用GAMLSS模型对马莲河流域洪德、贾桥、庆阳、雨落坪站的实测水文数据进行概率分布拟合，得到不同站点的不同月份不同频率下的径流量。
3) 通过对比GAMLSS模型设计值与P-III曲线设计值可知：不考虑径流序列本身的非一致性，直接应用P-III曲线进行拟合适线，所得到的不同频率的设计径流量与面向水文非一致性的GAMLSS方法差异较大。总体上表现为，特枯水年型(95%和90%)，GAMLSS模型的计算结果小于P-III曲线，而在其他年型，则表现为GAMLSS模型的计算结果大于P-III曲线。
4) 根据两种方法计算结果最大偏差的统计分析结果可知，特枯水期(95%, 90%)计算的最大偏差主要集中在−67%~−90%，即GAMLSS模型的计算结果小于P-III曲线，枯水期(75%)计算的最大偏差主要集中在70%~136%，GAMLSS模型的计算结果大于P-III曲线，平水期(50%)计算的最大偏差主要集中在86%~210%，丰水期(25%)计算的最大偏差主要集中在88%~223%，特丰水期(10%, 5%)计算的最大偏差主要集中在92%~296%。
5) 从年内时程分布上看，两种模型的计算偏差具有较强的聚集特征，主要集中于4~8月份。其中，95%和90%年型，主要偏差集中于7~8月，而在75%、50%、25%、10%和5%年型下，主要偏差集中于4~6月份。
综上所述，对于非一致性径流序列，如何选择正确的频率计算方法是至关重要的。频率分析计算所得设计值直接影响后续水利工程规模的确定。有效的水文分析是工程设计和流域区域水资源管理的根本。
基金项目
本研究获得国家自然科学基金项目(51979005)及陕西省自然科学基础研究计划项目(2020JM-250)支持。