杜 濤，熊立華，李 帥，邵 駿，許崇育，4，閆 磊

（1.長(zhǎng)江水利委員會(huì)水文局，湖北 武漢 430010；2.武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072；3.中國(guó)長(zhǎng)江三峽集團(tuán)有限公司，湖北 宜昌 443133；4.奧斯陸大學(xué) 地學(xué)系，挪威 奧斯陸 N-0316）

1 研究背景

設(shè)計(jì)洪水是水利工程規(guī)劃設(shè)計(jì)、防洪決策等工作的一個(gè)重要特征值［1-2］。由于氣候變化及人類活動(dòng)影響，水文極值事件存在非一致性已成為共識(shí)［3-6］，致使基于一致性假設(shè)所得的設(shè)計(jì)結(jié)果受到質(zhì)疑。因此，研究非一致性條件下設(shè)計(jì)洪水具有重要的實(shí)際應(yīng)用價(jià)值。

目前，關(guān)于水文極值序列非一致性的研究主要集中于時(shí)變矩法［7-8］，其基本思想是：假設(shè)水文變量所服從的分布類型不變，但分布的統(tǒng)計(jì)參數(shù)隨時(shí)間或其他氣象協(xié)變量變化。在此情況下，給定某一超過(guò)概率p，不難得出每年各自的設(shè)計(jì)值，即，其中為水文極值變量Z的逆分布函數(shù)，θt為第t年的統(tǒng)計(jì)參數(shù)向量。Strupczewski等［9-11］首次提出一個(gè)非一致性洪水頻率分析模型，將趨勢(shì)性成分引入到洪水頻率分布的一、二階矩中；Katz等［12］在前人基礎(chǔ)上直接建立了洪水頻率分布參數(shù)與時(shí)間的關(guān)系；Villarini等［13］基于考慮位置、尺度和形狀的廣義可加模型（Generalized additive Models for Location，Scale and Shape， GAMLSS）研究了Little Sugar Creek流域年最大洪水趨勢(shì)性變化問(wèn)題，結(jié)果發(fā)現(xiàn)傳統(tǒng)意義下100年一遇的設(shè)計(jì)洪水在非一致性條件下已不再是一個(gè)固定值，該設(shè)計(jì)值最小為1957年的2.1 m3s-1km-2，最大為2007年的5.1 m3s-1km-2。理論上，時(shí)變矩方法能夠較好地描述洪水極值序列年際間的非一致性變化情況，然而在實(shí)際應(yīng)用中，給定某一超過(guò)概率，每年都有一個(gè)設(shè)計(jì)值較難應(yīng)用于實(shí)際［14］。Schumann［15］指出，考慮水文風(fēng)險(xiǎn)的水文極值事件設(shè)計(jì)值對(duì)水利工程規(guī)劃設(shè)計(jì)顯得越發(fā)重要。梁忠民等［16］基于“等可靠度”的概念建立了一致/非一致性條件下設(shè)計(jì)洪水計(jì)算方法間的聯(lián)系，保證了非一致性條件下水文設(shè)計(jì)成果與現(xiàn)行工程采用的成果之間的銜接與協(xié)調(diào)。Rootzén等［17］將水文風(fēng)險(xiǎn)概念引入到非一致性極值事件設(shè)計(jì)值的研究中。他們指出，在變化環(huán)境下水利工程規(guī)劃設(shè)計(jì)中，應(yīng)該量化兩方面的基礎(chǔ)信息：①設(shè)計(jì)使用年限，即水利工程將發(fā)揮作用的年限；②水文風(fēng)險(xiǎn)，即設(shè)計(jì)使用年限內(nèi)水文極值事件超過(guò)某一特定設(shè)計(jì)值的概率。在給定設(shè)計(jì)使用年限及水文風(fēng)險(xiǎn)情況下，將有唯一的設(shè)計(jì)值與之對(duì)應(yīng)，稱之為基于風(fēng)險(xiǎn)的水文極值事件設(shè)計(jì)值。Rootzén等［17］應(yīng)用以上方法研究了變化環(huán)境下非一致性年最大日降水及年最低氣溫的設(shè)計(jì)值。該方法為非一致性極值事件設(shè)計(jì)值的研究提供了一種基于風(fēng)險(xiǎn)的新思路，但其中存在有待改進(jìn)之處，即在應(yīng)用時(shí)變矩法對(duì)非一致性極值序列進(jìn)行頻率分析時(shí)僅選取時(shí)間因子作為協(xié)變量。單純以時(shí)間為協(xié)變量既無(wú)法擬合統(tǒng)計(jì)參數(shù)波動(dòng)性變化較為復(fù)雜的情況，同時(shí)又缺乏物理意義，更重要的是進(jìn)一步推求設(shè)計(jì)使用年限內(nèi)極值事件的統(tǒng)計(jì)分布時(shí)，假設(shè)所擬合出的上升、下降趨勢(shì)或者跳躍突變將在未來(lái)時(shí)期無(wú)限延續(xù)下去。

綜上，有研究選取具有物理意義的氣象因子為協(xié)變量進(jìn)行非一致性水文頻率分析［7-8］，同時(shí)也有研究以時(shí)間為協(xié)變量，將水文風(fēng)險(xiǎn)概念引入到非一致性極值事件設(shè)計(jì)值推求當(dāng)中［17］，但未見有研究將兩者相結(jié)合，以推求更具物理意義的非一致性設(shè)計(jì)洪水成果。本文擬采用該方法探討變化環(huán)境下渭河流域洪水極值事件的設(shè)計(jì)值，并用自助抽樣（Bootstrap）方法［18］推求該設(shè)計(jì)值的不確定性，以期為水利工程規(guī)劃設(shè)計(jì)及防洪決策提供一定參考依據(jù)。

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

渭河發(fā)源于甘肅省渭源縣鳥鼠山，是黃河最大支流。干流總長(zhǎng)818 km，流域面積134 800 km2，界于黃土高原東南部 33°40′—37°26′N，103°57′—110°27′E之間（圖 1）。渭河流域地處半干旱半濕潤(rùn)的大陸性季風(fēng)氣候區(qū)，多年平均降水量540 mm（1954—2009年，下同），年平均氣溫6～14℃，年潛在蒸散發(fā)量660～1600 mm，多年平均天然徑流量100億m3，約占黃河的17%［19-20］。

圖1 渭河流域水系及水文氣象站點(diǎn)

本研究主要應(yīng)用3類數(shù)據(jù)：（1）觀測(cè)水文數(shù)據(jù)；（2）觀測(cè)氣象數(shù)據(jù)；（3）美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心（National Centres for Environmental Prediction，NCEP）再分析數(shù)據(jù)以及國(guó)際耦合模式比較計(jì)劃第5階段（Coupled Model Intercomparison Project Phase 5， CMIP5）提供的大氣環(huán)流模型（General Circulation Mod?el，GCM）輸出數(shù)據(jù)。

（1）華縣水文站1954—2009年日平均流量數(shù)據(jù)來(lái)自黃委水文局，選取其年最大日平均流量序列作為洪水極值事件進(jìn)行研究（圖2（a））。華縣站位于渭河與黃河接口上游70 km處，站點(diǎn)控制面積106 500 km2，約為渭河流域總面積的80%。

（2）氣溫和降水是與徑流密切相關(guān)的氣象變量，因此本研究選取其作為協(xié)變量對(duì)洪水序列進(jìn)行非一致性頻率分析。渭河流域及其周邊22個(gè)氣象站點(diǎn)1954—2009年日平均氣溫及日總降水來(lái)自中國(guó)氣象科學(xué)數(shù)據(jù)共享網(wǎng)（http：//cdc.cma.gov.cn）。兩個(gè)變量的面平均序列由泰森多邊形法獲得，進(jìn)一步提取年平均氣溫（Temp）及年總降水（Prep）統(tǒng)計(jì)量作為非一致性洪水頻率分析的協(xié)變量（圖2（b）（c））。理論上講，季節(jié)性氣象統(tǒng)計(jì)量應(yīng)更適合于描述洪水極值事件的非一致性，考慮到華縣站約有90%的實(shí)測(cè)洪水極值事件發(fā)生在7—10月，本研究同樣檢驗(yàn)了夏季平均氣溫（TempS）和夏季總降水（PrepS）與實(shí)測(cè)洪水序列的相關(guān)性（圖2（d）（e）），結(jié)果稍低于年值統(tǒng)計(jì)量。

（3）本文應(yīng)用NCEP再分析數(shù)據(jù)及CMIP5提供的GCM輸出數(shù)據(jù)結(jié)合統(tǒng)計(jì)降尺度方法［21］獲取氣溫及降水的未來(lái)情景。1954—2009年的NCEP預(yù)報(bào)因子來(lái)自NOAA地球系統(tǒng)研究實(shí)驗(yàn)室（ESRL）（http：//www.esrl.noaa.gov）［22-23］。GCM輸出數(shù)據(jù)選取IPCC第5次評(píng)估報(bào)告中RCP2.6、RCP4.5和RCP8.5三種排放情景下 9個(gè)模型（BCC-CSM1.1，BNU-ESM，CanESM2，CCSM4，CNRM-CM5，GFDL-ESM3M，HadGEM2-ES，MIROC-ESM-CHEM，NorESM1-M）2010—2099年與NCEP相應(yīng)的大尺度預(yù)報(bào)因子（http：//cmip-pcmdi.llnl.gov/cmip5）。由于NCEP和GCM數(shù)據(jù)均為網(wǎng)格數(shù)據(jù)且空間分辨率不同，本研究首先將兩組數(shù)據(jù)插值到氣象站點(diǎn)，再應(yīng)用泰森多邊形法獲取各預(yù)報(bào)因子的流域面平均序列。

圖2 華縣水文站實(shí)測(cè)洪水序列與氣象協(xié)變量相關(guān)關(guān)系

3 研究方法

3.1 水文風(fēng)險(xiǎn)及設(shè)計(jì)洪水量級(jí)隨機(jī)變量Z表示洪水極值事件。在一致性條件下，F(xiàn)Z(z，θ0)表示Z的累積概率分布函數(shù)，其中θ0為固定的統(tǒng)計(jì)參數(shù)向量。給定某一特定設(shè)計(jì)使用年限T，則該設(shè)計(jì)使用年限內(nèi)洪水極值事件超過(guò)某一設(shè)計(jì)值z(mì)T的概率，即一致性條件下zT所對(duì)應(yīng)的水文風(fēng)險(xiǎn)RT為［14，17］：

式中：LT為設(shè)計(jì)使用年限總年數(shù)；t1和tn分別為設(shè)計(jì)使用年限初年和末年。

此時(shí)，給定某一水文風(fēng)險(xiǎn)RT，可推求相應(yīng)的一致性設(shè)計(jì)洪水量級(jí)zT，即轉(zhuǎn)變式（1）為：

在非一致性條件下，F(xiàn)Z(z，θt)表示Z的累積概率分布函數(shù)，其中統(tǒng)計(jì)參數(shù)向量θt隨時(shí)間或其他物理因子變化。于是有非一致性情況下設(shè)計(jì)值z(mì)T所對(duì)應(yīng)的水文風(fēng)險(xiǎn)為：

此時(shí)，給定某一水文風(fēng)險(xiǎn)RT，相應(yīng)的非一致性設(shè)計(jì)洪水量級(jí)zT將無(wú)法像式（2）顯式表達(dá)，但可通過(guò)數(shù)值迭代求出其數(shù)值解。當(dāng)θt=θ0，即一致性條件仍成立，式（3）將退化為式（1）。

通過(guò)對(duì)實(shí)測(cè)洪水序列進(jìn)行非一致性頻率分析，進(jìn)一步求出式（3）中設(shè)計(jì)使用年限內(nèi)各年的統(tǒng)計(jì)參數(shù)θt。此時(shí)針對(duì)特定設(shè)計(jì)使用年限T可進(jìn)行兩方面研究：①對(duì)于某一量級(jí)的設(shè)計(jì)洪水zT，推求非一致性條件下的水文風(fēng)險(xiǎn)RT；②對(duì)于某一水文風(fēng)險(xiǎn)RT，推求非一致性條件下相應(yīng)的設(shè)計(jì)洪水量級(jí)zT。第①種情況易于求解，然而實(shí)際工程設(shè)計(jì)中通常更關(guān)心第②種情況的設(shè)計(jì)結(jié)果，因此本文將對(duì)給定水文風(fēng)險(xiǎn)條件下的非一致性設(shè)計(jì)洪水及其不確定性進(jìn)行研究。

3.2 非一致性洪水頻率分析針對(duì)特定設(shè)計(jì)使用年限T及相應(yīng)風(fēng)險(xiǎn)RT，研究非一致性設(shè)計(jì)洪水量級(jí)zT的一個(gè)重要步驟即確定式（3）中設(shè)計(jì)使用年限內(nèi)各年的統(tǒng)計(jì)參數(shù)θt，而θt的確定又依賴于其與協(xié)變量之間的函數(shù)關(guān)系，即非一致性洪水頻率分析。本文前述國(guó)內(nèi)外相關(guān)研究成果中通常僅選取時(shí)間因子為協(xié)變量進(jìn)行非一致性頻率分析并進(jìn)一步推求極值事件的設(shè)計(jì)值［14，17］，而氣象變量相比于時(shí)間具有更強(qiáng)的物理意義及解釋能力，因此本文將建立洪水概率分布統(tǒng)計(jì)參數(shù)與氣象因子的關(guān)系，并結(jié)合統(tǒng)計(jì)降尺度下GCM輸出的氣象數(shù)據(jù)獲得設(shè)計(jì)使用年限內(nèi)各年的統(tǒng)計(jì)參數(shù)θt，最終推求某一水文風(fēng)險(xiǎn)RT所對(duì)應(yīng)的非一致性設(shè)計(jì)洪水量級(jí)zT及其不確定性。本文同時(shí)選取了以時(shí)間為協(xié)變量的情況進(jìn)行比較研究。選取國(guó)內(nèi)外較為常用的兩參數(shù)Gumbel（GU）分布和三參數(shù)GEV分布、P-Ⅲ分布作為備選洪水頻率分布，詳見表1。由于GEV和P-Ⅲ分布的形狀參數(shù)非常敏感且較難估計(jì)，通常不考慮該參數(shù)的非一致性變化［12，24-26］，故本文同時(shí)檢查各備選分布位置參數(shù)μ及尺度參數(shù)σ的非一致性［27］，并應(yīng)用赤池信息準(zhǔn)則（Akaike Information Criterion， AIC）評(píng)價(jià)準(zhǔn)則［28］選取最優(yōu)非一致性模型。用worm圖［29］、分位圖、Filliben相關(guān)系數(shù)（Fr）［30］以及 Kolmogorov-Smirnov（KS）檢驗(yàn)統(tǒng)計(jì)量（DKS）［31］評(píng)價(jià)模型擬合優(yōu)度。

表1 非一致性洪水頻率分析備選概率分布

3.3 統(tǒng)計(jì)降尺度模型（SDSM）GCM輸出數(shù)據(jù)使得預(yù)估未來(lái)時(shí)期氣象因子成為可能，將其引入以氣象因子為協(xié)變量的最優(yōu)非一致性模型中便可獲取設(shè)計(jì)使用年限T內(nèi)各年的統(tǒng)計(jì)參數(shù)θt并進(jìn)一步計(jì)算與水文風(fēng)險(xiǎn)RT相應(yīng)的設(shè)計(jì)洪水量級(jí)zT。然而，較低的空間分辨率限制了GCM的應(yīng)用范圍，針對(duì)此問(wèn)題，Wilby等［21，32］提出了統(tǒng)計(jì)降尺度模型（SDSM）來(lái)建立大尺度預(yù)報(bào)因子與區(qū)域或站點(diǎn)尺度預(yù)報(bào)量間的關(guān)系，該模型融合了天氣發(fā)生器和多元線性回歸技術(shù)。本文應(yīng)用SDSM結(jié)合站點(diǎn)實(shí)測(cè)氣象數(shù)據(jù)、NCEP及GCM數(shù)據(jù)對(duì)渭河流域的氣溫和降水進(jìn)行降尺度。

綜上，本文基于風(fēng)險(xiǎn)的非一致性設(shè)計(jì)洪水及其不確定性研究流程如圖3所示。

圖3 基于風(fēng)險(xiǎn)的非一致性設(shè)計(jì)洪水及其不確定性研究流程

3.4 Bootstrap方法推求設(shè)計(jì)洪水不確定性本研究采用Bootstrap方法［18］推求由于樣本數(shù)據(jù)長(zhǎng)度有限所引起的設(shè)計(jì)洪水統(tǒng)計(jì)不確定性。

（1）一致性條件下，采用Bootstrap方法計(jì)算設(shè)計(jì)洪水不確定性具體步驟為：①由原始實(shí)測(cè)洪水序列計(jì)算一定水文風(fēng)險(xiǎn)RT下的設(shè)計(jì)洪水量級(jí)zT；②采用有放回方式從原始實(shí)測(cè)洪水序列中抽取與其等長(zhǎng)度的新樣本序列；③用新樣本序列重新計(jì)算風(fēng)險(xiǎn)RT下的設(shè)計(jì)洪水量級(jí)；④重復(fù)第②、③步N次，得到N個(gè)抽樣設(shè)計(jì)值，i=1，2，…，N ；⑤對(duì)N個(gè)抽樣設(shè)計(jì)值由小到大排序，zT的置信區(qū)間為。

（2）非一致性條件下，無(wú)法從原始樣本中直接抽樣，需要采用對(duì)一致性殘差進(jìn)行抽樣的方法來(lái)推求設(shè)計(jì)值的不確定性，具體步驟為：①由原始實(shí)測(cè)洪水序列計(jì)算一定水文風(fēng)險(xiǎn)RT下的設(shè)計(jì)洪水量級(jí)zT；②由非一致性洪水頻率分析得到的時(shí)變統(tǒng)計(jì)參數(shù)計(jì)算模型標(biāo)準(zhǔn)正態(tài)殘差③采用有放回方式從殘差序列ri，i=1，2，…，n中抽取與其等長(zhǎng)度的新殘差序列；④用新殘差序列計(jì)算洪水樣本序列，重新按照第①步計(jì)算風(fēng)險(xiǎn)RT下的設(shè)計(jì)洪水量級(jí)；⑤重復(fù)第③、④步N次，得到N個(gè)抽樣設(shè)計(jì)值；⑥對(duì) N個(gè)抽樣設(shè)計(jì)值由小到大排序，zT的100(1-α)%置信區(qū)間為。

4 實(shí)例分析

4.1 洪水序列非一致性識(shí)別及模型建立采用Mann-Kendall檢驗(yàn)法對(duì)實(shí)測(cè)洪水序列進(jìn)行初步非一致性檢驗(yàn)，結(jié)果表明序列存在明顯下降趨勢(shì)。GAMLSS模型擬合結(jié)果如圖4所示，結(jié)果表明，一致性情況下的最優(yōu)分布為P-Ⅲ分布；以時(shí)間t為協(xié)變量情況下的最優(yōu)分布同樣為P-Ⅲ分布，且位置參數(shù)為時(shí)間的線性函數(shù)為最優(yōu)非一致性模型。定性評(píng)價(jià)指標(biāo)（worm圖和分位圖）表明模型能夠通過(guò)顯著性水平為0.05的擬合優(yōu)度檢驗(yàn)（圖5（c）（d）），定量評(píng)價(jià)指標(biāo)同樣說(shuō)明選取的最優(yōu)非一致性模型能夠通過(guò)檢驗(yàn)（表2），然而如前文所述，以缺乏物理意義的時(shí)間因子為協(xié)變量，不能對(duì)洪水序列的波動(dòng)性給出較好的擬合（圖5（d））。

圖4 非一致性模型擬合結(jié)果比較

表2 渭河流域?qū)崪y(cè)洪水序列非一致性頻率分析結(jié)果及評(píng)價(jià)指標(biāo)統(tǒng)計(jì)量

以年平均氣溫Temp及年總降水Prep為協(xié)變量對(duì)實(shí)測(cè)洪水序列進(jìn)行頻率分析時(shí)，AIC評(píng)價(jià)準(zhǔn)則表明Gumbel分布為最優(yōu)分布，且μ和σ兩個(gè)統(tǒng)計(jì)參數(shù)分別為Temp、Prep以及Temp的線性函數(shù)為最優(yōu)非一致性模型。定性和定量評(píng)價(jià)結(jié)果同樣說(shuō)明所選取的最優(yōu)非一致性模型能夠通過(guò)擬合優(yōu)度檢驗(yàn)（圖5（e）（f）和表2）。而且，以Temp和Prep為協(xié)變量的最優(yōu)模型相比于以t為協(xié)變量具有更小的AIC值，同時(shí)其不僅能夠擬合出洪水序列的下降趨勢(shì)，還可以對(duì)其波動(dòng)性給出較好的擬合（圖5（f）），進(jìn)一步印證了選取具有物理意義的氣象因子作為協(xié)變量的必要性和有效性。本研究同樣分析了以季節(jié)性統(tǒng)計(jì)量為協(xié)變量的情況，結(jié)果表明，P-Ⅲ分布為最優(yōu)分布，且μ為TempS、PrepS的線性函數(shù)，σ和κ均為常數(shù)，相應(yīng)模型檢驗(yàn)結(jié)果見表2，最終發(fā)現(xiàn)該最優(yōu)非一致性模型效果遠(yuǎn)差于以Temp和Prep為協(xié)變量情況，甚至不及以時(shí)間t為協(xié)變量情況。

本文同樣給出了一致性情況下擬合優(yōu)度檢驗(yàn)結(jié)果，結(jié)果表明，模型的Filliben相關(guān)系數(shù)不能通過(guò)顯著性檢驗(yàn)，說(shuō)明渭河流域?qū)崪y(cè)洪水序列確實(shí)存在顯著的非一致性，若仍然以一致性假設(shè)為前提進(jìn)行工程規(guī)劃設(shè)計(jì)將增大事故風(fēng)險(xiǎn)。

圖5 渭河流域最優(yōu)非一致性模型擬合優(yōu)度檢驗(yàn)worm圖及分位曲線圖

依據(jù)最優(yōu)分布線型及參數(shù)變化情況，用1954—2000年資料重新擬合每種情況下的最優(yōu)非一致性模型，進(jìn)而將所得模型應(yīng)用于2001—2009年以驗(yàn)證模型效果，結(jié)果見圖6。結(jié)果表明，相比一致性模型和以時(shí)間t為協(xié)變量的非一致性模型，以氣象因子為協(xié)變量的非一致性模型具有更強(qiáng)的預(yù)測(cè)能力。

圖6 最優(yōu)洪水頻率分布模型驗(yàn)證

4.2 降水及氣溫的統(tǒng)計(jì)降尺度對(duì)于以氣象因子為協(xié)變量的最優(yōu)非一致性模型，需要提供協(xié)變量未來(lái)年份的預(yù)估值才能進(jìn)一步計(jì)算設(shè)計(jì)使用年限T內(nèi)各年的統(tǒng)計(jì)參數(shù)μt和σt以及特定水文風(fēng)險(xiǎn)RT下的非一致性設(shè)計(jì)洪水zT。參考Wilby等［21，32］，經(jīng)相關(guān)性分析，本文最終選取降水和氣溫的預(yù)報(bào)因子如表3所示。

表3 日降水和日平均氣溫降尺度所選大尺度預(yù)報(bào)因子及模型類型

本文實(shí)測(cè)日降水和日平均氣溫為1954—2009年共56年資料，在建立SDSM模型時(shí)，選取1954—1993年共40年的數(shù)據(jù)率定模型，1994—2009年共16年的數(shù)據(jù)驗(yàn)證模型。結(jié)果表明，所建立的SDSM模型在率定期和檢驗(yàn)期內(nèi)模擬值與觀測(cè)值擬合良好，結(jié)果見圖7。進(jìn)一步將3個(gè)典型濃度路徑（RCP2.6、RCP4.5和RCP8.5）下的9個(gè)GCM模式［33］輸出的大尺度氣候因子2010—2099年的未來(lái)預(yù)估值代入到所建立的SDSM模型當(dāng)中，生成每個(gè)GCM模式下2010—2099年日降水和日平均氣溫的未來(lái)預(yù)估值，進(jìn)一步統(tǒng)計(jì)渭河流域年總降水和年平均氣溫未來(lái)時(shí)期各模型預(yù)估結(jié)果如圖8所示。

圖7 渭河流域日降水（（a）、（b））和日平均氣溫（（c）、（d））SDSM模型模擬結(jié)果（（a）、（c）率定期，（b）、（d）檢驗(yàn)期）

為減小GCM模式選取的差異性，采用每個(gè)典型濃度路徑下多模型平均值推求非一致性設(shè)計(jì)洪水［34］。RCP2.6、RCP4.5和RCP8.5三種典型濃度路徑下多模型平均序列及其趨勢(shì)線如圖9所示。對(duì)于年總降水，每個(gè)RCP下差別并不大，隨著輻射強(qiáng)迫的增強(qiáng)，序列逐漸由輕微下降轉(zhuǎn)變?yōu)檩p微上升，但幅度都較?。▓D9（a））；對(duì)于年平均氣溫，每個(gè)RCP下都呈現(xiàn)出一定的上升趨勢(shì)，隨著輻射強(qiáng)迫的增加，上升越顯著，RCP2.6情景下年平均氣溫首先有較為緩慢的上升，到2040年左右趨于穩(wěn)定；RCP4.5情景下2070年以前上升較為明顯，之后趨于穩(wěn)定；在最高輻射強(qiáng)迫的RCP8.5情景下，呈現(xiàn)出明顯的直線上升趨勢(shì)（圖9（b））。

4.3 非一致性設(shè)計(jì)洪水推求本文選取設(shè)計(jì)使用年限T為2015—2064年的50年。以時(shí)間為協(xié)變量時(shí)，設(shè)計(jì)使用年限t=2015，2016，…，2064年內(nèi)時(shí)變統(tǒng)計(jì)參數(shù)μt和σt可由最優(yōu)非一致性模型（表2）直接計(jì)算。給定水文風(fēng)險(xiǎn)R2015-2064，進(jìn)而迭代推求相應(yīng)的非一致性設(shè)計(jì)洪水z2015-2064。以Temp和Prep為協(xié)變量時(shí)，將設(shè)計(jì)使用年限t=2015，2016，…，2064年內(nèi)Temp和Prep的降尺度數(shù)據(jù)（多模型平均值）代入最優(yōu)非一致性模型（表2）計(jì)算時(shí)變統(tǒng)計(jì)參數(shù)μt和σt，進(jìn)而推求與R2015-2064相應(yīng)的非一致性設(shè)計(jì)洪水z2015-2064。同時(shí)給出了一致性假設(shè)成立情況下不同水文風(fēng)險(xiǎn)對(duì)應(yīng)的設(shè)計(jì)洪水。

結(jié)果表明，以氣象因子Temp和Prep為協(xié)變量時(shí)，3種典型濃度路徑下非一致性設(shè)計(jì)洪水差別很小，隨著輻射強(qiáng)迫的增加，設(shè)計(jì)洪水有輕微的增大，但并不顯著，因此下面僅以RCP2.6情景下設(shè)計(jì)結(jié)果為代表進(jìn)行分析。無(wú)論以時(shí)間t或者Temp和Prep為協(xié)變量，所得的非一致性設(shè)計(jì)洪水均明顯偏低于一致性假設(shè)情況。意味著在考慮非一致性的情況下，設(shè)計(jì)使用年限內(nèi)洪水極值事件相比于一致性假設(shè)將明顯減弱，從觀測(cè)期的洪水序列也發(fā)現(xiàn)明顯的下降趨勢(shì)（圖2（a）），說(shuō)明考慮非一致性的設(shè)計(jì)結(jié)果較為合理。

圖8 渭河流域RCP2.6、RCP4.5和RCP8.5三種典型濃度路徑下9個(gè)GCM模式［33］（左）年總降水和（右）年平均氣溫未來(lái)時(shí)期（2010—2099年）預(yù)估值

圖9 渭河流域RCP2.6、RCP4.5和RCP8.5三種典型濃度路徑下9個(gè)GCM模式年總降水和年平均氣溫未來(lái)時(shí)期（2010—2099年）多模型平均值

同時(shí)，兩種非一致性設(shè)計(jì)結(jié)果間也存在較大差別，以t為協(xié)變量的設(shè)計(jì)結(jié)果明顯偏低于以Temp和Prep為協(xié)變量的情況。雖然觀測(cè)期的洪水序列表明渭河流域的洪水極值事件確實(shí)存在一定的減弱趨勢(shì)，然而由于時(shí)間協(xié)變量情況假設(shè)所擬合出的下降趨勢(shì)將在未來(lái)時(shí)期延續(xù)下去，基于此外延假設(shè)得到的設(shè)計(jì)洪水不免有些夸張。

當(dāng)R2015-2064=10%，以t為協(xié)變量的非一致性設(shè)計(jì)洪水量級(jí)為3303 m3/s，意味著在未來(lái)2015—2064的50年間，超過(guò)3303 m3/s的洪水極值事件發(fā)生概率為10%，然而在觀測(cè)期1954—2009的56年間，此事件卻發(fā)生18次，說(shuō)明以t為協(xié)變量的設(shè)計(jì)結(jié)果減弱趨勢(shì)確實(shí)有所偏大，若以該設(shè)計(jì)結(jié)果為依據(jù)可能會(huì)增加工程事故風(fēng)險(xiǎn)。相比之下以Temp和Prep為協(xié)變量的設(shè)計(jì)結(jié)果6126 m3/s顯得更為合理。應(yīng)該指出，在假設(shè)一致性條件仍成立情況下，該設(shè)計(jì)值達(dá)到13 168 m3/s，而觀測(cè)期內(nèi)最大洪水量級(jí)也僅為5770 m3/s，更加印證了一致性設(shè)計(jì)結(jié)果的不合理性。

圖10 不同水文風(fēng)險(xiǎn)R2015-2064對(duì)應(yīng)的設(shè)計(jì)洪水z2015-2064及其95%置信區(qū)間

表4 設(shè)計(jì)洪水結(jié)果統(tǒng)計(jì) （單位：m3/s）

4.4 非一致性設(shè)計(jì)洪水不確定性采用Bootstrap方法推求設(shè)計(jì)洪水不確定性，抽樣結(jié)果見圖11，相應(yīng)的95%置信區(qū)間見圖10及表4。通過(guò)對(duì)華縣站進(jìn)行歷史大洪水調(diào)查，1933年歷史洪水量級(jí)8340 m3/s，從R2015-2064=1%設(shè)計(jì)洪水及相應(yīng)不確定性可見，以Temp和Prep為協(xié)變量的設(shè)計(jì)結(jié)果最為合理。

5 結(jié)論與討論

在水文風(fēng)險(xiǎn)概念下，將氣象因子引入到非一致性洪水頻率分析中，并在給定設(shè)計(jì)使用年限情況下推求對(duì)應(yīng)某一目標(biāo)水文風(fēng)險(xiǎn)的非一致性設(shè)計(jì)洪水及其不確定性。該方法與常見的以時(shí)間為協(xié)變量的方法一起應(yīng)用于有實(shí)測(cè)洪水資料流域的設(shè)計(jì)洪水研究，所得結(jié)論如下：（1）以氣象因子為協(xié)變量時(shí)，Gumbel分布為最優(yōu)洪水頻率分布，其他情況下均以P-Ⅲ分布為最優(yōu)洪水頻率分布，比較符合我國(guó)設(shè)計(jì)洪水研究的實(shí)際情況。（2）以時(shí)間為協(xié)變量時(shí)，最優(yōu)模型中位置參數(shù)為時(shí)間的線性函數(shù)；以氣象因子為協(xié)變量時(shí)，最優(yōu)模型中位置和尺度參數(shù)分別為年平均氣溫、年總降水以及年平均氣溫的線性函數(shù)。并且，以氣象因子為協(xié)變量的最優(yōu)非一致性模型優(yōu)于以時(shí)間為協(xié)變量的情況。（3）以時(shí)間為協(xié)變量和以氣象因子為協(xié)變量?jī)煞N情況下的非一致性設(shè)計(jì)洪水及其不確定性相比于一致性假設(shè)情況存在顯著差異。同時(shí)，兩種非一致性設(shè)計(jì)洪水及其不確定性之間也存在著明顯不同，前者由于相對(duì)不太合理的假設(shè)及參數(shù)擬合使得設(shè)計(jì)結(jié)果有所失準(zhǔn)，而后者則提供了更具物理意義、更為可信的設(shè)計(jì)結(jié)果。

圖11 設(shè)計(jì)洪水Bootstrap抽樣圖

本研究成果可為工程規(guī)劃設(shè)計(jì)及防洪決策提供一定參考依據(jù)，但其中也存在幾點(diǎn)值得進(jìn)一步深入研究之處：（1）不同頻率分布模型由于其尾部的差異性，對(duì)設(shè)計(jì)洪水量級(jí)的影響較大，因此，下一步有必要重點(diǎn)研究頻率分布模型的選取對(duì)水文分析計(jì)算的影響。（2）本文綜合比較了氣象因子年統(tǒng)計(jì)量與季節(jié)性統(tǒng)計(jì)量作為協(xié)變量的擬合效果，最終選取了效果更優(yōu)的年平均氣溫和年總降水作為氣象協(xié)變量，下一步有必要進(jìn)一步研究該變量與場(chǎng)次洪水之間的關(guān)系以及是否還有較之更為適合的協(xié)變量。（3）本文僅考慮了由于樣本數(shù)據(jù)有限所引起的設(shè)計(jì)洪水統(tǒng)計(jì)不確定性，對(duì)于分布模型選取以及引入GCM統(tǒng)計(jì)降尺度數(shù)據(jù)所增加的新的不確定性沒(méi)有進(jìn)行研究。（4）GCM輸出的未來(lái)時(shí)期氣象因子預(yù)估結(jié)果最長(zhǎng)可到2300年，時(shí)間跨度完全可以滿足現(xiàn)行規(guī)范規(guī)定的水利水電工程合理使用年限要求。然而，對(duì)于設(shè)計(jì)期較長(zhǎng)的情況下，GCM輸出的氣象協(xié)變量的預(yù)測(cè)仍存在精度較低等不足之處，因此下一步工作中考慮引入人類活動(dòng)因素作為協(xié)變量，結(jié)合流域未來(lái)水利工程規(guī)劃等影響，建立非一致性洪水頻率分布模型，增加設(shè)計(jì)結(jié)果的可靠性。本文所選研究區(qū)域面積較大，支流眾多，下一步考慮選取一些中小流域進(jìn)行比較研究，以探究本文方法的普適性。

［1］郭生練，劉章君，熊立華 .設(shè)計(jì)洪水計(jì)算方法研究進(jìn)展與評(píng)價(jià)［J］.水利學(xué)報(bào)，2016，47（3）：302-314.

［2］丁晶，何清燕，覃光華，等.論水庫(kù)工程之管運(yùn)洪水［J］.水科學(xué)進(jìn)展，2016，27（1）：107-109.

［3］宋松柏，李揚(yáng)，蔡明科.具有跳躍變異的非一致分布水文序列頻率計(jì)算方法［J］.水利學(xué)報(bào)，2012，43（6）：734-739.

［4］謝平，張波，陳海健，等.基于極值同頻率法的非一致性年徑流過(guò)程設(shè)計(jì)方法——以跳躍變異為例［J］.水利學(xué)報(bào)，2015，46（7）：828-835.

［5］徐宗學(xué)，張楠.黃河流域近50年降水變化趨勢(shì)分析［J］.地理研究，2006，25（1）：27-34.

［6］DU T，XIONG L H，XU C Y，et al.Return period and risk analysis of nonstationary low-flow series under cli?mate change［J］.Journal of Hydrology，2015，527：234-250.

［7］XIONG L H，JIANG C，DU T.Statistical attribution analysis of the nonstationarity of the annual runoff series of the Weihe River［J］.Water Science and Technology，2014，70（5）：939-946.

［8］VILLARINI G，SMITH J A，NAPOLITANO F.Non-stationary modeling of a long record of rainfall and tempera?ture over Rome［J］.Advances in Water Resources，2010，33（10）：1256-1267.

［9］STRUPCZEWSKI WG，KACZMAREK Z.Non-stationary approach to at-site flood frequency modelling II.Weighted least squares estimation［J］.Journal of Hydrology，2001，248（1/4）：143-151.

［10］STRUPCZEWSKI W G，SINGH V P，F(xiàn)ELUCH W.Non-stationary approach to at-site flood frequency modelling I.Maximum likelihood estimation［J］.Journal of Hydrology，2001，248（1/4）：123-142.

［11］STRUPCZEWSKI W G，SINGH V P，Mitosek H T.Non-stationary approach to at-site flood frequency model?ling.III.Flood analysis of Polish rivers［J］.Journal of Hydrology，2001，248（1/4）：152-167.

［12］KATZ R W，PARLANG M B，NAVEAU P.Statistics of extremes in hydrology［J］.Advances in Water Resourc?es，2002，25（8）：1287-1304.

［13］VILLARINI G，SMITH J A，SERINALDI F，et al.Flood frequency analysis for nonstationary annual peak re?cords in an urban drainage basin［J］.Advances in Water Resources，2009，32（8）：1255-1266.

［14］SALAS J D，OBEYSEKERA J.Revisiting the concepts of return period and risk for non-stationary hydrologic ex?treme events［J］.Journal of Hydrologic Engineering，2014，19：554-568.

［15］SCHUMANN A H.Flood risk assessment and management［M］.London：Springer，2011.

［16］梁忠民，胡義明，王軍，等.基于等可靠度法的變化環(huán)境下工程水文設(shè)計(jì)值估計(jì)方法［J］.水科學(xué)進(jìn)展，2017，28（3）：398-405.

［17］ROOTZéN H，KATZ R W.Design life level：quantifying risk in a changing climate［J］.Water Resources Re?search，2013，49：5964-5972.

［18］DAVISON A C，HINKLEY D V.Bootstrap methods and their application［M］.Cambridge University Press，UK，1997.

［19］左德鵬，徐宗學(xué)，李景玉，等.氣候變化情景下渭河流域潛在蒸散量時(shí)空變化特征［J］.水科學(xué)進(jìn)展，2011，22（4）：455-461.

［20］XIONG L H，DU T，XU C Y，et al.Non-stationary annual maximum flood frequency analysis using the norming constants method to consider non-stationarity in the annual daily flow series［J］.Water Resources Management，2015，29（10）：3615-3633.

［21］WILBY R L，DAWSON C W，BARROW E M.SDSM-a decision support tool for the assessment of regional cli?mate change impacts［J］.Environmental Modelling and Software，2002，17（2）：147-159.

［22］NASSERI M，TAVAKOL-DAVANI H，ZAHRAIE B.Performance assessment of different data mining methods in statistical downscaling of daily precipitation［J］.Journal of Hydrology，2013，492：1-14.

［23］趙芳芳，徐宗學(xué).統(tǒng)計(jì)降尺度方法和Delta方法建立黃河源區(qū)氣候情景的比較分析［J］.氣象學(xué)報(bào)，2007，65（4）：653-662.

［24］COLES S G.An introduction to statistical modeling of extreme values［M］.London：Springer，2001.

［25］GILROY K L，MCCUEN R H.A non-stationary flood frequency analysis method to adjust for future climate change and urbanization［J］.Journal of Hydrology，2012，414：40-48.

［26］OBEYSEKERA J，SALAS J.Quantifying the uncertainty of design floods under nonstationary conditions［J］.Journal of Hydrologic Engineering，2014，19（7）：1438-1446.

［27］RIGBY R A，STASINOPOULOS D M.Generalized additive models for location，scale and shape［J］.Journal of the Royal Statistical Society Series C-Applied Statistics，2005，54（3）：507-554.

［28］AKAIKE H.A new look at the statistical model identification［J］.IEEE Transactions on Automatic Control，1974，19（6）：716-723.

［29］BUUREN S V，F(xiàn)REDRIKS M.Worm plot：a simple diagnostic device for modeling growth reference curves［J］.Statistical in Medicine，2001，20：1259-1277.

［30］FILLIBEN J J.The probability plot correlation coefficient test for normality［J］.Technometrics，1975，17（1）：111-117.

［31］MASSEY F J Jr.The Kolmogorov-Smirnov test for goodness of fit［J］.Journal of the American Statistical Associa?tion，1951，46（253）：68-78.

［32］WILBY R L，DAWSON C W.SDSM 4.2-a decision support tool for the assessment of regional climate change im?pacts［Z］.User Manual，2007.

［33］van VUUREN D P，EDMONDS J，KAINUMA M，et al.The representative concentration pathways：an overview［J］.Climatic Change，2011，109（1/2）：5-31.

［34］CHEN L，F(xiàn)RAUENFELD O W.A comprehensive evaluation of precipitation simulations over China based on CMIP5 multimodel ensemble projections［J］.Journal of Geophysical Research：Atmospheres，2014，119.doi：10.1002/2013JD021190.