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

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

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

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

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

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

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

1 研究區概況及數據資料

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

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

2 GAMLSS模型

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

2.1 模型的定義

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

gk(θk)=Xkβk

(1)

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

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

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

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

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

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

2.3 概率分布

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

3 結果與分析

3.1 Mann-Kendall趨勢檢驗

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

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

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

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

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

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

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

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

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

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

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

圖3 響應變量與各氣溫序列相關系數對比

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

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

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

表2 以氣溫因子為協變量的GAMLSS模型最優分布參數

模型擬合的好壞需要通過殘差分析進行評價。首先對模型的殘差進行正態標準化處理，然后進行PPCC檢驗[20]。表3為殘差分布統計結果，表明各模型殘差序列分布情況接近標準正態分布，Fillben相關系數也非常接近1，這表明模型擬合效果十分理想。

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

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

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

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

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

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

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

注：1表示分布參數為常數；t表示以時間為協變量擬合分布參數。

4 結 論

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

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

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

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