高 潔

(水電水利規(guī)劃設計總院，北京100120)

1 研究背景

水文系列的平穩(wěn)性是工程設計中頻率計算、抽樣統(tǒng)計的基礎。在工程實踐中尤其重視水文系列的三性檢驗，要求滿足可靠性、一致性和代表性。當氣候變化、人類活動等多種因素影響系列的一致性時，通常進行還原計算，以保證采用平穩(wěn)序列進行頻率適線。但是，隨著人類活動加劇，眾多江河水文情勢發(fā)生顯著變化，水文極端事件頻發(fā)[1]；區(qū)域性的非平穩(wěn)序列廣泛存在，非一致性問題并非個案，流域水文律情的“一致性”背景已不復存在?，F(xiàn)有的基于平穩(wěn)序列的工程水文分析計算方法將面臨變化環(huán)境帶來的設計頻率失真風險[2]。在普遍非平穩(wěn)的環(huán)境中，需要采用針對非平穩(wěn)序列的統(tǒng)計模型和計算流程[2]。

針對非平穩(wěn)序列統(tǒng)計參數(shù)的時變特征，Strupczewski等[3- 5]通過對時變一階矩、二階矩采用極大似然法、加權最小二乘法進行參數(shù)估計；并基于AIC準則優(yōu)選模型的原理和方法，以波蘭河流洪峰流量系列為例，在Mann-Kendall趨勢檢驗的基礎上，分析序列平穩(wěn)性及統(tǒng)計參數(shù)變化趨勢。葉長青等[6]采用時變矩模型，通過5種分布、8種趨勢選擇最優(yōu)擬合模式，分析武江坪石水文站1964年～2008年和東江龍川水文站1953年～2008年歷年最大日流量系列。研究表明，兩者均符合對數(shù)正態(tài)分布，分別滿足均值均方差線性變化、均值均方差拋物線型變化模式?；谄椒€(wěn)序列坪石站、龍川站100年一遇洪峰流量分別為4 170、6 020 m3/s；采用時變矩模型后，坪石站4 170 m3/s的洪水重現(xiàn)期從1970年前大于200年一遇變?yōu)?000年以后小于60年一遇。受水庫調(diào)蓄影響，龍川站6 020 m3/s洪水重現(xiàn)期從1962年前小于100年一遇變?yōu)?000年以后大于400年一遇，研究認為傳統(tǒng)平穩(wěn)序列的洪水重現(xiàn)期概念應修正。劉丙軍等[7]通過時變矩模型分析西江和北江水系的馬口水文站和三水水文站1960年～2009年歷年最大日流量系列，發(fā)現(xiàn)馬口站符合GLO分布，均值拋物線型變化，均方差與變差系數(shù)Cv相關；三水站符合GEV分布，均值拋物線型變化，均方差平穩(wěn)不變。馬口站100年一遇洪水洪峰洪量，1960年為56 000 m3/s，1975年為52 000 m3/s，2009年增至67 000 m3/s；三水站100年一遇洪水洪峰洪量1960年為15 000 m3/s，1990年為16 000 m3/s，2009年增至20 000 m3/s。同一頻率設計洪水在20世紀70年代較小，21世紀后明顯增大，差異化明顯。

考慮到城市化進程對小流域洪水的影響，Villarini等[8]以美國北卡羅來納州集水面積110 km2的Littele Sugar Creek為對象，通過GAMLSS(Generalized Additive Models for Location， Scale and Shape)研究發(fā)現(xiàn)，在83年時間序列中，傳統(tǒng)方法計算平穩(wěn)序列100年一遇洪水洪峰流量設計值為3.2 m3/s。非平穩(wěn)序列分析發(fā)現(xiàn)，設計流量3.2 m3/s的重現(xiàn)期從1957年的5 000年一遇減小為2007年的8年一遇，與傳統(tǒng)頻率計算成果差異較大。

本文詳細介紹基于開源的專業(yè)語言R及軟件平臺RStudio開發(fā)的GAMLSS[9]在非平穩(wěn)序列中的應用，以美國Little Sugar Creek的經(jīng)典案例為例，闡明GAMLSS在非平穩(wěn)時間序列分析中的流程方法。

2 研究基礎

2.1 廣義可加模型GAM

廣義可加模型(Generalized Additive Model，GAM)是在廣義線性模型(Generalized Linear Model，GLM)和可加模型(Additive Model)的基礎上，由Hastie和Tibshirani于1990年提出。模型的解釋變量不限于正態(tài)分布，響應變量形式靈活，采用非參數(shù)的方法進行擬合，將存在復雜非線性關系的解釋變量通過不同的函數(shù)形式，以加和的方式構(gòu)成模型，可反映變量間非單調(diào)、非線性關系[10- 12]，包括了可加項(additive component)、隨機項(random component)和聯(lián)系兩者的連接函數(shù)(link function)。即

廣義線性模型

g(θ)=η=α+β1X1+…+βJkXJk

(1)

可加模型

E(Y/X1，…，XJk)=α+h(x1)+…+h(xJk)

(2)

廣義可加模型

(3)

式中，η為預測值或觀測值；X為解釋變量；Jk為解釋變量的個數(shù)；α為常數(shù)項；β為解釋變量的線性參數(shù)；g(·)為連接函數(shù);θ=E(Y/X1，…，XJk)為Y的數(shù)學期望值；hj(·)為平滑函數(shù)，包括光滑樣條函數(shù)、核函數(shù)、局部回歸平滑函數(shù)等各種不同形式。其靈活非參數(shù)形式，可反映解釋變量的非線性效應[13]。

2.2 基于位置、尺度和形狀參數(shù)的廣義可加模型GAMLSS

GAMLSS模型是由Rigby和Stasinopoulos于2005年提出的(半)參數(shù)回歸模型。GAMLSS模型在擬合位置、尺度和形狀參數(shù)的基礎上，通過可加的半?yún)?shù)項或非參數(shù)項、隨機效應項，建立響應變量統(tǒng)計參數(shù)(位置、尺度、形狀等)與解釋變量的關系[14]。GAMLSS模型在醫(yī)學、經(jīng)濟學[15- 16]等領域已得到廣泛應用，對于超峰度、平頂峰度、高度正偏或負偏等高偏度和高峰度的隨機變量系列具有較好的擬合效果[17]。江聰和熊立華[17]將GAMLSS模型應用于宜昌站1882年～2009年歷年平均和年最小月流量系列趨勢分析。研究結(jié)果表明，該站的年最小月流量系列非平穩(wěn)，均值呈非線性變化，偏度系數(shù)呈線性變化。

2.2.1模型I——不考慮隨機項

GAMLSS模型假設獨立隨機變量第i時刻觀測值的概率密度函數(shù)為f(yi|θi)。其中，θi可通過位置參數(shù)(以均值、中位數(shù)等表征)、尺度參數(shù)(方差均方差等表征)和形狀參數(shù)(偏度系數(shù)，峰度系數(shù))反映。即

(4)

針對位置參數(shù)θ1和尺度參數(shù)θ2的兩參數(shù)模型，以時間t為解釋變量不考慮隨機效應的全參數(shù)模型

g1(θ1)=η1=X1β1=β11+β21t+…+βI1tI1-1

(5)

g2(θ2)=η2=X2β2=β12+β22t+…+βI2tI2-1

(6)

式中，可設定解釋變量矩陣Xk為時間線性相關或拋物線型相關矩陣，并通過RS法[3- 4]對β進行參數(shù)估計。

2.2.2模型II——半?yún)?shù)化可加項

GAMLSS還引入一些子模型，如半?yún)?shù)可加模型、非線性半?yún)?shù)模型、非線性參數(shù)模型等[9]。相同的分布函數(shù)有多種參數(shù)化方案，需要針對具體問題建模優(yōu)選。設Zjk=In，In是m×m的單位矩陣，γjk=hjk=hjk(xjk)，則GAMLSS半?yún)?shù)可加模型

(7)

式中，Xkβk是解釋變量的線性函數(shù)矩陣，hjk是解釋變量xjk的未知函數(shù)，xjk(j=1，2，…，Jk)均是n維向量，hjk(xjk)體現(xiàn)了對解釋變量模擬效果更為靈活的平滑函數(shù)。時間序列模擬如不考慮其他變量的線性可加項，僅保留時間變量的平滑函數(shù)，針對位置參數(shù)θ1和尺度參數(shù)θ2的半?yún)?shù)化模型可表示為

(8)

(9)

上述模型均可通過AIC(Akaike Information Criterion)或SBC(Schwartz Bayes Criterion)準則對擬合效果進行優(yōu)選，選出AIC值或SBC值最小者為最優(yōu)模型，最終求得非平穩(wěn)序列的統(tǒng)計參數(shù)分布函數(shù)及模型形式

(10)

AIC=-2lnML+2k

(11)

SBC=-2lnML+ln(n)

(12)

式中，n為序列長度(觀測數(shù)目)；ML為似然函數(shù)極大值；k為模型參數(shù)個數(shù)。

3 研究方法

3.1 研究步驟

通過GAMLSS模型對水文系列進行分析，可分為模型擬合、殘差判別和成果分析三個步驟。

(1)在模型擬合中，一是根據(jù)模型I的特點，以水文系列時間變量或大氣環(huán)流因子作為解釋變量，針對水文系列的位置參數(shù)、尺度參數(shù)等選擇不同的概率分布和時間變化趨勢進行組合，通過AIC準則選擇最優(yōu)模型。二是根據(jù)模型II的特點，在平滑函數(shù)中考慮時間變量或大氣環(huán)流因子，選擇不同的概率分布，通過AIC準則優(yōu)選位置參數(shù)、尺度參數(shù)等有效自由度，并確定模型。

(2)模型殘差判別，一是通過plot函數(shù)進行總體檢驗，二是通過worm plot檢驗正態(tài)標準化殘差序列是否符合標準正態(tài)分布。所有散點位于兩條橢圓線之間的95%置信區(qū)間，模型選擇合理，模擬結(jié)果可接受。

(3)在成果分析中，通過統(tǒng)計參數(shù)隨時間的變化趨勢，獲取不同時期各種頻率的流量值，以揭示同一流量不同時期測算的重現(xiàn)期之間的差異，反映水文頻率特征值是否隨時間變化及變化情況，為研究水文時間序列的非平穩(wěn)性提供基礎。

3.2 研究思路

(1)概率分布。根據(jù)河川年最大流量年際分布特點，本研究初選兩參數(shù)的Gamma、Lognormal(對數(shù)正態(tài))、Gumbel、Weibull和Normal(正態(tài))分布函數(shù)作為模型優(yōu)選的基礎。位置參數(shù)和尺度參數(shù)分別采用均值和均方差表征。

(2)模型優(yōu)選：①模型I，根據(jù)水文時間序列特點，綜合考慮模型擬合效果和預測外延性，選擇統(tǒng)計參數(shù)隨時間的變化函數(shù)。本研究將均值和均方差隨時間變化模式，簡化為三種情形：無趨勢(S)、線性變化(L)和拋物線變化(P)。根據(jù)AIC準則，對均值和均方差選擇最優(yōu)概率分布函數(shù)和時間變化模式。②模型II，根據(jù)水文時間序列特點，本研究選擇三次樣條插值函數(shù)作為平滑函數(shù)。通過AIC準則，分別選出均值、均方差擬合效果最優(yōu)的有效自由度。綜合考慮模型應用的預測外延性，在擬合效果基本不降低、模型殘差滿足要求的條件下，可對自由度進行微調(diào)和降維。

(3)殘差判別。通過殘差散點圖、QQ圖、worm plot圖等驗證模型擬合效果。

(4)成果分析。首先判斷水文時間序列是否平穩(wěn)，再基于統(tǒng)計參數(shù)隨時間的變化特征分析頻率設計值。

表1 概率分布

4 研究應用

本研究主要借鑒Villarini等[8]采用的美國Little Sugar Creek 1924年～2006年共83 a年最大洪峰流量數(shù)據(jù)，介紹GAMLSS模型進行非平穩(wěn)序列分析的方法和流程。

4.1 軟件基礎

采用開源的專業(yè)軟件平臺RStudio，下載并加載gamlss，gamlss.data及gamlss.dist軟件包和數(shù)據(jù)包，導入水文系列數(shù)據(jù)集。

4.2 概率分布

在gamlss.family中，Gamma、Lognormal、Gumbel、Weibull、Normal概率分布，分別對應于的GA()、LOGNO()、GU()、WEI()和NO()函數(shù)。

4.3 模型擬合

4.3.1模型I

(1)對位置參數(shù)(mu)和尺度參數(shù)(sigma)，采用不同的概率分布函數(shù)，分別按照平穩(wěn)過程(S)、隨時間線性變化(L)、隨時間拋物線型變化(P)進行模型擬合。

(2)選擇AIC值最小者為最佳模型，具體AIC值見表2。

由表2可知，采用Gamma分布和對數(shù)正態(tài)分布模擬效果均較好。Gamma分布中，位置參數(shù)拋物線型、尺度參數(shù)平穩(wěn)模型最優(yōu)，位置參數(shù)拋物線型、尺度參數(shù)拋物線型模型其次。經(jīng)多模型比較，尺度參數(shù)與時間不相關、線性相關、拋物線型相關的模擬效果差異不大，隨時間的變化規(guī)律不明顯。

表2 不同概率分布模式和變化趨勢模型AIC值

注：下劃線標識不同分布函數(shù)的最優(yōu)擬合值，斜體標識最優(yōu)分布函數(shù)的最優(yōu)擬合值。

4.3.2模型II

對位置參數(shù)和尺度參數(shù)，采用不同的概率分布函數(shù)，選取三次樣條平滑函數(shù)，通過迭代算法，推求均值和均方差最優(yōu)自由度，選擇AIC值最小者為最佳模型。

模型II仍以Gamma分布的擬合效果最優(yōu)(見表3)，位置參數(shù)和尺度參數(shù)的自由度分別為2.0和0.9。根據(jù)Villarini等[8]研究采用的正態(tài)分布模擬結(jié)果，其流量系列均值和均方差的自由度分別為1.71和5.02，與本文相應正態(tài)分布的結(jié)果(1.7，5.1)基本一致。Villarini等為簡化模型復雜度，在不顯著降低模擬效果的情況下，最終將尺度參數(shù)的自由度調(diào)整為1.7。根據(jù)AIC值初步判斷，本文所選用的Gamma分布更優(yōu)于Villarini采用的正態(tài)分布擬合方式。

表3 不同概率分布和自由度模型AIC值

注：下劃線和斜體標識分別為最優(yōu)分布函數(shù)和自由度的最優(yōu)擬合值。

4.4 殘差判別

通過plot函數(shù)統(tǒng)計殘差的均值、方差、偏度系數(shù)、峰度系數(shù)、Filliben相關系數(shù)，并根據(jù)核密度估計圖、QQ圖等檢驗正態(tài)標準化殘差是否服從標準正態(tài)分布。當模型殘差服從標準正態(tài)分布時，均值、方差、偏度系數(shù)和峰度系數(shù)分別接近0、1、0和3。經(jīng)檢驗，模型殘差滿足標準正態(tài)分布要求。

4.5 模型成果

(1)根據(jù)模擬效果較優(yōu)的模型I(mu-P-Sigma-S)和模型II(Gamma(2.0，0.9))，該流域年最大流量分布均值(位置參數(shù))隨時間先減小后增大，在1950年左右出現(xiàn)轉(zhuǎn)折。1950年之前，流量呈現(xiàn)緩慢減小趨勢，1950年后流量分布均值增幅逐漸增大，尤其在1965年以后增幅顯著增大(見圖2)。結(jié)合文獻[8]的背景介紹，在83年觀測期間，該區(qū)域城市化進程主要分為兩個時期。以1965年為界，后半期的城市化進程明顯加劇。關于流量分布均方差(尺度參數(shù))的時變特征，與時間無關、線性相關、拋物線型相關均適用，自由度變化范圍較大，對模擬效果影響不顯著，不是模型的敏感參數(shù)。

圖1 模型I均值和均方差的時變特征

圖2 模型II位置參數(shù)和尺度參數(shù)的時變特征

(2)結(jié)合水文系列實測點，分析模型模擬效果。實測水文系列的點據(jù)分散，存在一定帶寬(見圖3)，通過單一頻率的擬合線僅能從趨勢上逼近，很難反映散點的分布特征。對不同時期、不同量級的實測流量賦予時間和頻率雙重屬性，圖4和圖5可反映了流量隨時間的變化趨勢和分布的頻率特征，更全面和直觀的解釋實測流量數(shù)據(jù)的意義。

圖3 模型擬合效果

圖4 模型I流量分位數(shù)

圖5 模型II流量分位數(shù)

(3)根據(jù)模型成果，分析水文系列在不同時期的頻率特征值。通過centiles.pred函數(shù)可導出不同時間各種頻率相應流量。根據(jù)文獻[8]成果，經(jīng)自由度優(yōu)選模型，該流域100年一遇洪水洪峰流量，在計算之初的1924年為2.8 m3/s，20世紀50年代后期最小值為2.1 m3/s，截至2007年達到5.1 m3/s。在本研究中，根據(jù)模型I和II模擬結(jié)果，100年一遇洪水洪峰流量設計值在1924年分別為2.5 m3/s和2.7 m3/s，截至2006年分別達5.1 m3/s和4.6 m3/s。其中，模型I在1947年～1950年降至最小值2.1 m3/s，模型II在1951年～1956年最小值為1.9 m3/s。多種方法計算成果基本一致，均有效揭示了該區(qū)域水文時間序列的非平穩(wěn)特征。

由于水文系列的非平穩(wěn)特征，在不同時期，不同頻率設計值波動較大；因此，頻率值的時間屬性尤為重要。

在模型I中，1940年～1958年期間，1000年一遇洪水洪峰流量約2.6 m3/s，該流量級相當于1976年的100年一遇、1987年的20年一遇、1992年的10年一遇；在模型II中，2006年10年一遇洪峰流量約3.3 m3/s，相當于1999年的20年一遇，1989年的100年一遇，1926年和1980年的1000年一遇。隨著時間變化，由于下墊面條件不同，某個時期10年一遇洪水量級甚至可能超過其他時期1000年一遇洪水。GAMLSS模型揭示的相同流量重現(xiàn)期差異正反映了水文時間序列的非平穩(wěn)現(xiàn)象。

5 結(jié) 語

本研究通過GAMLSS模型及RStudio開發(fā)平臺研究水文時間序列的非平穩(wěn)特征?；贕AM模型及GAMLSS模型原理，根據(jù)對水文時間序列模型參數(shù)化方案不同的處理方式，在前人研究的基礎上，按照模型I和模型II兩種方案，對水文系列進行模擬。借鑒美國Little Sugar Creek 1924年至2006年共83年歷年最大洪峰流量數(shù)據(jù)的研究案例，介紹了GAMLSS模型研究水文系列非平穩(wěn)性的方法，通過概率分布時變模式選擇、模型優(yōu)選、殘差判別、成果分析的系列流程，多模型多途徑對比分析模型的有效性。

研究發(fā)現(xiàn)，GAMLSS模型研究方法和流程可有效揭示水文系列的非平穩(wěn)性特征。在氣候變化和人類活動加劇的區(qū)域，水文系列的非平穩(wěn)性已顯著影響統(tǒng)計參數(shù)和頻率設計值。以研究案例區(qū)域為例，受城市化影響，某個時期10年一遇洪水量級甚至可能超過其他時段1000年一遇洪水，由此造成了工程設計和環(huán)境評估的極大不確定性，應在推求頻率設計值時，綜合考慮設計值的頻率要素和時間趨勢。