，， ， ， ，

(1. 長江水利委員會 水文局， 武漢 430010；2. 中國長江三峽集團有限公司，湖北 宜昌 443133；3. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072)

1 研究背景

枯水流量是水文學(xué)中一個重要的特征值，研究其重現(xiàn)期及設(shè)計值對流域水資源管理、水環(huán)境保護等具有重要意義。然而，目前關(guān)于枯水的研究內(nèi)容、方法及深度遠低于洪水[1-4]。氣候變化及人類活動影響下，水文時間序列存在非一致性已成為共識[5-8]，基于一致性假設(shè)的設(shè)計結(jié)果受到質(zhì)疑。因此，研究非一致性條件下枯水流量設(shè)計值具有重要的實際應(yīng)用價值。

有學(xué)者選取具有物理意義的氣象因子為協(xié)變量進行非一致性水文頻率分析[9-10]，但將其與EWT或ENE概念下的重現(xiàn)期相結(jié)合并用于研究水文極值事件的重現(xiàn)期或設(shè)計值比較少見。本文應(yīng)用該方法探討渭河流域枯水極值事件的設(shè)計值，以期為該區(qū)域枯水期的水資源管理提供更為可靠的信息參考。

2 研究區(qū)域及數(shù)據(jù)

2.1 研究區(qū)域簡介

渭河發(fā)源于甘肅省渭源縣鳥鼠山，是黃河最大的支流。干流總長818 km，流域面積134 800 km2，介于黃土高原東南部33°40′N —37°26′N，103°57′E—110°27′E之間。渭河流域地處半干旱半濕潤大陸性季風(fēng)氣候區(qū)，平均多年總降水量540 mm，年平均氣溫6～14 ℃，年潛在蒸散發(fā)量為660～1 600 mm。渭河平均多年天然徑流量為10×109m3，約占黃河的17%[15-16]。渭河流域水系及水文氣象站點如圖1所示。

圖1渭河流域水系及水文氣象站點
Fig.1Location,topography,hydro-meteorologicalstationsandriversystemsoftheWeiheRiverbasin

2.2 研究數(shù)據(jù)

本研究主要應(yīng)用如下3類數(shù)據(jù)。

(1)觀測水文數(shù)據(jù)：華縣水文站1954—2009年日平均流量數(shù)據(jù)來自黃河水利委員會水文局，本文選取其連續(xù)30 d最小流量序列作為枯水極值事件進行研究(圖2(a))。

(2)觀測氣象數(shù)據(jù)：氣溫和降水是與徑流密切相關(guān)的氣象變量，因此本文選取其作為協(xié)變量對枯水徑流進行非一致性頻率分析。渭河流域及其周邊22個氣象站點1954—2009年日平均氣溫及日總降水來自中國氣象科學(xué)數(shù)據(jù)共享網(wǎng)(http://cdc.cma.gov.cn)。2個變量華縣站以上的面平均序列由泰森多邊形法獲得，最終提取年平均氣溫(Temp)及年總降水(Prep)統(tǒng)計量作為枯水流量頻率分析的協(xié)變量。2個統(tǒng)計量與枯水極值序列存在一定的線性相關(guān)(圖2(b)和圖2(c))。

(3)NCEP(National Centres for Environmental Prediction)再分析數(shù)據(jù)以及CMIP5(Coupled Model Intercomparison Project Phase 5)的GCM(General Circulation Model)輸出數(shù)據(jù)：本文應(yīng)用NCEP再分析數(shù)據(jù)及GCM輸出數(shù)據(jù)結(jié)合統(tǒng)計降尺度方法獲取氣溫及降水的未來情景[17]。由于NCEP和GCM數(shù)據(jù)均為網(wǎng)格數(shù)據(jù)且空間分辨率不同，本文先將2組數(shù)據(jù)插值到氣象站點再應(yīng)用泰森多邊形法獲取各預(yù)報因子的面平均序列。

3 研究方法

3.1 枯水極值事件超過制概率

隨機變量Z表示枯水極值事件。由于研究對象為枯水，即水資源的短缺特點，因此定義設(shè)計枯水流量zp0為在任一年里“小于”該流量的事件出現(xiàn)概率為p0。洪水量級高于某一設(shè)計洪水的概率通常被稱為超過概率[18]，本文仍采用“超過概率”一詞，但此時“超過”表示枯水流量“小于”某一設(shè)計值。

pt=Fz(zp0,θt),t=1,2,...,T。

(1)

圖2華縣水文站枯水流量序列與協(xié)變量相關(guān)圖
Fig.2Correlationanalysisbetweenthelow-flowseriesandcorrespondingcovariatesofHuaxianstation

當(dāng)θt=θ0，即假設(shè)一致性條件仍成立，式(1)將最終退化為pt=p0。

3.2 ENE概念下非一致性設(shè)計枯水流量

在非一致性條件下，M的期望可表示為[14]

(2)

根據(jù)ENE定義，E(M)=1，通過數(shù)值迭代式(3)可得非一致性條件下設(shè)計值zP0的重現(xiàn)期T，即

(3)

通過對實測枯水序列進行非一致性頻率分析，進一步求出式(1)中未來年份的θt(當(dāng)選取氣象因子為協(xié)變量時，需要代入GCM輸出的氣象因子推求θt)，對于某一目標(biāo)重現(xiàn)期T，即可推求非一致性條件下的初始設(shè)計值zp0。

3.3 非一致性枯水流量頻率分析

針對某一目標(biāo)重現(xiàn)期，研究非一致性設(shè)計枯水流量的一個重要步驟為確定式(1)中未來年份的超過概率pt，而pt的確定又依賴于統(tǒng)計參數(shù)與協(xié)變量之間的函數(shù)關(guān)系。前文介紹的研究中僅選取時間為協(xié)變量進行非一致性頻率分析并進一步推求極值事件的重現(xiàn)期或設(shè)計值[11-14]，而氣象變量相比于時間具有更強的物理意義及更可信的解釋能力，因此本文將建立枯水流量統(tǒng)計參數(shù)與氣象因子的關(guān)系，并結(jié)合統(tǒng)計降尺度下的GCM輸出氣象數(shù)據(jù)用以計算未來年份超過概率pt，最終推求目標(biāo)重現(xiàn)期T所對應(yīng)的非一致性設(shè)計枯水流量zP0。

考慮位置、尺度和形狀參數(shù)的廣義可加模型(Generalized Additive Models for Location，Scale and Shape, GAMLSS)[19]結(jié)合時變矩法對實測枯水流量進行非一致性頻率分析，選取Gamma，Weibull，Gumbel，Logistic，Lognormal 5種常用于擬合枯水流量的概率分布函數(shù)作為備選分布。同時檢查位置參數(shù)μ及尺度參數(shù)σ的非一致性并應(yīng)用AIC評價準(zhǔn)則[20]，選取最優(yōu)非一致性模型。分別用worm圖[21]、分位圖以及Filliben相關(guān)系數(shù)(Fr)[22]、Kolmogorov-Smirnov(KS)檢驗統(tǒng)計量(DKS)[23]定性及定量評價模型擬合優(yōu)度。

3.4 統(tǒng)計降尺度模型(SDSM)

應(yīng)用SDSM結(jié)合站點實測氣象數(shù)據(jù)、NCEP及GCM數(shù)據(jù)對渭河流域的氣溫和降水進行降尺度。參考Wilby等[17,24]的研究，選取平均海平面氣壓、500 hPa位勢高度場、500 hPa緯向風(fēng)速場、850 hPa氣溫以及近地表氣溫作為氣溫的預(yù)報因子，在此基礎(chǔ)上增加850 hPa相對濕度作為降水的預(yù)報因子。選取Nash-Sutcliffe效率系數(shù)(NSE)評價觀測期內(nèi)模擬效果。

4 結(jié)果與討論

4.1 非一致性枯水流量頻率分析

選取時間t為協(xié)變量對實測枯水流量進行頻率分析時，AIC評價準(zhǔn)則表明Weibull分布(μ和σ均為對數(shù)連接函數(shù))為最優(yōu)分布，且2個統(tǒng)計參數(shù)均為時間線性函數(shù)為最優(yōu)非一致性模型(圖3(a))。

選取年平均氣溫Temp及年總降水Prep為協(xié)變量對實測枯水流量進行頻率分析時，AIC評價準(zhǔn)則同樣表明Weibull分布(μ和σ均為對數(shù)連接函數(shù))為最優(yōu)分布，且2個統(tǒng)計參數(shù)分別為Temp和Prep的線性函數(shù)為最優(yōu)非一致性模型(圖3(b))。

圖3不同分布及非一致性模型擬合渭河流域枯水流量序列AIC值
Fig.3AICvaluesasfunctionsoftimeandmeteorologicalvariableswiththeobservedlow-flowseriesfittedbydifferentnonstationarymodelswithdifferentdistributions

以t為協(xié)變量的worm圖中，大部分點都落在95%置信區(qū)間內(nèi)(圖4(a))，同時大部分觀測點都在5%～95%分位曲線之間(圖4(b))，表明模型能夠較好地擬合出枯水流量的非一致性。Fr和DKS統(tǒng)計量同樣說明選取的最優(yōu)非一致性模型能夠通過檢驗(表1)。

Temp和Prep為協(xié)變量的worm圖中，所有點都落在95%置信區(qū)間內(nèi)(圖4(c))，同時大部分觀測點都在5%～95%分位曲線之間(圖4(d))，表明模型能夠充分地擬合出枯水流量的非一致性。Fr和DKS統(tǒng)計量同樣說明選取的最優(yōu)非一致性模型能夠通過檢驗(表1)。而且，以Temp和Prep為協(xié)變量的最優(yōu)模型相比于以t為協(xié)變量具有更小的AIC值(表1)及更優(yōu)的worm圖(圖4(a), 圖4(c))，進一步驗證了選取具有物理意義的氣象因子作為協(xié)變量的必要性和有效性。

圖4 最優(yōu)非一致性模型擬合優(yōu)度評價結(jié)果Fig.4 Diagnostic plots for assessing the performance ofthe optimal nonstationary model using time andmeteorological covariates

注：*統(tǒng)計參數(shù)與協(xié)變量間關(guān)系為ln(μt)=μa+μb(t+τ),ln(σt)=σa+σb(t+τ),τ為觀測期長度且τ=56；**統(tǒng)計參數(shù)與協(xié)變量間關(guān)系為ln(μt)=μa+μbTempt,ln(σt)=σa+σbPrept

4.2 降水及氣溫的統(tǒng)計降尺度

對于以Temp和Prep為協(xié)變量的最優(yōu)非一致性模型，需要提供這2個變量未來年份的預(yù)估值才能進一步計算統(tǒng)計參數(shù)μt和σt以及式(1)中的超過概率。觀測期內(nèi)模擬和實測Temp間的NSE=0.69，說明氣溫模擬結(jié)果較好。但是對于Prep，模擬和實測值間的NSE=0.45，效果一般。前人研究也存在同樣問題，即降水的降尺度結(jié)果相比于氣溫效果要差[24-25]，因此本文認(rèn)為該模擬結(jié)果尚可接受。最終得到7個GCM模型下2010—2099年Temp和Prep的預(yù)估情景(圖5)，進而通過已建的最優(yōu)非一致性模型推求未來時期極端枯水事件的統(tǒng)計分布，結(jié)合ENE概念計算枯水流量設(shè)計值[26]。

圖5不同GCM下2010—2099年氣象因子降尺度結(jié)果
Fig.5ProjectedmeteorologicalvariablesfromdifferentGCMsforthefutureperiod2010-2099

4.3 ENE概念下非一致性設(shè)計枯水流量

給定目標(biāo)重現(xiàn)期T，以時間為協(xié)變量時，未來時期t=1,2,...,T年內(nèi)時變統(tǒng)計參數(shù)μt和σt可由最優(yōu)非一致性模型(表1中的*)直接計算，進而迭代推求與T相應(yīng)的非一致性設(shè)計枯水流量zp0。以Temp和Prep為協(xié)變量時，將未來時期t=1,2,...,T年內(nèi)Temp和Prep的降尺度數(shù)據(jù)(7個模型的平均值)代入最優(yōu)非一致性模型(表1中的**)計算時變統(tǒng)計參數(shù)μt和σt，進而推求與T相應(yīng)的非一致性設(shè)計枯水流量zp0。本文同時給出了一致性假設(shè)成立情況下與T相應(yīng)的設(shè)計枯水流量zp0。不同目標(biāo)重現(xiàn)期T對應(yīng)的設(shè)計枯水流量zp0如圖6。

圖6不同目標(biāo)重現(xiàn)期T對應(yīng)的設(shè)計枯水流量zp0

Fig.6Nonstationarydesignlow-flowquantilecorrespondingtodifferentreturnperiods

由圖6可知，以時間為協(xié)變量時，對于任一目標(biāo)重現(xiàn)期T，非一致性設(shè)計枯水流量明顯低于一致性假設(shè)情況。意味著在考慮非一致性的情況下，枯水極值事件相比于一致性假設(shè)將更加嚴(yán)重。例如，當(dāng)T=10 a時，以時間為協(xié)變量的非一致性設(shè)計枯水流量為35.1 m3/s，遠低于基于一致性假設(shè)的設(shè)計枯水流量175.3 m3/s。

由圖6可知，以Temp和Prep為協(xié)變量時，非一致性設(shè)計枯水流量同樣與一致性假設(shè)情況存在差別，但又不同于以時間為協(xié)變量的情況。當(dāng)目標(biāo)重現(xiàn)期T≤25 a時，非一致性設(shè)計枯水流量高于一致性假設(shè)情況，而當(dāng)T>25 a時，非一致性設(shè)計枯水流量逐漸低于一致性假設(shè)情況。意味著在考慮非一致性的情況下，當(dāng)目標(biāo)重現(xiàn)期較短時，枯水極值事件嚴(yán)重程度相比于一致性假設(shè)有所減輕；但當(dāng)目標(biāo)重現(xiàn)期較長時，枯水極值事件仍有所加重。例如，當(dāng)T=10 a，以Temp和Prep為協(xié)變量的非一致性設(shè)計枯水流量為232.8 m3/s，高于基于一致性假設(shè)的設(shè)計枯水流量175.3 m3/s。當(dāng)T=50 a，兩者分別為24.2 m3/s和48.0 m3/s。

4.4 討 論

無論以時間或氣象因子為協(xié)變量，非一致性設(shè)計枯水流量相比于一致性假設(shè)情況均有明顯差別(見圖6)，說明考慮枯水流量的非一致性很有必要。同時，2種不同協(xié)變量情況下的非一致性設(shè)計枯水流量也存在顯著不同。以時間為協(xié)變量所得設(shè)計枯水流量表明，渭河流域未來的枯水極值事件將更加嚴(yán)重。然而，以氣溫和降水為協(xié)變量所得設(shè)計枯水流量顯示，未來枯水極值事件變緩或者加重取決于目標(biāo)重現(xiàn)期T的長短。

以時間為協(xié)變量時，非一致性枯水頻率分析所擬合的統(tǒng)計參數(shù)μ和σ隨時間單調(diào)遞減(圖7)，如前文所述，該情況既缺少物理意義又無法擬合統(tǒng)計參數(shù)實際變化較為復(fù)雜的情形，同時假設(shè)擬合出的單調(diào)下降趨勢在未來無限持續(xù)也存在一定的不合理性。從圖2(a)也可發(fā)現(xiàn)，1995年左右，枯水流量序列有一個緩慢的上升，同樣佐證了統(tǒng)計參數(shù)隨時間持續(xù)下降假設(shè)的不合理性。相反地，以氣象因子為協(xié)變量的非一致性模型擁有更好的模型性能(見圖3、圖 4和表1)、更合理的擬合參數(shù)(圖7)以及更可信的枯水流量設(shè)計值。同時，T=50 a所對應(yīng)的2種情況下的非一致性設(shè)計枯水流量分別為0.4 m3/s和24.2 m3/s，而在本文56 a觀測期內(nèi)最小實測枯水流量為25.6 m3/s，同樣說明了本文所提方法的合理性和有效性。

圖7位置和尺度參數(shù)擬合結(jié)果
Fig.7Resultsofthefittedlocationandscaleparameters

5 結(jié) 論

在期望超過次數(shù)(ENE)概念下，將氣象因子引入到枯水流量非一致性頻率分析中并最終推求對應(yīng)某一目標(biāo)重現(xiàn)期的非一致性設(shè)計枯水流量。該方法與傳統(tǒng)的以時間為協(xié)變量的方法一起應(yīng)用于渭河流域設(shè)計枯水流量研究，所得結(jié)論如下：

(1)Weibull分布(2個統(tǒng)計參數(shù)均為對數(shù)連接函數(shù))為模擬實測枯水流量的最優(yōu)分布，且其位置和尺度參數(shù)均存在顯著的非一致性。以時間為協(xié)變量時，最優(yōu)模型中2個參數(shù)均為時間的線性函數(shù)；以氣象因子為協(xié)變量時，最優(yōu)模型中2個參數(shù)分別為年平均氣溫和年總降水的線性函數(shù)。并且，以氣象因子為協(xié)變量的最優(yōu)非一致性模型優(yōu)于以時間為協(xié)變量的情況。

(2)以時間為協(xié)變量和以氣象因子為協(xié)變量2種情況下的非一致性設(shè)計枯水流量相比于一致性假設(shè)情況存在顯著差異。同時，2種非一致性設(shè)計枯水流量之間也存在著明顯不同，前者由于不太合理的參數(shù)擬合及假設(shè)使得設(shè)計結(jié)果有所失準(zhǔn)，而后者則提供了更具物理意義、更為可信的設(shè)計結(jié)果。

(3)本文應(yīng)用GCM的降尺度數(shù)據(jù)時仍存在較大的不確定性，其中降水較明顯。同時，僅選取了RCP8.5情景下的GCM數(shù)據(jù)，將來工作中將考慮其他不同典型濃度情景以及不同GCM。最后，非一致性設(shè)計枯水流量的不確定性將是下一步工作研究重點。

(編輯：姜小蘭)