李 婧, 閆 磊, 屈春來, 劉章君, 魯東陽

(1.河北工程大學水利水電學院，河北 邯鄲 056000； 2.江西省水利科學研究院，江西 南昌 330029)

洪水頻率分析對水利工程建設(shè)至關(guān)重要，一致性假設(shè)一直是洪水頻率分析的基本假設(shè)條件，但受到氣候變化以及人類活動的影響，降雨時空分配過程以及流域下墊面的產(chǎn)匯流條件發(fā)生了改變，導致洪水頻率分析的一致性假設(shè)受到挑戰(zhàn)[1-4],傳統(tǒng)頻率計算方法構(gòu)建模型的合理性受到質(zhì)疑。因此，在變化環(huán)境下進行洪水序列的非一致性頻率分析具有重要意義。非一致性頻率分析成為變化環(huán)境下洪水分析計算的研究熱點之一[2,5]。當前最常用的非一致性頻率分析方法是采用基于位置、尺度和形狀參數(shù)的廣義可加模型(generalized additive models for location scale and shape, GAMLSS)構(gòu)建非一致性概率分布模型。Villarini等[6-7]應用GAMLSS模型分別對羅馬降水序列以及美國Little Sugar Creek流域的年最大洪峰序列進行了趨勢分析，取得了較好的擬合效果。江聰?shù)萚8]認為運用傳統(tǒng)方法對水文序列均值的趨勢性進行分析相對片面，通過GAMLSS模型引入方差和偏態(tài)系數(shù)等參數(shù)對宜昌站多年年流量序列進行了趨勢分析。張冬冬等[9]采用GAMLSS模型研究了大渡河流域年最大日降水極值序列的非一致性，發(fā)現(xiàn)GAMLSS模型可以模擬水文序列的方差、偏態(tài)系數(shù)等統(tǒng)計參數(shù)的非線性變化趨勢。

然而，GAMLSS模型構(gòu)建過程較復雜，首先需要選定概率分布模型，然后假定所選概率分布模型統(tǒng)計參數(shù)與協(xié)變量的關(guān)系，模型統(tǒng)計參數(shù)多且結(jié)構(gòu)較復雜。為解決這一問題，董潔等[10-12]嘗試引入非參數(shù)方法進行非一致性水文頻率分析。Koenker等[13]首先提出了分位數(shù)回歸方法研究樣本點在不同的分位數(shù)下隨時間變化情況，隨著該方法的進一步研究開發(fā)，目前廣泛應用于水文分析中。Mazvimavi[11]利用分位數(shù)回歸方法對津巴布韋的降雨序列變化進行分析。Barbosa[12]等運用分位數(shù)回歸對波羅的海海平面變化進行分析，發(fā)現(xiàn)極大值處的斜率趨勢性更大。Wang等[14]利用分位數(shù)回歸方法對美國南部月降雨變化情況進行分析。馮平等[15]利用分位數(shù)回歸方法分析灤河流域年降雨、徑流序列的變化特性，發(fā)現(xiàn)流域年降雨序列總體保持一致性，而年徑流序列受人類活動等影響表現(xiàn)為減小趨勢，且高分位數(shù)的減小速率大于低分位數(shù)。本文選取位于我國不同氣候區(qū)的渭河流域和珠江流域的5個代表站點的年最大洪水序列作為研究對象，綜合定性分析和定量分析，對比分位數(shù)回歸模型與GAMLSS模型在進行非一致性洪水頻率分析時的表現(xiàn)。

1 流域概況和數(shù)據(jù)資料

渭河流域發(fā)源于甘肅省渭源縣鳥鼠山，流經(jīng)甘肅、陜西兩省，至渭南市潼關(guān)縣匯入黃河，河長818 km，流域面積134 766 km2，位于北緯33°40′～37°26′，東徑103°57′～110°27′之間[16-17]。渭河中下游徑流的年際變化表現(xiàn)為南部小、北部大。渭河干流秋季徑流量最大，約占年徑流量的38%～40%，夏季占32.8%～34.2%，春季占17.7%～19.1%，冬季為8.3%～9.9%。渭河中下游降雨集中于7—9月，而且多大暴雨，洪澇災害較多。

珠江發(fā)源于云南省曲靖市馬雄山，流經(jīng)云南、貴州、廣西、廣東、湖南、江西6省及越南北部，在下游從8個入海口注入南海，全長2 320 km，流域面積45 369 km2。珠江流域由西江、北江、東江和珠江三角洲諸河等4個水系組成[18-19]，流域多年平均降水量1 200～2 200 mm，年徑流量超3 300億m3，汛期4—9月徑流量約占年徑流量的80%，夏季(6—8月)則占年徑流量的50%以上。

數(shù)據(jù)資料選取渭河的咸陽站、華縣站、張家山站以及珠江流域的高道站和大湟江口站共5個重要站點的多年年最大洪水序列。5個站點的基本概況詳見表1及圖1。

圖1 流域水文站點示意圖

表1 各站點基本概況

2 研究方法

2.1 Mann-Kendall趨勢檢驗方法

采用Mann-Kendall非參數(shù)趨勢檢驗法檢測降水和徑流序列的長期變化趨勢[20]。Mann-Kendall的檢驗方法[15,21]如下：

設(shè)Yi(i=1,2,…,n)為假設(shè)檢驗的隨機變量，n為樣本的觀測長度，標準化檢驗統(tǒng)計量Zc為

(1)

式中：P′為序列中所有對偶觀測值(Yi,Yj,i0時則表示該樣本序列呈現(xiàn)上升趨勢，反之則表現(xiàn)為下降趨勢。

2.2 分位數(shù)回歸法

分位數(shù)回歸是依據(jù)因變量Y的條件分位數(shù)對自變量X進行回歸，得到所有分位數(shù)下的回歸模型；分位數(shù)回歸對于模型中的隨機干擾項不需要做概率分布假定，且估計量基本不受異常值的影響而更加穩(wěn)定，可擬合出一簇曲線[15]。

假設(shè)隨機變量Y的分布函數(shù)為F(y)=P(Y≤y)，則Y的第τ分位數(shù)為

Q(τ)=inf{y:F(y)≥τ} (0<τ<1)

(2)

根據(jù)不同的分位數(shù)τ，可得到不同的分位數(shù)函數(shù)QY(τ)=ωiβ(τ)(QY(τ)為第τ分位數(shù)時Y的分位數(shù)函數(shù)，β(τ)為參數(shù)值,i=1,2,…,n)。分位數(shù)回歸的目標函數(shù)是一個加權(quán)的絕對殘差和，目的是使樣本值與擬合值之間的距離最短。給定觀測數(shù)據(jù)(ω1,y1),(ω2,y2),…,(ωn,yn)，由此求解τ分位數(shù)的回歸估計：

(3)

其中ρτ(·)是一個非對稱的損失函數(shù)，定義為

可通過下式得到參數(shù)估計值，根據(jù)一個τ值確定一組β(τ)值：

(5)

2.3 GAMLSS模型

GAMLSS模型是基于位置、尺度、形狀的廣義可加模型，可以描述統(tǒng)計參數(shù)與協(xié)變量之間的線性和非線性關(guān)系，模擬隨機變量分布函數(shù)類型廣泛，適用于高偏度、高峰度、超峰度、平頂峰度等的隨機變量序列[8,22]。

2.3.1模型的定義

在GAMLSS模型中，假設(shè)某一時刻t(t=1,2,…,n)的相對獨立隨機變量觀測值yt的分布函數(shù)可以表示為F(yt|θt)，其中θt=(θt1,θt2,…,θtf)為時刻t對應的分布統(tǒng)計參數(shù)向量，f為分布參數(shù)的個數(shù)，n為觀測值的個數(shù)[8]。記gk(θk)為θk與相應的協(xié)變量Yk的單調(diào)函數(shù)關(guān)系，一般表示為

(6)

式中：ηk為長度n的向量；βk=(β1k,β2k,…,βIkk)T，為長度為Ik的回歸參數(shù)向量；Yk為n×Ik的協(xié)變量矩陣；Zjk為已知的n×qjk矩陣；γjk為qjk維的正態(tài)分布隨機變量向量；Zjkγjk為第j項隨機效應[8]。

若不考慮隨機效應的影響，則gk(θk) =ηk=Ykβk，針對位置參數(shù)μ、尺度參數(shù)σ和形狀參數(shù)υ，以時間t為協(xié)變量的不考慮隨機效應的全參數(shù)模型為

(7)

式中：μt、σt、υt分別為時變的位置參數(shù)、尺度參數(shù)、形狀參數(shù)，可體現(xiàn)非平穩(wěn)序列統(tǒng)計參數(shù)隨時間的變化特征，因為形狀參數(shù)對模型影響較大，所以一般不做考慮，可通過RS法[23]對β進行參數(shù)估計。

2.3.2模型的評價準則

選取GAIC準則(generalized akaike information criterion，簡稱GAIC)作為GAMLSS模型擬合評價指標，其值VGAIC的計算公式為

VGAIC=VGD+#df

(8)

2.4 擬合優(yōu)度檢驗

2.4.1正態(tài)QQ圖

通過正態(tài)QQ圖來分析最優(yōu)模型的殘差分布狀況。正態(tài)QQ圖是在平面坐標系中以經(jīng)驗殘差為縱坐標、理論殘差為橫坐標繪制而成的，數(shù)據(jù)點與1∶1的線偏差越小，說明模型的擬合優(yōu)度越好。

2.4.2模型概率覆蓋率檢驗

通過模型概率覆蓋率偏差數(shù)值的大小定性分析分位數(shù)回歸模型與GAMLSS模型的擬合優(yōu)度。模型概率覆蓋率是各條分位數(shù)曲線覆蓋范圍內(nèi)的實測樣本點數(shù)量與總樣本點數(shù)量的比值。計算模型概率覆蓋率與對應分位數(shù)的差值，差值越小說明模型擬合優(yōu)度越好。以50%分位數(shù)曲線為例，理論上講，如果50%分位數(shù)曲線覆蓋范圍內(nèi)的實測樣本點數(shù)量與總樣本點數(shù)量的比值越接近50%，則模型概率覆蓋率與50%差值越小，模型擬和優(yōu)度越好。

2.4.3Filliben檢驗

可通過Filliben相關(guān)系數(shù)來確定洪水序列最優(yōu)分布形式，通過Filliben相關(guān)系數(shù)的大小來判定序列最優(yōu)擬合分布，F(xiàn)illiben值越大則分布擬合優(yōu)度越好[24]。

(9)

3 結(jié)果和討論

3.1 Mann-Kendall趨勢分析

使用R語言中trend程序包對各站點歷史資料進行Mann-Kendall檢驗，顯著性水平設(shè)定為5%，得出華縣站、高道站、大湟江口站、咸陽站、張家山站5個站點P、Zc和S值見表2。

表2 各站點趨勢性分析結(jié)果

華縣站、高道站、咸陽站、張家山站4個站點通過顯著性水平5%的檢定，而大湟江口站通過顯著性水平為10%的檢定，同時根據(jù)Zc值正負關(guān)系可得出高道和大湟江口2個站點呈上升趨勢，其他3個站點則呈下降趨勢。故華縣站、咸陽站、張家山站3個站點顯著下降，高道站顯著上升，而大湟江口站呈現(xiàn)不顯著上升趨勢。

3.2 GAMLSS最優(yōu)模型擬合效果分析

3.2.1GAMLSS最優(yōu)模型判定

本算例為選定最優(yōu)概率分布類型，選取我國不同氣候區(qū)的渭河、珠江流域的5個站點最大洪水序列為研究對象，時間t為協(xié)變量，AIC為評價準則，分別對比了P-Ⅲ分布、Gamma分布(簡稱GA分布)、Lognormal分布(簡稱LN分布)和GEV分布4種分布形式，通過正態(tài)QQ圖判別擬合優(yōu)度。每個站點每種分布類型根據(jù)位置參數(shù)和尺度參數(shù)分別隨時間變化可形成4種不同的模型，總計16種，相應AIC取值詳見表3。在表3模型代碼中，以“L”表示位置參數(shù)，“S”表示尺度參數(shù)，“0”表示參數(shù)不變，“t”表示參數(shù)隨時間協(xié)變量變化，此外粗體AIC值表示該站點的最佳模型。

表3 各站點模型AIC值對比

由此可得出華縣站、高道站、咸陽站的最優(yōu)模型為位置參數(shù)和尺度參數(shù)均隨時間變化的Gamma分布，AIC值分別為1 051.84、1 057.64、955.81；大湟江口站最優(yōu)模型則為位置參數(shù)隨時間變化而尺度參數(shù)保持不變的Gamma分布，AIC值為1 155.26；張家山站最優(yōu)模型則為位置參數(shù)隨時間變化而尺度參數(shù)保持不變的Lognormal分布，AIC值為894.26。

圖2給出了最優(yōu)模型的正態(tài)QQ圖的擬合優(yōu)度檢驗結(jié)果。結(jié)果顯示，根據(jù)AIC值選擇的最優(yōu)模型點均勻分布在1∶1直線附近，經(jīng)驗殘差和理論殘差基本一致，表明模型具有良好的擬合效果。

圖2 各站點QQ圖

3.2.2GAMLSS最優(yōu)模型的定性分析

對于該算例的實際樣本點來說，通過對比各站點GAMLSS模型的分位數(shù)曲線圖，定性分析模型概率覆蓋率，可得出華縣站GAMLSS模型的概率覆蓋率偏差最小為0.16%，最大為7.26%，高道站、大湟江口站、咸陽站、張家山站GAMLSS模型的偏差范圍分別為1.56%～6.97%、1.43%～5.36%、0.17%～5.17%和0.08%～2.97%，具體見圖3與表4。

圖3 各站點基于GAMLSS模型的分位數(shù)曲線

表4 各站點不同分位數(shù)GAMLSS模型擬合優(yōu)度對比

3.2.3GAMLSS最優(yōu)模型的定量分析

通過對比Filliben相關(guān)系數(shù)進行定量對比分析，基于式(9)進行相應計算可得GAMLSS模型的Filliben相關(guān)系數(shù)：華縣站為0.989 96；高道站為0.987 843；大湟江口站為0.986 808；咸陽站為0.994 089；張家山站為0.991 817。

3.3 分位數(shù)回歸模型擬合效果分析

3.3.1分位數(shù)回歸模型的定性分析

通過對比各站點分位數(shù)回歸模型的分位數(shù)曲線圖，可定性分析模型概率覆蓋率，可得出華縣站分位數(shù)回歸模型下模型概率覆蓋率偏差最小為0.81%，最大為4.84%；高道站、大湟江口站、咸陽站、張家山站分位數(shù)回歸模型下偏差范圍分別為0.08%～5.74%、0～0.36%、、0～2.59%和0.08%～1.61%，具體見圖4與表5。

表5 各站點不同分位數(shù)回歸模型擬合優(yōu)度對比

圖4 各水文站基于分位數(shù)回歸模型的分位數(shù)線

3.3.2分位數(shù)回歸模型的定量分析

可通過對比Filliben相關(guān)系數(shù)進行定量對比分析，基于式(9)進行相應計算可得分位數(shù)回歸模型的Filliben相關(guān)系數(shù)：華縣站，0.983 076；高道站，0.983 149；大湟江口站，0.983 534；咸陽站，0.983 022；張家山站，0.984 437。

3.4 分位數(shù)回歸模型與GAMLSS模型擬合效果比較

根據(jù)模型概率覆蓋率進行定性分析，可發(fā)現(xiàn)各站點分位數(shù)回歸模型均優(yōu)于GAMLSS模型，華縣站顯示優(yōu)度為0.3%～6.4%，高道站為1.6%～3.3%，大湟江口站為1.1%～5.4%，咸陽站為1.4%～5.2%，張家山站為0～1.8%；根據(jù)Filliben相關(guān)系數(shù)進行定量分析，兩種模型Filliben相關(guān)系數(shù)較高，均通過顯著性水平5%時，F(xiàn)illiben相關(guān)系數(shù)臨界值為0.975的檢驗[25-28]，兩種模型的擬合效果GAMLSS模型優(yōu)于分位數(shù)回歸模型的范圍基本位于0.3%～0.9%之間?；贔illiben相關(guān)系數(shù)的定量分析結(jié)果顯示，GAMLSS模型擬合效果更好但優(yōu)勢不明顯；基于模型概率覆蓋率的定性分析結(jié)果表明，分位數(shù)回歸模型的擬合效果更好。因此，綜合上述定性分析和定量分析結(jié)果，總體來說，分位數(shù)回歸模型的擬合效果更優(yōu)。

4 結(jié) 論

a. 通過Mann-Kendall趨勢檢驗可以得出5%顯著性水平下華縣站、咸陽站、張家山站3個站點顯著下降，高道站顯著上升，而大湟江口站呈不顯著上升趨勢。

b. 各站點在GAMLSS模型中擬合效果最好的分布，以Gamma分布為主、Lognormal分布為輔。華縣站、高道站、咸陽站表現(xiàn)為Gamma分布下的位置參數(shù)和尺度參數(shù)均隨時間變化的模型效果最佳，大湟江口站則為位置參數(shù)隨時間變化而尺度參數(shù)保持不變的Gamma分布效果最佳，張家山站表現(xiàn)為Lognormal分布下的位置參數(shù)隨時間變化而尺度參數(shù)保持不變的模型最佳。

c. 通過定性分析對比各站點兩種模型概率覆蓋率，發(fā)現(xiàn)分位數(shù)回歸模型整體擬合效果優(yōu)于GAMLSS模型。定量分析Filliben相關(guān)系數(shù)，分位數(shù)回歸模型與GAMLSS模型擬合效果基本一致。故綜合定性分析和定量分析可得出分位數(shù)回歸模型的擬合效果優(yōu)于GAMLSS模型。

d. GAMLSS模型可以考慮多種因素的影響以及更多的分布類型，進行更加復雜的模擬運算，但是統(tǒng)計參數(shù)較多，且統(tǒng)計參數(shù)與協(xié)變量的關(guān)系可能不只是線性，實際變化情況更加復雜，使得模型構(gòu)建過程較為煩瑣，而基于非參數(shù)方法的分位數(shù)回歸模型構(gòu)建更加簡單且易于操作，模型擬合更加靈活，無須確定各項參數(shù)以及參數(shù)與協(xié)變量之間的復雜關(guān)系。