李 敏，李建柱，馮 平，陳 亮

（天津大學(xué) 水利工程仿真與安全國家重點實驗室，天津 300350）

1 研究背景

20世紀(jì)以來，干旱成為了亞洲地區(qū)影響力最廣泛的自然災(zāi)害之一，超過6百萬的人死于干旱，11億人受到影響，并能夠引起糧食危機、生態(tài)環(huán)境惡化等［1-3］。美國氣象學(xué)會將干旱分為4種類型；氣象干旱、水文干旱、農(nóng)業(yè)干旱與社會經(jīng)濟干旱［4］。水文干旱是自然干旱中的一種重要類型，關(guān)系著水文循環(huán)和水量平衡，并且對人類用水有著直接的影響。為了分析水文干旱特征，國內(nèi)外學(xué)者提出了多種干旱評估指數(shù)，如標(biāo)準(zhǔn)化徑流指數(shù)SRI（standardized runoff index）、累計流量距平、地表供水指數(shù)SWSI（surface water supply index）、河川徑流干旱指數(shù)SDI（Streamflow Drought Index）、標(biāo)準(zhǔn)化徑流指數(shù)SSI（Standardized Streamflow Index）等［5-10］。通過計算干旱評估指數(shù)，可以客觀地確定干旱歷時、干旱烈度與干旱頻率，進而準(zhǔn)確地反映干旱的特征及變化規(guī)律。

Milly等［11］指出，在氣候變化和人類活動的共同作用下，基于一致性假設(shè)的水文頻率計算理論和方法將面臨失真。因此，研究變化環(huán)境下的水文干旱演變規(guī)律，成為目前干旱領(lǐng)域研究的熱點與難點。Russo等［12］提出了基于GAMLSS模型的非一致性干旱指數(shù)SnsPI（Stationary and non-stationary Standardized Precipitation Index），該指數(shù)考慮了變化環(huán)境下降水序列服從的Gamma分布的參數(shù)μ為時間的函數(shù)，并應(yīng)用該指數(shù)分析了變化環(huán)境下歐洲極端干濕狀況。應(yīng)用結(jié)果也表明，SnsPI的計算原理更合理，能夠較好地反映實際干濕狀況。Li等［13］建立了以氣候因子為協(xié)變量的非一致性模型，提出了非一致性標(biāo)準(zhǔn)化降水指數(shù)（NSPI），結(jié)果表明，NSPI指數(shù)比傳統(tǒng)的指數(shù)SPI更好的反映出灤河流域的干旱特征。Wang等［14］假設(shè)灤河流域夏季降水序列服從參數(shù)μ隨時間變化的Gamma分布，計算了基于GAMLSS模型的SPIt指數(shù)。結(jié)果表明以時間為協(xié)變量的指數(shù)SPIt對灤河流域的氣象干旱評價效果優(yōu)于傳統(tǒng)的SPI指數(shù)。

目前針對水資源、洪水頻率進行非一致性的研究較多［15-18］，而對變化環(huán)境條件下形成的非一致性干旱特征分析的研究較少。Shukla等［5］類比SPI指數(shù)［19］的理論建立了標(biāo)準(zhǔn)化徑流指數(shù)SRI來評價干旱。

SRI在實際應(yīng)用中具有和SPI相同的優(yōu)勢，該指數(shù)計算簡單，所需輸入數(shù)據(jù)單一且容易獲得，適用于多時間尺度計算，可以滿足不同地區(qū)、不同應(yīng)用的需求，能夠為不同時間尺度的干旱監(jiān)測服務(wù)。本文以SRI指數(shù)為基礎(chǔ)，提出了變化環(huán)境下時變標(biāo)準(zhǔn)化徑流指數(shù)SRIt，即假設(shè)徑流序列為二階非平穩(wěn)序列，序列服從Gamma分布的兩參數(shù)μ、σ均隨時間變化，充分考慮了徑流序列分布參數(shù)隨時間變化的可能性，以研究區(qū)域歷史干旱事件為研究對象，對比分析傳統(tǒng)SRI和SRIt的評價結(jié)果，以驗證該指數(shù)的合理性。

2 研究區(qū)域概況

潘家口水庫流域位于河北省唐山市與承德地區(qū)交界處，面積約為3.37萬km2，占灤河流域總面積的75%。流域?qū)俑睙釒Ъ撅L(fēng)區(qū)，春季少雨，夏季多雨，冬季雨雪稀少，易旱災(zāi)。流域多年均氣溫為5～12℃，多年平均降水量在400～700 mm之間，年內(nèi)分配不均，主要集中在夏季，夏季降水量在200～560 mm之間。年水面蒸發(fā)量在950～1150 mm之間。流域徑流年際豐枯變化比較懸殊，常出現(xiàn)連續(xù)豐水年和枯水年，多年平均徑流量為24.5億m3。根據(jù)歷史資料記載［20］，在 1949—2010年期間，灤河流域的主要干旱災(zāi)害發(fā)生于1961、1963、1968、1972、1980—1984、2000、2007及2009年。近年來，由于全球氣候變化和流域人類活動的影響日益加劇，灤河流域降水量和徑流量均明顯趨于減少，干旱現(xiàn)象頻發(fā)，特別是2000年以后，流域極端干旱和連年干旱頻繁出現(xiàn)。干旱事件的頻發(fā)造成流域大面積受災(zāi)，損失嚴(yán)重，制約了流域經(jīng)濟和社會的發(fā)展。研究區(qū)域位置及區(qū)域內(nèi)各水文站點分布如圖1。

圖1 潘家口水庫流域地理位置及9個水文站分布

3 研究方法

3.1 標(biāo)準(zhǔn)化徑流指數(shù)SRI 標(biāo)準(zhǔn)化徑流指數(shù)SRI是由Shukla和Wood［5］于2008年首次提出的水文干旱指數(shù)，SRI的計算就是將給定時間尺度的累積徑流量的分布通過等概率變換轉(zhuǎn)化為標(biāo)準(zhǔn)正態(tài)分布。一定時間尺度的累計徑流量的計算公式如下：

式中：v表示年份；τ表示月份；k為給定的時間尺度；為給定時間尺度k下的累積徑流量，m3；x為τ-i月的月徑流量，m3。本文選用兩參數(shù)的Gamma分布作為一定時間尺度的累積徑流量的分布線型。兩參數(shù)Gamma分布概率密度函數(shù)的表達式為：

時變標(biāo)準(zhǔn)化徑流指數(shù)SRIt本文基于GAMLSS模型構(gòu)建以時間為協(xié)變量的非一致性模型。在假設(shè)水文序列所服從的分布線型不變時，該模型能夠靈活地模擬概率分布的參數(shù)隨協(xié)變量的變化［21］。模型的參數(shù)估計方法采用極大似然法，模型的評價準(zhǔn)則采用全局?jǐn)M合偏差GD（Global Deviance，GD）、AIC（Akaike Information Criterion）準(zhǔn)則［22］和 SBC（Schwarz Bayesian Criterion）［23］準(zhǔn)則，判斷準(zhǔn)則指數(shù)越小，模型模擬效果越好。為了避免回歸方程太過復(fù)雜，本文將多項式次數(shù)的上限設(shè)定為二次，即n范圍為0～2。

在時變矩模型中，假設(shè)實測徑流序列的分布參數(shù)是時間t的函數(shù)，即：

本文采用n次多項式來描述分布參數(shù)μ、σ與時間t之間的關(guān)系，表示為：

式中，a0、b0為常數(shù)項，ai、bi為多項式系數(shù)，i=0，1，2，…，n。

為了比較兩參數(shù)隨時間的變化，本文分別建立兩參數(shù)μ、σ均為常數(shù)，即均不隨時間變化、單個參數(shù)μ隨時間變化、兩參數(shù)μ、σ均隨時間變化的三種模型（i=0，1，2），并且基于三種模型采用全局?jǐn)M合偏差GD、AIC準(zhǔn)則和SBC準(zhǔn)則分別選出擬合效果最好的最優(yōu)模型，分別記為模型1、模型2與模型3。

本文通過計算模型殘差的均值、方差、偏態(tài)系數(shù)、峰態(tài)系數(shù)和Filliben系數(shù)以及Worm圖來評價最優(yōu)模型的構(gòu)建是否合理。若最優(yōu)模型殘差序列的均值接近于0，方差接近于1，偏態(tài)系數(shù)接近于0，峰態(tài)系數(shù)接近于3，F(xiàn)illiben系數(shù)大于等于標(biāo)準(zhǔn)值，則認(rèn)為殘差序列近似服從標(biāo)準(zhǔn)正態(tài)分布，最優(yōu)模型構(gòu)建合理。

表1 SRI與SRIt分級標(biāo)準(zhǔn)及相應(yīng)閾值

由于標(biāo)準(zhǔn)化徑流指數(shù)SRI的計算類似于SPI指數(shù)，因此，選用與SPI相同的干旱等級劃分標(biāo)準(zhǔn)［13］。指數(shù)的正值表示濕潤狀況，負(fù)值表示干旱狀況。

4 結(jié)果與分析

4.1 徑流序列非一致性識別 本文采用Mann-Kendall［24］檢驗法，對研究區(qū)域9個水文站1960—2011年的12個月時間尺度的累計徑流量序列進行趨勢分析，分析結(jié)果如表2。從表2得出，研究區(qū)域9個水文站的徑流序列均呈現(xiàn)下降趨勢，且除溝臺子站與韓家營站以外，其余站點的的徑流序列呈顯著下降趨勢。

表2 潘家口水庫流域9個水文站徑流序列趨勢檢驗結(jié)果

4.2 模型評價 以潘家口水庫流域波羅諾站為代表，以逐月實測徑流數(shù)據(jù)為基礎(chǔ)分別建立一致性模型1、以單參數(shù)μ隨時間變化的非一致性模型2、以兩參數(shù)μ、σ均隨時間變化的非一致性模型3，并分別通過GD、AIC和SBC準(zhǔn)則確定各自的最優(yōu)模型，由于徑流序列一般具有年尺度的周期性，具有代表和實踐意義，因此，本文以12個月時間尺度展開分析。波羅諾站12個月時間尺度的逐月徑流系列的3種模型的GD、AIC、SBC見表3。由表3可知，對于波羅諾站12個月尺度的逐月累計徑流分析序列，模型3的3個指數(shù)整體上最小，擬合效果最優(yōu)，其次為模型2。

表3 波羅諾站12個月時間尺度的逐月徑流系列模型的GD、AIC、SBC

表4 波羅諾站逐月的12個月時間尺度最優(yōu)模型參數(shù)

表4列出了波羅諾站逐月的12個月時間尺度最優(yōu)模型參數(shù)，從表4中看出，每個月的參數(shù)均與時間t有關(guān)。對于所采用的徑流樣本序列長度為52，當(dāng)Filliben系數(shù)≥0.978時才通過顯著性水平為95%的檢驗［25］。波羅諾站逐月的的12個月時間尺度的最優(yōu)模型殘差的均值、方差、偏態(tài)系數(shù)、峰態(tài)系數(shù)和Filliben系數(shù)均符合標(biāo)準(zhǔn)，認(rèn)為各站的最優(yōu)模型構(gòu)建合理。

圖2展示了波羅諾站1960—2011年累計徑流序列的殘差worm圖及分位數(shù)灰度圖，由于篇幅的限制，以12月份的12個月時間尺度為例說明。圖2（a）為累計徑流序列的最優(yōu)模型（模型3）下的殘差worm圖，所有殘差點均處在兩條橢圓曲線包圍的“可接受”區(qū)內(nèi)，表明Gamma分布對徑流序列的擬合效果較好。由圖2（b）—（d）中的分位數(shù)灰度圖可以看出，實測點據(jù)在最優(yōu)時變矩模型（模型3）中的分位數(shù)灰度圖中的分布狀況明顯更合理。模型1在0～50%分位數(shù)區(qū)間覆蓋了53%的實測點，25%～75%分位數(shù)灰度區(qū)間覆蓋了73%的實測點據(jù)；模型2在0～50%分位數(shù)區(qū)間覆蓋了54%的實測點，25%～75%分位數(shù)灰度區(qū)間覆蓋了76%的實測點據(jù)；而模型3在0～50%分位數(shù)區(qū)間覆蓋了58%的實測點，25%～75%分位數(shù)灰度區(qū)間覆蓋了77%的實測點據(jù)。因此，模型3對徑流實測點據(jù)的擬合效果更好，較好地描述了徑流數(shù)據(jù)的非線性隨機變化。此外，從模型3可以看出，波羅諾站徑流量序列整體上呈下降趨勢，與Mann-Kendall趨勢分析結(jié)果一致，且約在1960—1980年為增加趨勢，在1981—2011年為下降趨勢，1980年左右為徑流量序列的可能變異點，尤其是上淺灰色的變化趨勢尤為明顯，這與Li［26］，Wang［27］分析的結(jié)果一致，也證明了模型3具有較好的擬合效果。

圖2 波羅諾水文站1960—2011年12月份12個時間尺度的累計徑流序列模型3的Worm

通過比較模型1、模型2、模型3的GD、AIC、BIC值，選出3種模型中GD、AIC、BIC值最小的模型作為研究區(qū)域各水文站的最優(yōu)模型。表5展示了潘家口水庫流域各水文站12個月時間尺度的累計徑流量的最優(yōu)模型參數(shù)，考慮到篇幅的限制，只列出12月份最優(yōu)模型參數(shù)。可以看出，溝臺子站與韓家營站徑流量序列的Gamma分布的參數(shù)μ、σ均為常數(shù)，即最優(yōu)模型為模型1，該結(jié)果與Mann-Kendall趨勢檢驗的結(jié)果一致，認(rèn)為兩個站的徑流序列非一致性特征不明顯；三道河子站與波羅諾站的參數(shù)μ、σ均表現(xiàn)出與時間t的相依關(guān)系，即最優(yōu)模型均為模型3；其余站點的參數(shù)μ與時間t存在相依關(guān)系，σ為常數(shù)。此外，表5中各水文站的最優(yōu)模型殘差的均值、方差、偏態(tài)系數(shù)、峰態(tài)系數(shù)和Filliben系數(shù)均符合標(biāo)準(zhǔn)，認(rèn)為各站的最優(yōu)模型構(gòu)建合理。

4.3 指數(shù)的驗證 以歷史干旱事件為研究對象，對比分析12個月時間尺度的SRI1（模型1）、SRI2（模型2）、SRI3（模型3）的評價結(jié)果，從而驗證了指數(shù)的合理性和適用性。SRI1、SRI2、SRI3的計算結(jié)果見圖3。從圖中可以看出，傳統(tǒng)指數(shù)SRI1與非一致性指數(shù)SRI2、SRI3的波動變化與歷史記載的干旱年份均較吻合。

1972年的旱災(zāi)為灤河流域歷史上罕見的大旱，流域內(nèi)多條河流斷流，灤河年徑流量僅為多年平均值的44%，導(dǎo)致區(qū)域水資源嚴(yán)重短缺，受災(zāi)面積高達1907 km2，成災(zāi)1567 km2，旱死禾苗為173 km2，干旱面積占耕地面積的80%［28］。在1972年6月至1973年7月期間的14個月內(nèi)，根據(jù)12個月時間尺度的SRI1指數(shù)計算結(jié)果，均屬于中度干旱；而在這期間，根據(jù)12個月時間尺度的SRI2的計算結(jié)果，1972年12月與1973年6月均為重旱，值分別為-1.9、-1.5，且1972年12月接近于極度干旱，分別反映出了1972年全年與1972年6月至1973年6月的干旱情況，與實際較為接近；根據(jù)12個月時間尺度的SRI3的計算結(jié)果，1972年12月SRI3值為-2.1，為特旱，反映出1972年全年的干旱情況，與實際情況更加接近。

表5 潘家口水庫流域各水文站12月份的12個月時間尺度的累計徑流量的最優(yōu)模型參數(shù)

圖3 典型站波羅諾站1960—2011年12個月時間尺度的逐月標(biāo)準(zhǔn)化徑流指數(shù)

1980—1984年間，灤河流域各地均出現(xiàn)長時期降水不足的現(xiàn)象，致使干旱災(zāi)害大范圍連續(xù)發(fā)生，旱情較為嚴(yán)重，給流域農(nóng)業(yè)生產(chǎn)和民居正常生活都帶來嚴(yán)重影響［20］。在這一時期，12個月時間尺度的SRI1計算的結(jié)果在1981年10月至1982年6月為重度干旱，其余時段均為中度干旱；12個月時間尺度的SRI2計算的結(jié)果在1980年12月，1981年9、10月，1982年6、9月，1983年9、11月，1984年7、9、10、12月均為重旱；SRI3計算的結(jié)果在1981年10、12月，1983年9月，1984年7、10、12月均為重旱。SRI2與SRI3計算的結(jié)果能夠反映出整個時期的干旱較為嚴(yán)重，與實際更加接近。

4.4 干旱特征分析 本文分別以波羅諾水文站12個月時間尺度的SRI1、SRI2、SRI3的干旱等級評價結(jié)果為基礎(chǔ)，將整個研究期劃分為1960—1972、1973—1985、1986—1998、1999—2011年時間長度比較均勻的4個階段，排除時間長度不同造成的干擾，分別計算4個時間段的各級水文干旱等級頻率，從而比較SRI1、SRI2、SRI3在不同時間段的干旱等級頻率的變化過程。計算結(jié)果見圖4。

圖4 波羅諾站12個月時間尺度的SRI1、SRI2與SRI3在4個時間段的水文干旱等級頻率

由圖4可知，（a）時段，在輕微干旱至極度干旱情況下（等級5～8級），SRI1、SRI2與SRI3對應(yīng)的頻率分別為48.08%、87.82%、84.62%，可以看出，SRI2與SRI3均比SRI1的頻率高，而在極端濕潤至正常情況下（1～4級），結(jié)果則相反。說明SRI2與SRI3計算的干旱程度比SRI1嚴(yán)重；在（b）時段，頻率的變化規(guī)律與前一個時間段（a）基本相同；在（c）時段，基于SRI1的等級主要集中于極度濕潤至輕微干旱的情況（1～5級），而基于SRI2與SRI3的干旱等級主要發(fā)生于正常至中等干旱的情況（4～6級）；在（d）時段，基于SRI1的干旱等級全部分布在輕微至中等干旱（5～6級），而SRI2與SRI3的干旱等級主要發(fā)生于輕微至重度干旱（5～8級）的情況，對應(yīng)的頻率分別為91.84%、89.8%。（c）與（d）時段也表現(xiàn)為SRI2與SRI3的干旱程度比SRI計算的更加嚴(yán)重。

總結(jié)以上分析，4個時間段在整體上均表現(xiàn)為SRI3與SRI2計算得到的干旱程度比SRI1嚴(yán)重，并且隨著時間的發(fā)展，3個指數(shù)計算得到的干旱情況在整體上均趨向于加重，主要原因是在非一致性模型中，以時間t為協(xié)變量能夠反映出變化環(huán)境下徑流序列變化的趨勢，該趨勢是影響徑流的一系列因素共同導(dǎo)致的，如氣候變化，人類活動等。而在受這些因素的影響下，潘家口水庫流域的水文干旱事件頻率增加，嚴(yán)重程度增強。

圖5展現(xiàn)出了波羅諾站12個月時間尺度的逐月累計徑流序列對應(yīng)的干旱等級。可以看出，使用傳統(tǒng)的SRI1計算得到的干旱等級只與累計徑流量有關(guān)，一定范圍的徑流量對應(yīng)一定的干旱等級，而SRI2與SRI3計算得到的干旱等級不只與累計徑流量有關(guān)，還與時間有關(guān)，同一個徑流量值因為發(fā)生時間的不同，可以對應(yīng)幾個干旱等級。例如，25×106～68×106m3范圍內(nèi)的徑流量，根據(jù)SRI1指數(shù)，得到干旱等級均為輕度干旱，而根據(jù)SRI2、SRI3指數(shù)，計算得到的干旱等級均有5種，分別是無旱、輕旱、中旱、重旱、特旱。此外，3個指數(shù)SRI3、SRI2與SRI1在5～8級干旱中分布的干旱事件依次減少。以上分析證明了在變化環(huán)境下，研究區(qū)域的徑流序列發(fā)生了顯著的非一致性變化以及干旱指數(shù)構(gòu)建中考慮非一致性的必要性。

圖5 波羅諾站12個月時間尺度的逐月累計徑流序列干旱等級統(tǒng)計

5 結(jié)論

本文以傳統(tǒng)的標(biāo)準(zhǔn)化徑流指數(shù)SRI的計算原理為基礎(chǔ)，基于以時間為協(xié)變量的非一致性模型構(gòu)建了時變標(biāo)準(zhǔn)化徑流指數(shù)SRIt，得出以下結(jié)論：（1）在受氣候與人類活動影響較大的背景環(huán)境下，1960—2011年間，潘家口水庫流域的9個水文站，除溝臺子站與韓家營站以外，其余站點的徑流序列呈顯著下降趨勢，即存在顯著的非一致性特征。（2）基于全局?jǐn)M合偏差GD、AIC準(zhǔn)則和SBC準(zhǔn)則計算出最優(yōu)模型參數(shù)。結(jié)果表明，徑流序列非一致性特征不顯著時，對應(yīng)的最優(yōu)模型參數(shù)均為常數(shù)，而當(dāng)徑流序列存在非一致性特征時，對應(yīng)的最優(yōu)模型參數(shù)為時間t的函數(shù)。以波羅諾站的徑流序列（存在非一致性特征）為例，比較一致性模型（模型1），一個參數(shù)μ隨時間變化的非一致性模型（模型2），兩個參數(shù)μ、σ均隨時間變化的非一致性模型（模型3）。結(jié)果表明，模型3的GD、AIC和SBC值整體上最小，其次為模型2。通過分位數(shù)灰度圖可以看出，模型3更能捕捉到徑流序列的變化情況，模擬效果更好，從而驗證了在非一致情況下，假設(shè)Gamma分布參數(shù)隨時間變化的合理性。（3）以歷史干旱事件為研究對象，對比分析波羅諾站3個指數(shù)的評價結(jié)果。結(jié)果表明，在變化環(huán)境下，SRI3最能反映出實際的干旱情況，驗證了時變標(biāo)準(zhǔn)化徑流指數(shù)構(gòu)建的合理性與必要性。（4）隨著時間的發(fā)展，SRI1與SRIt指數(shù)表示的干旱情況在整體上均趨向于加重，且在同一時間段內(nèi)，SRIt比SRI1表示的干旱情況整體上更加嚴(yán)重，這與Mann-Kendall非一致分析中得出的徑流序列存在顯著下降趨勢相對應(yīng)。此外，在時變矩模型中，同一范圍內(nèi)的徑流量隨時間的變化可能對應(yīng)不同的干旱等級。

為了降低非一致性模型計算的復(fù)雜程度，本文選擇較為簡單的Gamma分布作為徑流序列的理論分布線型。然而，模型的分布函數(shù)有多種選擇，如Lognormal分布、Pearson-Ⅲ型分布等，后續(xù)研究可通過模擬優(yōu)選出最優(yōu)的分布函數(shù)，進一步提高SRIt指數(shù)的適用性。