吳紹飛，羅 文，熊凡迪，歐陽芬，黃彬彬

(南昌工程學院 鄱陽湖流域水工程安全與資源高效利用國家地方聯(lián)合工程實驗室，江西 南昌 330099)

1 研究背景

在氣候變化與人類活動的影響下，暴雨、洪水等水文極值事件常表現(xiàn)為顯著非一致性[1-3]，基于獨立同分布假設的一致性設計洪水計算結果將增加工程風險或導致經(jīng)濟上的浪費[4-5]，非一致性水文頻率計算方法成為水文科學研究的熱點[3]。

在眾多研究中，基于非一致性水文系列直接進行水文頻率計算越來越受到國內外學者的重視，此類方法主要包含兩類，即分布函數(shù)加權綜合法(如混合分布法、條件概率模型法等)和變參數(shù)模型法(如時變矩法)[6]。王軍等[7]重建了淮河流域530 a夏季降水系列，引入混合分布構建了一個包含4個子平穩(wěn)系列的綜合概率分布，結果表明，混合分布模型對觀測系列具有較好的擬合效果，在一定程度上解決了水文要素變異前后的非同分布問題，但分布函數(shù)參數(shù)較多，采用常規(guī)方法估計困難；宋松柏等[8]發(fā)展了適用于具有跳躍變異的非一致性概率分布頻率計算公式，此類方法適用于因數(shù)據(jù)缺失或氣候差異等造成的非一致性水文系列頻率計算問題；He等[9]以新疆馬斯納河流域為例，系統(tǒng)研究了P-Ⅲ分布、混合分布及條件概率模型法在非一致性洪水頻率分析中的適用性；Strupczewski等[10]探索了在非一致性洪水頻率分析中應用時變矩模型的可行性；Rigby等提出了基于位置、尺度和形狀的廣義可加GAMLSS(generalized additive models for location scale and shape)模型[11]，該法是時變矩法的進一步發(fā)展，在非一致性水文頻率計算中得到廣泛應用[11-13]。然而，在上述基于時變矩的非一致性水文頻率計算框架下，給定的設計標準或重現(xiàn)期所對應的設計值并不是一個常數(shù)，將隨著協(xié)變量的變化而變化，其結果很難直接應用到實際工程水文設計中[14]，重現(xiàn)期的期望超過次數(shù)法(expected number of exceedances,ENE)是解決上述問題的有效方法之一[15]。

本文以鄱陽湖修河虬津站1956-2014年最大1、3、7、15 d和1個月最大洪水系列(QMax 1d、QMax 3d、QMax 7d、QMax 15d和QMax 1m)為例，引入ENE方法，系統(tǒng)研究虬津站多時間尺度的設計洪水問題，以期為該區(qū)域洪災防治和水利工程管理提供依據(jù)。

2 數(shù)據(jù)來源與研究方法

2.1 研究區(qū)概況

修河(圖1)為鄱陽湖“五河”之一，地處亞熱帶季風氣候區(qū)，多年平均降水量約為1 660 mm，其中4-6月降水量約占全年降水量的50%。修河干流下游虬津水文站以上流域面積為9 914 km2，約占修河流域總面積的70%，該水文站實測最大洪峰流量為3 420 m3/s(1998年8月1日)。

圖1 鄱陽湖修河流域概況

2.2 數(shù)據(jù)來源

本研究主要采用以下3種數(shù)據(jù)：(1)實測水文數(shù)據(jù)；(2)研究區(qū)內實測降水量、氣溫等數(shù)據(jù)；(3)CMIP5(Coupled Model Intercomparison Project 5)大氣環(huán)流模型GCM(geheral circulation model)輸出數(shù)據(jù)，其來源分別如下。

(1)虬津站1956-2014年日平均流量數(shù)據(jù)來自江西省水文監(jiān)測中心，分別選取該時段年最大1 d、最大連續(xù)3 d、最大連續(xù)7 d，最大連續(xù)15 d和最大連續(xù)1個月的流量系列作為洪水極值進行研究(限于篇幅，僅列出最大連續(xù)1個月洪水總量系列Max 1m，見圖2(a)，下同)。

(2)選取影響徑流的兩個主要氣象要素(降水量和氣溫)為協(xié)變量，構建洪水系列非一致性分布模型。其中降水數(shù)據(jù)為利用泰森多邊形法計算的流域內有代表性的氣象站日降水系列面平均值，氣溫數(shù)據(jù)采用修河氣象站日氣溫數(shù)據(jù)，并進一步選取年均降水量(P)和年均氣溫(T)作為非一致性洪水頻率分析的協(xié)變量(圖2(b)、2(c))。

圖2 虬津水文站實測洪水系列及其與降水量、氣溫的相關關系

(3)GCM 輸出數(shù)據(jù)選取政府間氣候變化專門委員會(Intergovernmental Panel on Climate Change,IPCC)第5次評估報告中RCP4.5(representative concentration pathway 4.5)和RCP8.5(representative concentration pathway 8.5)兩種排放情景下40個模型2010-2099 年與美國國家環(huán)境預報中心(National Centers for Environmental Prediction,NCEP)相應的大尺度預報因子，具體處理過程見文獻[16]。

2.3 研究方法

2.3.1 基于位置、尺度和形狀參數(shù)的廣義可加模型 GAMLSS模型主要用于刻畫隨機變量統(tǒng)計參數(shù)與協(xié)變量之間的線性或非線性關系[10-11]，其主要原理如下：

假設GAMLSS模型在時刻t(t=1,2,…,n)的獨立隨機變量觀測值zt服從概率密度函數(shù)f(zt|θt)，其中，θt=(θt1,θt2,…,θtp)為t時刻的p個分布統(tǒng)計參數(shù)向量，θk=(θ1k,θ2k,…,θnk)T為所有時刻的第k個(k=1,2,…,p)統(tǒng)計參數(shù)所組成的向量，而gk(·)為一個單調可微的連接函數(shù)。則在忽略隨機效應項的前提下，θk為Yk的單調函數(shù)：

gk(θk)=ηk=Ykβk

(1)

式中：ηk為長度為η的向量;Ik為協(xié)變量個數(shù);Yk為n×Ik的矩陣;βk=(β1k,β2k,…,βIkk)T為長度為Ik的向量。

在實際進行非一致性水文頻率分析時，常用參數(shù)p≤3的分布函數(shù)，于是參數(shù)統(tǒng)計向量可以進一步表示為θt=(μt,σt,τt)，其中μt，σt，τt分別為GAMLSS模型中的位置參數(shù)、尺度參數(shù)和形狀參數(shù)。將各參數(shù)帶入函數(shù)式(1)，則gk(θk)可表達為：

(2)

其中，協(xié)變量矩陣Yk為：

(3)

將公式(3)代入公式(2)，得到如下函數(shù)關系(僅以μt為例)：

(4)

公式(1)～ (4)給出了以時間t為協(xié)變量的非一致性模型構建框架。由于以時間t為協(xié)變量的GAMLSS模型缺乏物理意義，本研究引入年均降水量(P)、年均氣溫(T)等氣象因子構建洪水系列協(xié)變量。同時，引入水文學中3種常用的分布函數(shù)：對數(shù)正態(tài)分布LOGNO、伽馬分布GA、耿貝爾分布GU，作為洪水極值系列備選理論分布(表1)，各模型擬合效果采用蟲圖和AIC值綜合判定。

表1 常用的兩參數(shù)分布函數(shù)表

(5)

由上式可知，K的期望值E(K)為1，因此E(K)的表達式為：

(6)

在一致性條件下，超過初始設計值的ENE重現(xiàn)期為T=1/p0。

在非一致性條件下，初始設計值Zp0的超過概率Pt會隨時間發(fā)生變化，因此期望E(K)的表達式為：

3 結果與分析

3.1 虬津站考慮氣象要素的非一致性GAMLSS模型

基于水文變異綜合診斷系統(tǒng)，周秋紅等[17]系統(tǒng)研究了虬津站1956-2014年多時間尺度徑流變異情況，結果表明，虬津站各時間尺度徑流系列均表現(xiàn)出顯著跳躍變異。在引起該站洪水系列變異的各類因素中，文獻[17]進一步指出氣象因素的要素主要是年總降水量，年均氣溫對研究區(qū)域短歷時洪水極值系列有顯著影響。因此，本文以流域內年均降水量P和年均氣溫T為協(xié)變量，以LOGNO，GA和 GU作為各列備選分布線型，并以 AIC值大小比較各模型擬合的優(yōu)劣，建立了虬津站上述各系列最優(yōu)非一致性GAMLSS統(tǒng)計模型，如表2所示。

以QMax 1d系列為例，由表2可知，考慮P和T為協(xié)變量的統(tǒng)計模型，AIC評價準則表明GA分布為最優(yōu)分布，且統(tǒng)計參數(shù)μ為P、T的線性函數(shù)(統(tǒng)計參數(shù)σ為常數(shù))時為最優(yōu)非一致性統(tǒng)計模型。進一步計算表明，在僅考慮線性函數(shù)關系的前提下，與以時間為協(xié)變量的模型相比，以P和T為協(xié)變量的GAMLSS模型具有更小的AIC值。研究區(qū)域QMax 3d和QMax 1m系列最優(yōu)分布函數(shù)也為GA分布，QMax 7d和QMax 15d系列最優(yōu)分布函數(shù)為LOGNO分布。杜濤等[18]指出，引入具有物理意義的氣象因子作為協(xié)變量，不僅能夠擬合出洪水系列的下降趨勢，還可以更好地描述洪水系列的周期波動。由表2還可以看出，各時間尺度洪水流量系列位置參數(shù)均存在顯著的非一致性，QMax 1d、QMax 3d和QMax 1m系列最優(yōu)分布模型位置參數(shù)均為年總降水量和年均氣溫的線性函數(shù)，QMax 7d和QMax 15d系列最優(yōu)分布模型位置參數(shù)均為年總降水量的線性函數(shù)。各時間尺度洪水流量系列尺度參數(shù)不存在顯著的非一致性(QMax 7d系列除外)。此外，為了與不考慮系列非一致性的結果進行對比，表2還給出了相應的一致性統(tǒng)計模型參數(shù)。

表2 虬津站實測各洪水系列一致性和非一致性頻率統(tǒng)計模型

3.2 修河流域降水量和氣溫統(tǒng)計降尺度成果

本研究降水量和氣溫統(tǒng)計降尺度數(shù)據(jù)來源于歐陽芬[16]的研究成果，對于修河流域歷史時期(1961-2005年)的年平均降水量和氣溫模擬預估值與實測值采用各自的均值和標準差來進行對比，結果如圖3、4所示。

從圖3、4可知，研究區(qū)域歷史時期(1961-2005年)多年月均降水量與多年月均氣溫的預測模擬值與實測值在大部分月份下均較為接近，但年內不同月份的數(shù)值波動較大，且降水量的擬合效果比氣溫相對更優(yōu)。

圖3 1961-2005年修河流域多年月平均降水量及其標準差的模擬值與實測值對比

圖5給出了兩種RCP典型濃度路徑(RCP4.5、RCP8.5)下輸出的40個GCM模式未來時期(2020-2099年)的年均降水量P和年均氣溫T的各模型預估結果，多模型平均系列如圖6所示。

圖4 1961-2005年修河流域多年月平均氣溫及其標準差的模擬值與實測值對比

通過對兩種典型濃度路徑下GCM模式輸出的40種氣候模式在(2020-2099年)未來時期降尺度的降水量和氣溫預估值比較分析，發(fā)現(xiàn)對于年均降水量，在兩種RCP下均呈現(xiàn)微小的增大趨勢(圖5(a)、5(c))，而對于年均氣溫，在兩種RCP下均呈現(xiàn)明顯的上升趨勢(圖5(b)、5(d))。對于多模型平均系列而言，年均降水量和年均氣溫均呈現(xiàn)明顯的上升趨勢，2020-2100年RCP8.5下年均氣溫顯著高于RCP4.5下年均氣溫(圖6)。

圖5 兩種RCP典型濃度路徑(RCP4.5、RCP8.5)下輸出的40種氣候模式的未來時期(2020-2099年)年平均降水量P和年平均氣溫T的各模型預估結果

圖6 兩種RCP典型濃度路徑(RCP4.5、RCP8.5)下未來時期(2020-2099年)年均降水量P和年均氣溫T多模型平均值

3.3 基于ENE的非一致性設計洪水研究

即可得出某一特定重現(xiàn)期下非一致性設計值(本文僅給出重現(xiàn)期為1～50 a的情形)，如圖7所示。

以重現(xiàn)期20和50 a為例，對圖7中兩種典型濃度路徑在考慮和不考慮非一致性條件下的設計洪水流量結果進行分析：(1)虬津站各系列不考慮非一致性的設計洪水流量均遠大于考慮非一致性的設計洪水流量。以Max 1d系列為例，當重現(xiàn)期為20 a時，一致性的設計洪水流量為7 394.8 m3/s，而RCP4.5和RCP8.5情景下非一致性設計洪水流量分別為4 981.3及4 505.9 m3/s；當重現(xiàn)期為50 a時，一致性的設計洪水流量為9 224.8 m3/s，相應的RCP4.5和RCP8.5情景下非一致性設計洪水流量分別為5 262.9及4 598.2 m3/s。虬津站其他時間尺度的設計洪水流量結果也有類似的規(guī)律，在選定的兩種典型排放情景下，考慮非一致性的設計洪水流量均小于不考慮非一致性的設計洪水流量。(2)RCP4.5和RCP8.5情景下，各時間尺度20和50 a重現(xiàn)期非一致性設計洪水流量與一致性設計洪水流量相比分別減小13.5%～32.6%(平均約26.4%)和15.3%～39.1%(平均約31.0%)。(3)各重現(xiàn)期下，虬津站在RCP4.5排放情景下的設計洪水流量均較RCP8.5情景相應值更大，這與杜濤[11]在清江和西江流域的研究結果類似。(4)由于修河干流虬津站上游柘林水庫對洪水的調節(jié)作用，虬津站考慮非一致性的設計洪水流量顯著小于不考慮非一致性的設計洪水流量，這與相關研究[11,18]關于未來非一致性設計洪水流量顯著高于不考慮非一致性條件下設計洪水流量的結果正好相反。

圖7 基于ENE法的虬津水文站各系列非一致性和一致性條件下不同重現(xiàn)期的設計洪水值對比

4 結 論

本文以鄱陽湖修河干流虬津水文站多時間尺度洪水系列為例，將影響徑流的氣象因子作為非一致性洪水概率分布函數(shù)參數(shù)的協(xié)變量，在重現(xiàn)期期望超過次數(shù)(ENE) 概念下，推求對應某一目標重現(xiàn)期的非一致性設計洪水流量，主要結論如下：

(1)GA分布為模擬實測洪水流量QMax 1d、QMax 3d和QMax 1m的最優(yōu)分布，LOGNO分布為模擬實測洪水流量QMax 7d和QMax 15d的最優(yōu)分布，各時間尺度洪水流量系列位置參數(shù)均存在顯著的非一致性。

(2)與一致性情形相比，以氣象因子為協(xié)變量的最優(yōu)非一致性統(tǒng)計模型擬合效果更佳，且兩者設計洪水流量值存在顯著差異。在選定的RCP4.5和RCP8.5兩種典型排放情景下，后者設計洪水結果均小于前者設計洪水結果。

(3)本文應用GCM的降尺度數(shù)據(jù)時仍存在較大的不確定性，其中氣溫更為明顯。同時，僅選取了RCP4.5和RCP8.5兩種中、高排放情景下的GCM數(shù)據(jù)，下一步將重點考慮其他不同典型排放情景以及非一致性設計洪水流量的不確定性。