胡義明,楊晨曦,梁忠民,唐 超,陳 鈺,郭舉坤

(河海大學(xué)水文水資源學(xué)院,江蘇 南京 210098)

在分析工程對(duì)下游的防洪作用以及進(jìn)行水庫群聯(lián)合調(diào)洪計(jì)算時(shí),需要解決設(shè)計(jì)洪水的地區(qū)組成問題,以厘清當(dāng)設(shè)計(jì)斷面發(fā)生指定標(biāo)準(zhǔn)設(shè)計(jì)洪水時(shí),上游水庫及其區(qū)間的洪水空間組成情況。設(shè)計(jì)洪水地區(qū)組成的計(jì)算方法主要有地區(qū)組成法、頻率組合法和隨機(jī)模擬法[1]。在實(shí)踐應(yīng)用中,地區(qū)組成法中的同頻率組合法最為常用,即根據(jù)工程防洪要求,指定某一局部地區(qū)的洪量與下游控制斷面的洪量同為設(shè)計(jì)頻率,再根據(jù)水量平衡原則推求其余分區(qū)的洪量,如典型的“上游和下游同頻、區(qū)間相應(yīng)”或“下游和區(qū)間同頻、上游相應(yīng)”方法[2-3]。設(shè)計(jì)洪水地區(qū)組成分析是一種多維洪水特征變量聯(lián)合求解問題,諸多研究也采用Copula函數(shù)構(gòu)建不同區(qū)域洪水特征變量間的聯(lián)合分布函數(shù),以此推求設(shè)計(jì)洪水的地區(qū)組成[4-6]。應(yīng)用上述各類方法分析設(shè)計(jì)洪水地區(qū)組成時(shí),都需要洪水樣本系列滿足平穩(wěn)性要求。然而,由于氣候變化和人類活動(dòng)的影響,諸多站點(diǎn)的洪水樣本系列呈現(xiàn)非平穩(wěn)性變化,理論上,現(xiàn)行的基于平穩(wěn)性假定的工程水文計(jì)算理論與方法已不再適用于變化環(huán)境下的非平穩(wěn)性情形,給設(shè)計(jì)洪水地區(qū)組成分析帶來了挑戰(zhàn)[7-11]。目前,變化環(huán)境下水文設(shè)計(jì)值的計(jì)算主要有3種途徑:①基于水文極值系列重構(gòu)途徑;②基于分布函數(shù)加權(quán)綜合途徑;③基于變參數(shù)概率分布模型途徑[9]。諸多研究主要集中在單一洪水變量設(shè)計(jì)值計(jì)算或不同時(shí)段洪量組合設(shè)計(jì)值計(jì)算方面[12-16]。關(guān)于非平穩(wěn)性條件下,如何分析設(shè)計(jì)洪水的地區(qū)組成的研究較少,是工程水文領(lǐng)域亟須解決的難題。

本文基于Copula函數(shù)構(gòu)建不同區(qū)域非平穩(wěn)性洪量系列間的變參數(shù)多維聯(lián)合分布模型,結(jié)合條件期望組合法和等可靠度法,研究非平穩(wěn)性條件的設(shè)計(jì)洪水地區(qū)組成計(jì)算問題,并以寸灘站和宜昌站年最大15d洪量為例進(jìn)行了應(yīng)用分析。

1 理論與方法

1.1 非平穩(wěn)性邊緣分布模型

以我國設(shè)計(jì)洪水計(jì)算規(guī)范中推薦的皮爾遜三型分布(PE3)為基礎(chǔ),假定PE3分布中的某些參數(shù)隨其他因素(協(xié)變量)變化,構(gòu)建非平穩(wěn)性邊緣分布模型,以描述洪水特征的非平穩(wěn)性分布規(guī)律。非平穩(wěn)性PE3密度函數(shù)為[14]

(1)

其中αt=AZ+a0βt=BH+b0γt=CW+c0

式中:x為洪水特征變量;αt、βt、γt分別為PE3分布中t時(shí)刻的位置參數(shù)、尺度參數(shù)和形狀參數(shù);Γ(γt)為Gamma函數(shù);Z、H、W分別為影響位置參數(shù)、尺度參數(shù)和形狀參數(shù)變化的環(huán)境狀態(tài)變量因子集;A、B、C、a0、b0、c0為系數(shù)。

1.2 非平穩(wěn)性多維聯(lián)合分布模型

考慮到不同洪量自身的非平穩(wěn)性及其相依結(jié)構(gòu)的非平穩(wěn)性,基于Copula函數(shù)構(gòu)建可綜合分析上述兩類非平穩(wěn)性的變參數(shù)Copula模型[14],用于描述不同年份多變量聯(lián)合分布規(guī)律的演變特征。在變參數(shù)Copula模型中,各變量邊緣分布參數(shù)及Copula結(jié)構(gòu)參數(shù)隨協(xié)變量變化,不再是常數(shù)[14-17]。

本文以常用的Gumbel Copula函數(shù)為基礎(chǔ),構(gòu)建非平穩(wěn)性多維聯(lián)合分布模型,以分析非平穩(wěn)性條件下不同區(qū)域洪水特征變量聯(lián)合分布特征隨時(shí)間的演變規(guī)律。非平穩(wěn)性Gumbel Copula函數(shù)可表示為[14,18]:

C(μt,νt)=exp{-[(-lnμt)θct+(-lnvt)θct]1/θct} (θct≥1)

(2)

其中θct=exp(θD+θ0)+1

式中:C(·)為時(shí)變Gumbel Copula函數(shù);μt、vt為第t年變量xt和yt對(duì)應(yīng)的概率值;θct為第t年Gumbel Copula模型的結(jié)構(gòu)參數(shù);D為影響參數(shù)θct變化的協(xié)變量;θ和θ0為系數(shù)和截距。

1.3 條件期望組合法

基于構(gòu)建的不同區(qū)域間洪量(如上、下斷面洪量間或下斷面與區(qū)間洪量間)的非平穩(wěn)性Gumbel Copula函數(shù),可以推求給定下斷面設(shè)計(jì)洪量X=xp(xp為超過頻率p對(duì)應(yīng)的下斷面設(shè)計(jì)流量)條件下,上斷面洪量或區(qū)間洪量Y的條件概率密度函數(shù)?；谠摋l件概率密度函數(shù),可進(jìn)一步推求xp對(duì)應(yīng)的Y的條件期望值yt,E。(xp,yt,E)即為變量X和變量Y的條件期望組合[14]。

1.4 等可靠度法

等可靠度方法是由梁忠民等[19]提出,可用于推求非平穩(wěn)性條件下指定重現(xiàn)期的洪水設(shè)計(jì)值。以RS、RNS分別表示一致性、非一致性條件下,工程設(shè)計(jì)年限L和設(shè)計(jì)重現(xiàn)期T對(duì)應(yīng)的水文設(shè)計(jì)可靠度,其計(jì)算公式為

(3)

式中:Ft(XT,NS)為第t年水文極值的概率分布函數(shù);XT,NS為指定重現(xiàn)期T對(duì)應(yīng)的下斷面洪量設(shè)計(jì)值。

根據(jù)等可靠度方法:令RNS=RS,即可獲得非平穩(wěn)性條件下的XT,NS。通過變量X和變量Y的條件期望組合(xp,yt,E),可計(jì)算XT,NS條件下Y對(duì)應(yīng)的條件期望值yT,E、XT,NS和yT,E的差值為區(qū)間(或上斷面)的洪量;進(jìn)而獲得指定設(shè)計(jì)標(biāo)準(zhǔn)條件下,上、下斷面及區(qū)間的設(shè)計(jì)洪量組合。

1.5 模型評(píng)價(jià)

采用赤池信息量準(zhǔn)則(AIC)和貝葉斯信息準(zhǔn)則(BIC)對(duì)模型進(jìn)行評(píng)價(jià)。AIC和BIC是評(píng)估模型性能的常用指標(biāo)[20]。AIC建立在熵的概念基礎(chǔ)上,可以權(quán)衡所估計(jì)模型的復(fù)雜度和此模型擬合數(shù)據(jù)的優(yōu)良性 。模型參數(shù)個(gè)數(shù)越少,模型越簡(jiǎn)單,AIC值越小;對(duì)數(shù)似然值越大,模型精度越高,AIC值也越小。因此AIC準(zhǔn)則對(duì)模型的評(píng)價(jià)兼顧了模型的簡(jiǎn)潔性和精確性。AIC值越小,模型越好。

BIC準(zhǔn)則是對(duì)AIC準(zhǔn)則的改進(jìn),其將未知參數(shù)個(gè)數(shù)的懲罰權(quán)重由AIC準(zhǔn)則中的常數(shù)2變成了樣本容量。BIC值越小,模型越優(yōu)。

2 實(shí)例分析

以長江流域寸灘站和宜昌站1946—2014年的年最大15d洪量系列為例進(jìn)行實(shí)例分析。寸灘站位于長江干流上游,是三峽水庫入庫重要控制站,集水面積約86.66萬km2,宜昌站位于長江干流中游,集水面積約100.55萬km2(圖1)。

圖1 宜昌站和寸灘站位置信息

由圖2可知,寸灘站、宜昌站的洪量系列均存在下降趨勢(shì)。Mann-Kendall統(tǒng)計(jì)量的計(jì)算值分別為-2.6261和-3.319,其絕對(duì)值均大于0.05顯著性水平下對(duì)應(yīng)的閾值1.96,表明兩站年最大15d洪量系列減少趨勢(shì)顯著。

圖2 寸灘站、宜昌站年最大15d洪量系列

為了描述寸灘站和宜昌站年最大15d洪量的非平穩(wěn)性分布特征,采用了非平穩(wěn)性PE3分布模型。在模型構(gòu)建中,考慮到PE3分布函數(shù)中的參數(shù)必須大于0這一約束,參數(shù)均采用指數(shù)回歸形式。通過設(shè)定PE3分布函數(shù)中的位置參數(shù)和尺度參數(shù)隨時(shí)間變化,構(gòu)造了2種形式的非平穩(wěn)性PE3模型:分布函數(shù)中位置參數(shù)隨時(shí)間變化(PE3(L))、分布函數(shù)的位置參數(shù)和尺度參數(shù)隨時(shí)間變化(PE3(LS))(表1)。

表1 時(shí)變邊緣概率分布模型

關(guān)于非平穩(wěn)性邊際分布函數(shù)和變參數(shù)Copula模型中的參數(shù)估計(jì),均采用貝葉斯方法并結(jié)合馬爾科夫鏈蒙特卡洛方法(MCMC)抽樣技術(shù)進(jìn)行估計(jì),以最大后驗(yàn)估計(jì)作為模型使用參數(shù)[12,16-17]。在抽樣過程中,平行運(yùn)行5條鏈,每條鏈上抽取30000個(gè)樣本(每條鏈都能滿足收斂要求),去掉預(yù)熱的29000個(gè)樣本,每條鏈上僅采用最后的1000個(gè)樣本,5條鏈共計(jì)5000個(gè)樣本,其中使參數(shù)后驗(yàn)密度值達(dá)到最大的為參數(shù)的最大后驗(yàn)估計(jì)。

由表2可知,PE3(L)模型在寸灘站表現(xiàn)最優(yōu);而對(duì)于宜昌站,PE3(L)和PE3(LS)模型對(duì)應(yīng)的AIC和BIC值較為接近,但考慮到PE3(LS)模型具有更大的不確定性,最終選取PE3(L)模型描述宜昌站年最大15d洪量的非平穩(wěn)性分布特征[13]。表3給出了優(yōu)選的非平穩(wěn)性邊際分布函數(shù)和變參數(shù)Copula函數(shù)中的參數(shù)估計(jì)值。

表2 非平穩(wěn)性分布模型對(duì)應(yīng)的AIC和BIC值

基于非平穩(wěn)性Gumbel Copula函數(shù)模型分析了寸灘站和宜昌站年最大15d洪量系列的聯(lián)合分布隨時(shí)間的演變特征。圖3(a)給出了1990年、2000年、2010年、2020年和2030年5個(gè)典型年不同聯(lián)合概率(0.4、0.6、0.8、0.95)條件下,兩站年最大15d洪量聯(lián)合分布規(guī)律的演變特征。由圖3(a)可知,對(duì)于指定的聯(lián)合概率,5個(gè)典型年份的聯(lián)合分布是不同的,這就導(dǎo)致2個(gè)站點(diǎn)年最大15d洪量的最可能組合在不同年份不同。

圖3 寸灘站和宜昌站年最大15d洪量聯(lián)合分布特征

圖3(b)給出了1990年、2000年、2010年和2030年4個(gè)典型年不同聯(lián)合重現(xiàn)期(10a、20a、50a、100a)條件下,宜昌站和寸灘站年最大15d洪量聯(lián)合分布特征變化特征。由圖3(b)可知,對(duì)于給定的聯(lián)合重現(xiàn)期,2個(gè)站點(diǎn)年最大15d洪量的最可能組合在不同年份不同。

基于非平穩(wěn)性分布模型PE3(L),計(jì)算了0.01和0.02超過概率對(duì)應(yīng)的宜昌站年最大15d洪量的分位數(shù)(2015—2114年)。結(jié)合非平穩(wěn)性Copula模型,計(jì)算了不同超過概率條件下,與宜昌站年最大15d洪量對(duì)應(yīng)的寸灘站年最大15d洪量的條件期望值,進(jìn)而獲得了不同概率條件下宜昌站和寸灘站年最大15d洪量的組合樣本系列?；谠摻M合樣本系列,構(gòu)建了不同超過概率條件下兩站年最大15d洪量之間的條件期望組合關(guān)系曲線,如圖4所示。

圖4 不同超過概率條件下寸灘站和宜昌站年最大15d洪量的期望組合關(guān)系

根據(jù)等可靠度方法計(jì)算了不同重現(xiàn)期(50a和100a)條件下,不同工程設(shè)計(jì)壽命(10～100a)對(duì)應(yīng)的宜昌站年最大15d洪量的設(shè)計(jì)值。隨后基于不同超過概率下宜昌站和寸灘站年最大15d洪量的條件期望組合關(guān)系曲線,推求了不同重現(xiàn)期條件下,與宜昌站年最大15d洪量設(shè)計(jì)值對(duì)應(yīng)的寸灘站年最大15d洪量設(shè)計(jì)值。由于宜昌站的來水與寸灘站來水及區(qū)間來水要滿足水量平衡約束,為此宜昌站和寸灘站設(shè)計(jì)值之差即為區(qū)間來水設(shè)計(jì)值,結(jié)果如圖5所示。

圖5 不同重現(xiàn)期條件下宜昌站、寸灘站及區(qū)間的年最大15d洪量設(shè)計(jì)值

由圖5可知,在同一重現(xiàn)期條件下,宜昌站和寸灘站的洪量設(shè)計(jì)值隨著設(shè)計(jì)壽命的增加呈現(xiàn)減少趨勢(shì),這主要是由兩站的洪量系列本身呈現(xiàn)減少趨勢(shì)導(dǎo)致。如在宜昌站設(shè)計(jì)重現(xiàn)期為100a、工程設(shè)計(jì)壽命為10a和50a條件下,宜昌站、寸灘站及相應(yīng)區(qū)間洪量的設(shè)計(jì)值分別為(652.8,581.9,70.9)和(636.5,576.4,60.1)。在指定重現(xiàn)期條件下,宜昌站、寸灘站及區(qū)間的年最大15d洪量設(shè)計(jì)值隨著工程設(shè)計(jì)壽命增加而減小。

3 結(jié) 語

本文構(gòu)建了可綜合考慮不同區(qū)域洪量邊緣分布非平穩(wěn)性及區(qū)域洪量間相關(guān)結(jié)構(gòu)非平穩(wěn)性的時(shí)變動(dòng)態(tài)Copula模型,可用于分析不同區(qū)域洪水間的多維聯(lián)合分布規(guī)律隨時(shí)間的非平穩(wěn)性演變特征。綜合運(yùn)用等可靠度法和條件期望組合法,提出了變化環(huán)境下設(shè)計(jì)洪水地區(qū)組成分析方法,可用于變化環(huán)境下指定重現(xiàn)期對(duì)應(yīng)不同區(qū)域洪水組合設(shè)計(jì)值的推求。宜昌站和寸灘站的應(yīng)用結(jié)果表明,宜昌站、寸灘站及兩站區(qū)間洪水設(shè)計(jì)值均隨著工程使用年限的增加呈減小趨勢(shì),這與宜昌站、寸灘站的洪量系列本身的減少趨勢(shì)相吻合。

在示例分析中,關(guān)于非平穩(wěn)性Copula模型的構(gòu)建,限于資料條件,只考慮了邊緣分布參數(shù)和Copula結(jié)構(gòu)參數(shù)隨時(shí)間因子變化,而沒有考慮其他因素對(duì)模型參數(shù)的影響。在后續(xù)研究中,可考慮基于不同類型Copula函數(shù)、不同因素驅(qū)動(dòng)模型參數(shù)變化等方式,構(gòu)建多套不同的非平穩(wěn)性Copula模型,并通過優(yōu)選方式確定優(yōu)勢(shì)聯(lián)合分布模型。