GAMLSS模型在黃河流域蒸發(fā)皿蒸發(fā)量序列統(tǒng)計建模中的應用研究

張彩玲，宋昕熠，肖偉華，方 文，剌美如

（1.中國水利水電科學研究院，北京 100038；2.北京恒華偉業(yè)科技股份有限公司，北京 100011）

1 研究背景

在氣候變化和人類活動影響下，水文要素序列越來越復雜，部分序列不再服從同一分布，出現(xiàn)突變趨勢，使序列不再滿足一致性要求?；谝恢滦约僭O理論下的水文頻率計算結(jié)果受到質(zhì)疑，為進一步研究水文頻率計算方法和變化環(huán)境對水文序列的影響，需對水文序列進行非一致性分析，它對明確和減少水文頻率分析計算結(jié)果中的不確定性具有重要意義［1-2］。目前，在水文序列非一致性研究方面應用較為廣泛的是GAMLSS模型，國內(nèi)外諸多學者對其進行了研究和應用。

GAMLSS （generalised additive models for location，scale and shape）模型是由Rigby等［3］和Stasinopoulos等［4］在2005年提出的應用于非一致性分析的（半）參數(shù)回歸模型。隨著應用需求的增加，出現(xiàn)了以R軟件為平臺的GAMLSS工具包［5-6］，為時變矩法在非一致性水文頻率分析中的應用提供了強有力的工具，可根據(jù)水文序列特征選擇對應的函數(shù)，極大方便了水文非一致性研究工作。Gabriele等［7-8］應用Gumbel、Gamma等5種兩參數(shù)分布模型研究了羅馬地區(qū)長期降水和氣溫序列的趨勢和成因，并采用GAMLSS模型分析了高度城市化的Little Sugar Creek流域的最大洪峰流量的非平穩(wěn)性，取得較好的效果。張冬冬等［9］采用Gumbel、Gamma、Lognormal、Weibull、和Logistics5種分布，以GAIC為GAMLSS模型擬合評價指標，研究了大渡河流域降水頻率的非一致性。江聰?shù)龋?0］采用GAMLSS模型對宜昌站的年平均流量和年最小月流量序列的趨勢進行了研究，結(jié)果表明年最小月流量是非平穩(wěn)序列。

本文利用GAMLSS模型，通過對比分析黃河流域蒸發(fā)序列一致性與非一致性的擬合結(jié)果，分析黃河流域三級區(qū)不同時間尺度上蒸發(fā)序列的變化特征，為今后黃河流域水文分析和水資源評價等工作提供理論依據(jù)。

2 研究方法

2.1 模型定義 GAMLSS模型提供了多種方式產(chǎn)生的不同分布，極大地擴展了分布的種類；引入了更為復雜的半?yún)?shù)或非參數(shù)部分及隨機效應部分，不但可以建立均值與解釋變量間的模型，也可以建立分布的其他參數(shù)與解釋變量間的模型，從而可以研究因變量統(tǒng)計特征隨解釋變量變化而變化的非一致性特征。GAMLSS模型［11］定義gk(·)（k=1，2，…，p）作為參數(shù)向量θk與解釋變量和隨機效應項之間關系的單調(diào)連接函數(shù)（Monotonic link function）。本文使用本體對應連接（Identity link function）和對數(shù)對應連接（Logarithm link function）這兩種連接函數(shù)，連接函數(shù)統(tǒng)一的表達式為［9］：

式（1）中的參數(shù)向量和隨機效應參數(shù)可以通過最大化懲罰似然函數(shù)來估計。若假定隨機變量服從三參數(shù)概率分布，則GAMLSS模型的表達式為：

若分布函數(shù)表現(xiàn)出時間的非一致性，則GAMLSS模型的表達式為：

2.2 模型擬合優(yōu)度準則與殘差評價 采用GAIC準則（Generalized Akaike Information Criterion）在對序列進行擬合時，自由度取值為1～3，忽略懲罰因子帶來的影響。本文以全局擬合偏差GD（Global Deviance）作為GAMLSS模型評價指標，最小GD值對應的模型為最優(yōu)模型。以AIC準則值和SBC準則值［5］為評價指標評價非一致性的擬合效果，AIC和SBC準則值越小，說明非一致性擬合效果越好。本文采用概率點據(jù)相關系數(shù)法［12］，檢驗正態(tài)標準化的殘差序列是否服從標準正態(tài)分布，F(xiàn)illiben值越接近1，則待檢驗序列越接近標準正態(tài)分布。

2.3 概率分布函數(shù)類型 選擇線性函數(shù)、三次樣條函數(shù)、冪次分段多項式函數(shù)作為參數(shù)解釋變量之間的聯(lián)系函數(shù)，選擇Gumbel（GU）、Gamma（GA）、Logistic（LO）、Generalized Gamma（GG）4種分布函數(shù)，通過計算每次擬合的全局擬合偏差（Global Deviance）選取最優(yōu)分布。

3 數(shù)據(jù)來源與處理

本文選取90個氣象站作為黃河流域的代表站點。金堤河和天然文巖渠、花園口以下干流區(qū)間2個三級區(qū)沒有氣象監(jiān)測站，取金堤河和天然文巖渠附近3個氣象站，花園口以下干流區(qū)間附近2個氣象站，1963—2013年的蒸發(fā)皿蒸發(fā)數(shù)據(jù)作為2個三級區(qū)的日蒸發(fā)量。以95個氣象站的20cm口徑蒸發(fā)皿日蒸發(fā)量作為研究蒸發(fā)序列非一致性演變規(guī)律的基礎數(shù)據(jù)。氣象站點分布情況如圖1所示。

本文中蒸發(fā)量的計算通過ArcGIS實現(xiàn)，按照不同時間尺度上三級區(qū)各站點的點蒸發(fā)數(shù)據(jù)×面積權重，取和作為相應三級區(qū)的面蒸發(fā)數(shù)據(jù)，面積權重通過泰森多邊形法實現(xiàn)。通過Matlab將不同時間尺度逐日蒸發(fā)數(shù)據(jù)轉(zhuǎn)為年、汛期和非汛期數(shù)據(jù)，站點缺測數(shù)據(jù)通過Matlab軟件使用臨近站點的數(shù)據(jù)進行插補。本文中關于GAMLSS模型的計算是基于R軟件平臺的GAMLSS程序包，在RStudio軟件中完成計算。

4 黃河流域三級區(qū)蒸發(fā)序列擬合結(jié)果分析

圖1 黃河流域氣象站點分布及三級區(qū)蒸發(fā)趨勢結(jié)果

表1 趨勢檢驗結(jié)果

4.1 趨勢性分析 黃河流域29個三級區(qū)年、汛期（6—10月）、非汛期（11—次年5月）蒸發(fā)序列變化趨勢采用Kendall檢驗法進行檢驗。統(tǒng)計量U≥1.98，表明蒸發(fā)序列顯著增加；統(tǒng)計量U≤-1.98，表明序列顯著減少；統(tǒng)計量U介于兩者之間，表明序列無顯著變化。因篇幅有限，蒸發(fā)序列組數(shù)較多，這里Kendall檢驗結(jié)果如表1所示。年、汛期、非汛期蒸發(fā)序列均呈顯著下降趨勢的三級區(qū)有9個，主要集中在黃河流域下游；有4個三級區(qū)年、汛期、非汛期蒸發(fā)序列均呈顯著增加趨勢，主要集中在黃河流域中游區(qū)域；5個三級區(qū)年、汛期、非汛期蒸發(fā)序列無明顯變化趨勢，主要集中于黃河流域上游中部和中游下部區(qū)域。另外，龍羊峽-蘭州干流區(qū)間、石嘴山-河口鎮(zhèn)南岸、內(nèi)流區(qū)3個三級區(qū)的汛期蒸發(fā)序列無明顯變化趨勢，但年蒸發(fā)和非汛期蒸發(fā)序列均呈明顯的下降趨勢。

4.2 一致性擬合結(jié)果分析 黃河流域29個水資源三級分區(qū)在不同時間尺度上的蒸發(fā)序列采用一致性模型進行擬合。將概率密度函數(shù)中的參數(shù)設為固定值，根據(jù)GU、GA、LO、GG 4種分布函數(shù)擬合黃河流域29個三級分區(qū)的蒸發(fā)序列，根據(jù)4種備用分布函數(shù)計算所得的全局擬合偏差GD值，以GD值最小為準則優(yōu)選三級區(qū)年、汛期、非汛期蒸發(fā)序列的分布函數(shù)類型，結(jié)果見圖2。

圖2 各三級區(qū)蒸發(fā)最優(yōu)空間分布

由圖2可看出，不同時間尺度上各三級區(qū)的最優(yōu)分布有明顯的差別，但多數(shù)三級區(qū)的蒸發(fā)序列以GG分布最優(yōu)。針對年蒸發(fā)序列，以GG為最優(yōu)分布函數(shù)的三級區(qū)最多，LO次之，GU最少，無以GA為最優(yōu)分布的三級區(qū)；汛期蒸發(fā)序列以LO為最優(yōu)分布函數(shù)的三級區(qū)最多，GG次之，GU最少，無以GA為最優(yōu)分布的三級區(qū)；非汛期蒸發(fā)序列以GG為最優(yōu)分布函數(shù)的三級區(qū)最多，龍羊峽-蘭州干流區(qū)間以LO分布為最優(yōu)分布函數(shù)，大夏河與洮河蒸發(fā)序列以GU為最優(yōu)分布函數(shù)，渭河寶雞峽-咸陽、石嘴山-河口鎮(zhèn)南岸非汛期蒸發(fā)序列的GG分布和GA分布擬合結(jié)果一致，本文選擇GA為最優(yōu)分布對這兩個三級區(qū)的蒸發(fā)序列進行擬合。

4.3 非一致性擬合結(jié)果分析 對于黃河流域各三級區(qū)年蒸發(fā)、汛期蒸發(fā)和非汛期蒸發(fā)序列的擬合結(jié)果，選擇適合序列的最優(yōu)分布函數(shù)，以時間t作為分布函數(shù)位置參數(shù)的解釋變量，不同時間尺度上各三級分區(qū)蒸發(fā)序列的擬合結(jié)果，如表2所示。

表2 非一致性擬合結(jié)果

由表2可看出，位置參數(shù)與時間的函數(shù)關系種類較多。其中，年蒸發(fā)序列最優(yōu)表達式為線性函數(shù)的三級區(qū)最多，自由度為3的三次樣條函數(shù)次之；汛期蒸發(fā)序列最優(yōu)表達式為自由度為3的冪次分段多項式函數(shù)的三級區(qū)最多；非汛期蒸發(fā)序列最優(yōu)表達式為自由度為3的三次樣條函數(shù)的三級區(qū)最多，自由度為1的三次樣條函數(shù)次之。

對殘差序列進行標準化處理，計算黃河流域各三級區(qū)蒸發(fā)序列最優(yōu)分布殘差的前四階中心矩值及Filliben系數(shù)值，由于篇幅原因，這里只給出Filliben計算結(jié)果，如表3所示。三級區(qū)概率點據(jù)相關系數(shù)在0.926～0.996之間，越接近于1，說明殘差序列越服從標準正態(tài)分布。對于年蒸發(fā)序列非一致

性擬合結(jié)果來說，小浪底-花園口干流區(qū)間、大通河享堂以上、大夏河與洮河、蘭州-下河沿、汾河、湟水等6個三級區(qū)的Filliben系數(shù)均小于0.97；汛期蒸發(fā)序列擬合結(jié)果顯示，大通河享堂以上、湟水、北洛河狀頭以上等3個三級區(qū)的Filliben系數(shù)均小于0.97；非汛期蒸發(fā)序列的非一致性擬合結(jié)果顯示各三級區(qū)的Filliben系數(shù)均大于0.97。相對來說，非汛期蒸發(fā)序列擬合結(jié)果最好，年蒸發(fā)序列較差。

表3 模型Filliben相關系數(shù)

表4 年蒸發(fā)序列非一致性模型和一致性模型結(jié)果對比

4.4 擬合結(jié)果對比分析 年、汛期和非汛期蒸發(fā)序列通過非一致性模型和一致性模型兩個模型的計算，得到兩組GD、AIC和SBC準則值，其擬合效果對比情況分別見表4、表5和表6。

由表4可以看出，除了渭河寶雞峽-咸陽、清水河與苦水河、下河沿-石嘴山、石嘴山-河口鎮(zhèn)北岸、河口鎮(zhèn)-龍門左岸、大汶河6個三級區(qū)外，其余三級區(qū)的擬合效果均有所改善，改善幅度不同，但相對較大。除石嘴山-河口鎮(zhèn)南岸的SBC值沒有改善外，22個三級區(qū)的全局擬合偏差GD、AIC和SBC值擬合效果比不考慮統(tǒng)計參數(shù)非一致性特征的情況下有不同程度改善。SBC值改善幅度最大的三級區(qū)是龍門-三門峽干流區(qū)間，為58.965，改善幅度最少的三級區(qū)是龍羊峽-蘭州干流區(qū)間，為0.679。以AIC值為例，AIC值減少最多的是龍門-三門峽干流區(qū)間，為66.693，渭河咸陽-潼關減少最少，為3.354。可見龍門-三門峽干流區(qū)間是改善最大的三級區(qū)，說明其年蒸發(fā)序列的非一致性特征最為明顯?？傮w上來說，在考慮統(tǒng)計參數(shù)非一致性特征情況下擬合效果有一定改善的三級區(qū)占79.31%，說明黃河流域年蒸發(fā)序列存在較為明顯的非一致性特征。

表5 汛期蒸發(fā)序列非一致性模型和一致性模型結(jié)果對比

表6 非汛期蒸發(fā)序列非一致性模型和一致性模型結(jié)果對比

汛期蒸發(fā)序列的擬合結(jié)果對比情況如表5所示。考慮統(tǒng)計參數(shù)非一致性特征情況下，河源-瑪曲、瑪曲-龍羊峽、涇河張家山以上、渭河咸陽-潼關、清水河與苦水河、下河沿-石嘴山、內(nèi)流區(qū)、吳堡以下右岸、河口鎮(zhèn)-龍門左岸、大汶河10個三級區(qū)的擬合效果沒有改善，且有變差趨勢。龍羊峽-蘭州干流區(qū)間、渭河寶雞峽以上、汾河3個三級區(qū)僅有GD值和AIC值有所改善，受自由度影響，函數(shù)曲線復雜度增加，SBC值也有所增加。蘭州-下河沿等16個三級區(qū)的全局擬合偏差GD、AIC和SBC值擬合效果比不考慮統(tǒng)計參數(shù)非一致性特征的情況下有一定善。SBC值改善最為明顯的三級區(qū)金堤河和天然文巖渠，為45.327，改善幅度最小的三級區(qū)是蘭州-下河沿，為0.456。同樣以AIC準則為例，金堤河和天然文巖渠的AIC值減少最多，為47.259，渭河寶雞峽以上的AIC值減少最少，為0.667。金堤河和天然文巖渠的AIC值和SBC值改善幅度最大，表明其汛期蒸發(fā)序列的非一致性特征最為明顯。綜合來說，除少部分三級區(qū)外，非一致性模型的擬合效果均比一致性模型好，且改善幅度較大，其占比為65.5%，說明黃河流域大部分三級區(qū)汛期蒸發(fā)序列存在非一致性特征。

非汛期蒸發(fā)序列的擬合結(jié)果對比情況如表6所示，龍羊峽-蘭州干流區(qū)間、渭河寶雞峽-咸陽、渭河咸陽-潼關、清水河與苦水河4個三級區(qū)的擬合結(jié)果沒有改善。在基于時變的非一致性模型擬合后，內(nèi)流區(qū)僅有全局擬合偏差GD和AIC值有所改善，AIC值減少1.18，幅度較小，改善不明顯。河源-瑪曲、瑪曲-龍羊峽等24個三級區(qū)的全局擬合偏差GD、AIC和SBC值擬合效果均比不考慮統(tǒng)計參數(shù)非一致性特征的情況下有所改善，不同三級區(qū)的改善值跨度較大。以SBC值為例，花園口以下干流區(qū)間減少最多，為43.684，渭河寶雞峽以上減少最少，僅為0.109。以AIC值為例，AIC值減少最多的三級區(qū)是花園口以下干流區(qū)間，為47.548，減少最少的三級區(qū)是內(nèi)流區(qū)，為1.18?；▓@口以下干流區(qū)間的AIC值和SBC值改善幅度最大，表明其非汛期蒸發(fā)序列的非一致性特征最為明顯。綜合來說，86.2%的三級區(qū)在考慮統(tǒng)計參數(shù)非一致性特征情況下擬合效果有一定改善，說明黃河流域非汛期蒸發(fā)序列存在較為明顯的非一致性特征。

5 結(jié)論

（1）黃河流域各三級區(qū)蒸發(fā)序列趨勢檢驗結(jié)果顯示：大通河享堂以上、汾河、蘭州-下河沿、清水河與苦水河、河口鎮(zhèn)-龍門左岸5個水資源三級區(qū)的年、汛期、非汛期的蒸發(fā)序列無顯著變化趨勢；大夏河與洮河、渭河寶雞峽-咸陽、北洛河狀頭以上、吳堡以下右岸4個水資源三級區(qū)的年、汛期、非汛期的蒸發(fā)序列呈顯著增加趨勢；湟水、龍門-三門峽干流區(qū)間、石嘴山-河口鎮(zhèn)北岸、伊洛河、沁丹河、三門峽-小浪底區(qū)間、小浪底-花園口干流區(qū)間、金堤河和天然文巖渠、花園口以下干流區(qū)間9個水資源三級區(qū)的年、汛期、非汛期的蒸發(fā)序列呈顯著減少趨勢。

（2）黃河流域各三級區(qū)蒸發(fā)序列一致性擬合結(jié)果顯示：年蒸發(fā)序列最優(yōu)分布為3種分布函數(shù)，以GG為最優(yōu)分布的三級區(qū)最多，LO分布次之，GU分布最少；汛期蒸發(fā)序列最優(yōu)分布為3種分布函數(shù)，以LO為最優(yōu)分布的三級區(qū)最多，GG次之，GU分布最少；非汛期蒸發(fā)序列最優(yōu)分布為4種，以GG為最優(yōu)分布的三級區(qū)最多，渭河寶雞峽-咸陽、石嘴山-河口鎮(zhèn)南岸兩個三級區(qū)蒸發(fā)序列的GG分布和GA分布擬合結(jié)果一致，LO分布和GU分布最少。

（3）黃河流域各三級區(qū)蒸發(fā)序列非一致性擬合結(jié)果顯示：年蒸發(fā)序列最優(yōu)表達式為線性函數(shù)最多；汛期以自由度為3的冪次分段多項式函數(shù)最多；非汛期以自由度為3的三次樣條函數(shù)最多。另外各三級區(qū)的Filliben系數(shù)均接近于1，說明不同時間尺度蒸發(fā)序列擬合效果較好，以非汛期蒸發(fā)序列擬合效果最好。

（4）黃河流域各三級區(qū)蒸發(fā)序列一致性與非一致性擬合結(jié)果對比分析：年蒸發(fā)序列擬合效果較好的三級區(qū)有23個，其中22個水資源三級區(qū)的三個指標值均有改善；汛期蒸發(fā)序列擬合效果較好三級區(qū)為19個，其中16個三級區(qū)的三個擬合指標值均有改善；非汛期蒸發(fā)序列有25個三級區(qū)擬合效果較好，其中24個三級區(qū)的三個指標值均有所改善。綜合來說，非汛期蒸發(fā)序列擬合效果相對最好，年、汛期和非汛期的蒸發(fā)序列非一致性擬合結(jié)果表明大部分三級區(qū)具有明顯的非一致性特征。