基于GAMLSS模型的三峽庫區(qū)主汛期降雨非一致性分析

高 斌, 肖偉華, 魯 帆, 崔 豪, 侯保燈, 張星娜, 楊 恒,3

(1.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點實驗室, 北京 100038; 2.北京交通大學(xué) 計算機與信息技術(shù)學(xué)院, 北京 100044; 3.中國長江三峽集團有限公司 科學(xué)技術(shù)研究院, 北京 100038)

水文序列的一致性是指水文樣本在過去、當前和未來均滿足獨立隨機同分布假設(shè)。但是，很多地區(qū)由于人類活動及氣候變化的影響使水文序列不再具有一致性特征[1-4]。Milly等[5]指出，變化環(huán)境下的“一致性”已經(jīng)不再滿足水文風險評估的假設(shè)，采用現(xiàn)有的一致性方法將會面臨由變化環(huán)境帶來的設(shè)計頻率“失真”風險。

極端降水往往會引發(fā)自然災(zāi)害，例如洪水、干旱、山體滑坡和泥石流，對國家和人民群眾的生命財產(chǎn)安全造成了嚴重威脅[6]。近年來，在全球范圍內(nèi)多地區(qū)已經(jīng)監(jiān)測到極端降水事件的頻率和強度有所增加[7-9]。過去幾十年來，人們廣泛討論了造成極端降水變化的原因。IPCC第五次會議指出，極端降水事件不可避免地可歸因于自然和人為因素的綜合影響[10]。具體而言，全球變暖和人為氣溶膠被認為是影響極端降水變化的主要因素[11-13]。這些因素是如何干擾降水序列的一致性，還在不斷討論中。

長江流域是中國最容易受洪澇災(zāi)害的流域之一，主汛期降雨量占全年降雨總量的41%左右，高強度的降水是導(dǎo)致流域內(nèi)主汛期洪澇災(zāi)害的主要原因[14]，且一些年份降雨的異常偏少也會引起旱災(zāi)[15]。三峽庫區(qū)位于三峽大壩的上游、長江流域的腹地，對其主汛期降雨的非一致性研究，對變化環(huán)境下庫區(qū)水利工程校核、水資源開發(fā)利用和長江中下游防洪安全等具有重要意義。

近年來，基于位置、尺度和形狀參數(shù)的廣義可加模型(GAMLSS)，在水文非一致性研究中應(yīng)用較為廣泛[16-19]，它可以基于R軟件平臺的程序包方便地對水文序列的多種分布參數(shù)(如均值、均方差和偏態(tài)系數(shù)等)進行非一致性建模。因此，本文采用GAMLSS模型，分別以時間t，4種氣溫形式(年平均氣溫Tave、年平均最高氣溫Tmax、年平均最低氣溫Tmin和年平均溫差ΔT)為協(xié)變量對三峽庫區(qū)的主汛期降雨序列進行非一致性研究，從時間上和空間上分別分析庫區(qū)Pmfs的分布參數(shù)(均值、均方差)的非一致性特征，為三峽庫區(qū)防洪減災(zāi)工作提供理論依據(jù)。

1 研究區(qū)概況及數(shù)據(jù)資料

三峽庫區(qū)(28°10′—31°50′N，105°10′—110°50′E)，位于長江上游末端，庫區(qū)面積約8.1萬km2，是我國典型暴雨中心之一。地形特點是溝壑縱橫，谷深水急，山體破碎，容易發(fā)生滑坡，生態(tài)環(huán)境十分脆弱。從地貌上看，三峽庫區(qū)處于川鄂湘黔隆起褶皺帶、川東平行嶺谷三大構(gòu)造單元和大巴山褶皺帶交匯處，地形以山地(約占75%)和丘陵(約占22%)為主。庫區(qū)雨量充沛，年降雨量為1 000～1 300 mm，主汛期(按6—8月)降雨量占全年降雨量的41%左右。多年平均氣溫17℃，氣溫垂直差異明顯。庫區(qū)由庫首(包含興山、秭歸、巴東、奉節(jié)氣象站)、庫腹(包含鎮(zhèn)坪、梁平、萬州、利川氣象站)和庫尾(包含江津、合川、沙坪壩、長壽氣象站)三部分組成。

本文所用數(shù)據(jù)資料來源于中國氣象局國家氣象信息中心提供的1961—2016年逐日降水、平均氣溫、最高氣溫和最低氣溫數(shù)據(jù)。根據(jù)研究需要，由原數(shù)據(jù)統(tǒng)計得到三峽庫區(qū)1961—2016年12個站點的主汛期(6—8月)降雨量、年平均氣溫Tave、年平均最高氣溫Tmax、年平均最低氣溫Tmin和年平均溫差ΔT作為備用數(shù)據(jù)。

2 GAMLSS模型

GAMLSS模型是由Rigby和Stasinopoulos(2005年)提出的(半)參數(shù)回歸模型。與傳統(tǒng)統(tǒng)計模型相比，GAMLSS模型建模更靈活[20]，主要體現(xiàn)在：(1) 分布的所有參數(shù)都可以建模為解釋變量的函數(shù)；(2) 可以靈活地對響應(yīng)變量使用多種函數(shù)分布類型，包括一系列高偏度和高峰度的離散分布和連續(xù)分布；(3) 在模型中它允許分布參數(shù)使用各種附加項。基于上述特征，使得GAMLSS模型可以有效地解決傳統(tǒng)指數(shù)分布族無法擬合具有超峰度和平頂峰度、高度正偏或負偏的隨機變量序列的難題。

2.1 模型的定義

在GAMLSS模型中，假設(shè)某一時刻t的觀測值yt服從概率分布函數(shù)F(Yt|θt)，其中θ=(θt1,θt2,…,θtp)是時刻t對應(yīng)的p個分布參數(shù)的向量。記任意時刻的第k個分布參數(shù)組成的向量θk=(θ1k,θ2k,…,θnk)T,(k=1,2, …,p)。若不考慮隨機效應(yīng)項，θk與解釋變量Xk的函數(shù)關(guān)系：

gk(θk)=Xkβk

(1)

式中:Xk為解釋變量矩陣,維度為n×Ik;βk=(β1k,β2k,…，βnk)是回歸參數(shù)向量,長度為Ik。

在GAMLSS模型中，位置參數(shù)向量被定義為第一參數(shù)向量θ1(用均值mu表示)，尺度參數(shù)向量被定義為第二參數(shù)向量θ2(用均方差或變差系數(shù)sigma表示)，如果分布中存在其他參數(shù)，就被定義為形狀參數(shù)(用nu表示)。

2.2 模型評價準則及殘差評價

GAMLSS采用全局擬合偏差GD來篩選最優(yōu)分布、分布參數(shù)與解釋變量的最優(yōu)函數(shù)類型，如下定義：

(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為樣本序列長度。

模型殘差評價是衡量模型擬合好壞的一個重要標準，其可以通過殘差正態(tài)QQ圖、worm plot圖、分位圖等定性驗證模型擬合效果，殘差評價原理及標準參見文獻[20]。

2.3 概率分布

根據(jù)1960—2016年三峽庫區(qū)主汛期降雨總量系列頻率分布特點，這里選用比較常見的Gumbel(GU)、Weibull(WEI)、Gamma(GA)、Logistic(LO)和Normal(NO)分布作為GAMLSS模型的備用分布函數(shù)。

3 結(jié)果與分析

3.1 Mann-Kendall趨勢檢驗

采用Mann-Kendall趨勢檢驗法對主汛期降雨量序列進行趨勢性分析，庫區(qū)12個站點有5個站點主汛期降雨量呈不顯著減小趨勢，7個站點主汛期降雨量呈增加趨勢，見圖1。從空間上看，三峽庫區(qū)主汛期降雨量由庫首秭歸地區(qū)的顯著增加(顯著性水平α=0.05)過渡到庫腹的不顯著減少，再到庫尾的不顯著增加。

圖1 庫區(qū)多年平均主汛期降雨量空間變化趨勢

3.2 基于協(xié)變量時間t的非一致性頻率分析

以時間t為協(xié)變量，建立主汛期降雨量序列分布參數(shù)與時間t的關(guān)系。使用9種組合方式和5種分布函數(shù)類型(共計45種擬合方式)對分布參數(shù)進行擬合，組合方式如下：(1) 均值(mu)、均方差(sigma)均為常數(shù)；(2) 均值(mu)為常數(shù)，均方差(sigma)隨時間t線性變化；(3) 均值(mu)為常數(shù)，均方差(sigma)隨時間t拋物線型變化；(4) 均值(mu)隨時間t線性變化，均方差(sigma)為常數(shù)；(5) 均值(mu)、均方差(sigma)均隨時間t線性變化；(6) 均值(mu)隨時間t線性變化，均方差(sigma)隨時間t拋物線型變化；(7) 均值(mu)隨時間t拋物線型變化，均方差(sigma)為常數(shù)；(8) 均值(mu)隨時間t拋物線型變化，均方差(sigma)隨時間t線性變化；(9) 均值(mu)、均方差(sigma)均隨時間t拋物線型變化。將AIC和SBC值作為擇優(yōu)的標準，AIC和SBC值越小，表示模型的模擬效果越好。前文模型評價準則及殘差評價中提到，AIC值為懲罰因子為2時的GAIC值，而SBC值為懲罰因子為ln(n)時的GAIC值[ln(56)=4.025，AIC值恒小于SBC值，差值為2.025·df]，因此擇優(yōu)標準AIC值和SBC值尋優(yōu)效果一致，下文選用AIC值作為典型指標進行迭代尋優(yōu)。

三峽庫區(qū)12個站點主汛期降雨量GAMLSS模型的AIC值見圖2。結(jié)果顯示：除長壽站最優(yōu)分布類型為Weibull外，其他11個站點的最優(yōu)分布均為Gamma分布；江津、萬州、奉節(jié)、秭歸、鎮(zhèn)坪和沙坪壩的AIC值擬合效果在考慮分布參數(shù)具有非一致性特征的情況下比認為分布參數(shù)為常數(shù)時要有所減小，6個站點AIC值依次減小了4.8，1.0，1.6，7.1，2.0，0.9，其中秭歸的改善幅度最大；其他6個站點AIC值在分布參數(shù)為常數(shù)時較小，可認為這些序列仍然是一致的。

上述根據(jù)AIC值準則，選出了各站點以時間t為協(xié)變量的最優(yōu)GAMLSS模型，需進一步驗證模型擬合的好壞。根據(jù)分布參數(shù)的回歸方程繪制各站點相應(yīng)的分位圖,統(tǒng)計結(jié)果見表1，結(jié)果表明：僅有長壽站、梁平站和利川站的實測點據(jù)頻率與相應(yīng)理論概率最大差值略微超過5%，其他的站點實測點據(jù)頻率與相應(yīng)理論概率最大差值均在5%以內(nèi)，表明統(tǒng)計參數(shù)方程與實測點的擬合關(guān)系良好。從分位圖中還可以得到，江津地區(qū)主汛期降雨量的方差不斷增大。秭歸、沙坪壩地區(qū)主汛期降雨量均值不斷增大，其中秭歸地區(qū)增長速度約為29.6 mm/10 a；萬州地區(qū)主汛期降雨量方差、奉節(jié)和鎮(zhèn)坪地區(qū)主汛期降雨量的均值在1990年以前越來越大，之后逐漸恢復(fù)。其他地區(qū)主汛期降雨量的均值和方差均為常數(shù)。

表1 以時間t為協(xié)變量的實測點據(jù)頻率與理論分位曲線概率對比 %

注：橫坐標1—9分別表示9種組合方式，其中1表示統(tǒng)計參數(shù)均為常數(shù)，為傳統(tǒng)的一致性模型，而2—9是分布參數(shù)與時間t相關(guān)，為非一致性模型。

雖然主汛期降雨序列采用了多種分布函數(shù)以及以時間t為協(xié)變量的多種組合方式得到了較優(yōu)的GAMLSS模型，但在氣候變化的背景下，缺乏相應(yīng)的物理意義。于是在這里進一步用與氣候變化密切相關(guān)的氣溫要素作為協(xié)變量擬合模型，研究其對主汛期降雨量搭建的GAMLSS模型是否具有修正意義。

3.3 基于協(xié)變量氣溫的非一致性頻率分析

在對年平均氣溫Tave、年平均最高氣溫Tmax、年平均最低氣溫Tmin和年平均溫差ΔT作為協(xié)變量分析前，先計算了它們分別與主汛期降雨量的線性相關(guān)系數(shù)r，見圖3。

圖3 響應(yīng)變量與各氣溫序列相關(guān)系數(shù)對比

結(jié)果顯示：主汛期降雨量與4種氣溫形式均呈負相關(guān)，其中長壽、梁平、萬州、鎮(zhèn)坪、利川、合川和沙坪壩共7個站點的主汛期降雨量與年均日溫差ΔT相關(guān)性很強，相關(guān)系數(shù)絕對范圍為0.31～0.42；奉節(jié)、興山和秭歸的主汛期降雨量序列與年平均氣溫Tave相關(guān)性最強，相關(guān)系數(shù)分別為-0.22，-0.29和-0.53；江津、巴東站的類似相關(guān)系數(shù)絕對值均沒有超過0.2，說明從線性趨勢上來看，這兩個站的主汛期降雨量序列受氣溫影響不明顯。

與前文以時間t為協(xié)變量時的9種組合方式相同(各分布參數(shù)與氣溫的組合方式：常數(shù)、一次函數(shù)、拋物線關(guān)系)，分別以Tave，Tmax，Tmin和ΔT為協(xié)變量對各站點的Pmfs進行非一致性頻率分析。

根據(jù)AIC值最小原則，在GAMLSS模型里擇優(yōu)得到最優(yōu)分布函數(shù)回歸方程，見表2。結(jié)果表示：除了長壽站以ΔT為協(xié)變量擬合時最優(yōu)分布類型為Logistic分布外，其他組合方式的最優(yōu)分布均為Gamma分布(11/12)；長壽、梁平、萬州、鎮(zhèn)坪、利川、合川和沙坪壩在以平均日溫差為協(xié)變量擬合響應(yīng)變量時，AIC值最小，模擬效果最優(yōu)；奉節(jié)、巴東、興山、秭歸在以年平均氣溫為協(xié)變量擬合響應(yīng)變量時，模擬效果最優(yōu)，均與相關(guān)系數(shù)分析結(jié)果對應(yīng)。

表2 以氣溫因子為協(xié)變量的GAMLSS模型最優(yōu)分布參數(shù)

模型擬合的好壞需要通過殘差分析進行評價。首先對模型的殘差進行正態(tài)標準化處理，然后進行PPCC檢驗[20]。表3為殘差分布統(tǒng)計結(jié)果，表明各模型殘差序列分布情況接近標準正態(tài)分布，F(xiàn)illben相關(guān)系數(shù)也非常接近1，這表明模型擬合效果十分理想。

表3 以氣溫因子為協(xié)變量的殘差分布計算

以氣溫因子為協(xié)變量的回歸模型分位圖見圖4，江津地區(qū)主汛期降雨量的均方差增大；秭歸地區(qū)在1998年左右發(fā)生了均值突變，主汛期降雨量由從原來的480 mm突變至580 mm左右；其他站點的均值也出現(xiàn)了以氣溫為協(xié)變量的非一致性，但方差為常數(shù)(由于篇幅有限，其他站點未列出)。

3.4 模型擬合結(jié)果對比及不確定性分析

從GAMLSS模型運算結(jié)果中，得到位置參數(shù)μ、尺度參數(shù)σ為常數(shù)或在不同協(xié)變量下的AIC值，見圖5。(1) 長壽、梁平、巴東、興山、利川和合川的最優(yōu)AIC值比較：AIC(時間t為協(xié)變量的最優(yōu)分布組合)>AIC(分布參數(shù)為常數(shù)的最優(yōu)分布)>AIC(氣溫為協(xié)變量的最優(yōu)分布組合)，說明這6個站點在以時間t為協(xié)變量擬合時分布參數(shù)是固定不變的，但考慮氣溫因子后，其表現(xiàn)出分布參數(shù)隨氣溫變化的非一致性。(2) 江津、萬州、秭歸、鎮(zhèn)坪和沙坪壩最優(yōu)AIC值比較：AIC(分布參數(shù)為常數(shù)的最優(yōu)分布)>AIC(時間t為協(xié)變量的最優(yōu)分布組合)>AIC(氣溫為協(xié)變量的最優(yōu)分布組合)，這說明這5個站點的分布參數(shù)不僅對時間t變化，也隨氣溫變化，但考慮以氣溫為協(xié)變量擬合模型不僅具有物理意義，擬合效果也最優(yōu)。(3) 奉節(jié)站的最優(yōu)AIC值規(guī)律：AIC(t為協(xié)變量的最優(yōu)分布組合)>AIC(氣溫為協(xié)變量的最優(yōu)分布組合)>AIC(分布參數(shù)為常數(shù)的最優(yōu)分布)，奉節(jié)站相應(yīng)變量的分布參數(shù)既沒隨時間變化也沒隨氣溫變化，可能人類活動或者其他氣象因素對其的干擾掩蓋了氣溫對主汛期降雨量的影響，有待進一步驗證。

注：A代表以時間為協(xié)變量；B代表氣溫為協(xié)變量。

根據(jù)分布參數(shù)非一致特征的空間分布情況，可以得出：(1) 以t為協(xié)變量時，庫尾江津地區(qū)主汛期降雨量的方差不斷增大，主汛期降雨量的量級的可能增大；庫首秭歸地區(qū)、庫尾沙坪壩地區(qū)主汛期降雨量均值不斷增大，其中秭歸地區(qū)增長速度約為29.6 mm/10 a，理論上這兩地主汛期降雨量未來會越來多。(2) 以氣溫為協(xié)變量時，庫首各站主汛期降雨與年均日溫差ΔT相關(guān)系數(shù)最高(-0.42～-0.31)，并在以ΔT為協(xié)變量時模型擬合效果最好；庫區(qū)中部和西南角各站點主汛期降雨與年平均日氣溫Tave相關(guān)系數(shù)最高(-0.53～0.22)，并在以Tave為協(xié)變量時模型擬合效果最好。庫尾江津地區(qū)在以Tave為協(xié)變量時模型擬合效果最好。

本研究以氣溫要素為協(xié)變量對主汛期降雨量進行了非一致性分析，分位圖中部分站點的經(jīng)驗頻率與理論概率的差值仍在5%以上，說明還有其他因素或?qū)χ餮雌诮涤甑姆且恢滦杂绊戄^大，例如：土地利用變化因素[21]、溫室氣體排量因素[22]、自然氣候變化因素[23](北太平洋濤動(NPO)、太平洋年代際濤動(PDO)、北極濤動(AO)、南方濤動(SO)和東部熱帶太平洋海表溫度異常等)，這些因素是否能夠改善模型，有待進一步驗證。

注：1表示分布參數(shù)為常數(shù)；t表示以時間為協(xié)變量擬合分布參數(shù)。

4 結(jié) 論

(1) 從趨勢檢驗上看，從1961—2016年，三峽庫區(qū)12個氣象站點中有5個站點的主汛期降雨量(Pmfs)呈不顯著減小趨勢，7個站點主汛期降雨量呈不顯著增加趨勢。從空間上看，三峽庫區(qū)主汛期降雨量由庫首的顯著增加(顯著性水平α=0.05)到庫腹的不顯著減少，再過渡到庫尾的不顯著增加。

(2) 以t為協(xié)變量擬合時，主汛期降雨量(Pmfs)的均值、均方差具有非一致性的站點個數(shù)分別為4，2。其中庫尾江津地區(qū)主汛期降雨量的均方差不斷增大，主汛期降雨量的量級的在未來可能增大；庫首的秭歸地區(qū)、庫尾的沙坪壩地區(qū)主汛期降雨量均值不斷增大(秭歸地區(qū)增長速度約為29.6 mm/10 a)，理論上這兩地主汛期降雨量逐漸增加。這些地區(qū)在未來可能引發(fā)洪澇、滑坡、干旱等自然災(zāi)害。

(3) 以氣溫為協(xié)變量擬合時，庫區(qū)站點主汛期降雨量均呈現(xiàn)非一致性(均值、均方差具有非一致性特征的站點個數(shù)為10，2)，模型總體擬合效果優(yōu)于前者且具有物理意義。庫首以Tave為協(xié)變量時，模型擬合效果最好；庫腹、庫尾總體上以ΔT為協(xié)變量時，擬合效果最優(yōu)。值得注意的是，庫首秭歸地區(qū)主汛期降雨量在1998年左右發(fā)生了均值突變，均值由從原來的480 mm突增至580 mm左右，可能是由于人類活動或氣候變化導(dǎo)致年平均氣溫Tave變化引起的；庫尾的江津地區(qū)主汛期降雨量的均方差呈增大趨勢，這兩地出現(xiàn)自然災(zāi)害的風險同以t為協(xié)變量的擬合結(jié)果。

(4) 3種模型(分布參數(shù)均為常數(shù)、以時間t為協(xié)變量和以氣溫為協(xié)變量)的比較表明，有近92%的站點最優(yōu)分布函數(shù)類型為Gamma分布，且以氣溫為協(xié)變量的非一致模型更好地反映了觀測數(shù)據(jù)的變化。這些結(jié)果證實，在主汛期降雨序列中非一致分析考慮物理協(xié)變量是必要且有效的，可以根據(jù)氣溫因素來減小三峽庫區(qū)主汛期降雨量非一致性分析的不確定性。