梁小青，紀(jì)昌明，張驗(yàn)科，俞洪杰，閻曉冉

(1.華北電力大學(xué)可再生能源學(xué)院，北京102206；2.甘肅民族師范學(xué)院，甘肅合作747000)

0 引 言

洪水頻率分析及確定設(shè)計(jì)洪水是進(jìn)行水利工程防洪規(guī)劃及洪水管理的重要依據(jù)。隨著時(shí)間的推移，自然因素和人為因素共同作用，改變了天然洪水的分布規(guī)律，不能滿足傳統(tǒng)洪水頻率計(jì)算中對洪水一致性的要求。一般情況下，當(dāng)水文序列的變化規(guī)律在一定時(shí)期內(nèi)比較平穩(wěn)時(shí)，認(rèn)為水文序列是一致的；當(dāng)水文序列的變化規(guī)律在一定時(shí)期內(nèi)發(fā)生了變異(具有跳躍性或趨勢性的變化)，從一種平穩(wěn)狀態(tài)變化到另一種平穩(wěn)狀態(tài)時(shí)，認(rèn)為水文序列是非一致的。目前非一致水文序列頻率分析的途徑主要有還原(還現(xiàn))和直接計(jì)算兩種。還原是將變異后的序列修正到變異前狀態(tài)，還現(xiàn)是將變異前的序列修正到變異后的狀態(tài)，這兩種方法均存在一定缺陷[1- 2]。直接計(jì)算的途徑有時(shí)變矩法、混合分布法、條件概率法等。Strupczewski等[3- 4]最早提出時(shí)變矩法，該方法將序列的分布參數(shù)表示成時(shí)間t的某種函數(shù)，是一種從變化趨勢角度來分析序列非一致性的方法?；旌戏植挤ê蜅l件概率法主要從跳躍變異角度來分析序列的非一致性?；旌戏植挤ǖ幕舅枷胧鞘紫雀鶕?jù)序列變異時(shí)間點(diǎn)將整個(gè)研究序列劃分為若干滿足一致性的子序列，然后對各子序列的分布函數(shù)進(jìn)行加權(quán)綜合，得到一個(gè)混合分布函數(shù)，用以描述整個(gè)序列的分布特征。Singh等[5]第一次將混合分布用于非一致性洪水的頻率分析，取得了較好的效果。條件概率法[6]的原理是按照洪水物理成因的不同將年內(nèi)洪水劃分為互不重疊的若干時(shí)段，由各時(shí)段內(nèi)洪水最大值構(gòu)成一個(gè)樣本空間，根據(jù)全概率公式可得到發(fā)生某一量級洪水的概率。我國學(xué)者通過直接計(jì)算途徑也做了相關(guān)方面的研究。王軍等[7]應(yīng)用混合分布估計(jì)了淮河流域長系列降雨量的總體分布。劉丙軍等[8]利用時(shí)變矩法研究了西北江三角洲地區(qū)的洪水頻率。李新等[9]應(yīng)用條件概率法和混合分布法分別對王快水庫洪峰及三個(gè)時(shí)段洪量進(jìn)行了非一致性頻率分析。

在洪水頻率分析中，非一致性是需要考慮的一個(gè)方面；另一方面，洪水是一個(gè)復(fù)雜的水文事件，具有高度非線性、非平穩(wěn)性，兼有多個(gè)特征屬性且各屬性之間又有一定的相依性。而傳統(tǒng)頻率分析是針對單變量進(jìn)行的；這就使得分析結(jié)果不能夠全面、可靠地反映洪水事件，給水利工程規(guī)劃設(shè)計(jì)或運(yùn)行管理帶來了一定風(fēng)險(xiǎn)。因此，在研究洪水序列非一致性的同時(shí)，還需構(gòu)建以各特征屬性為變量的多變量聯(lián)合分布以較全面地評價(jià)洪水事件。傳統(tǒng)多變量聯(lián)合分布構(gòu)建的前提要求是邊緣分布必須為同一類型，而且對于聯(lián)合分布的類型一般為正態(tài)分布，較難求解非正態(tài)、多維變量的聯(lián)合分布。隨著1959年現(xiàn)代Copula理論的形成，利用Copula函數(shù)求解多變量聯(lián)合分布的方法被逐漸引入到水文領(lǐng)域。Copula函數(shù)不僅對用于構(gòu)建聯(lián)合分布的單變量的邊緣分布形式?jīng)]有要求，而且結(jié)構(gòu)清晰，能有效描述水文事件的內(nèi)在規(guī)律和特征屬性之間的相互關(guān)系[10]，在國內(nèi)外得到了廣泛應(yīng)用。Grimaldi S等[11]利用非對稱Archimedean copula函數(shù)建立了洪峰、洪量、歷時(shí)三變量的聯(lián)合分布。Papaioannou G等[12]采用Copula函數(shù)族求得多瑙河年最大洪峰與多個(gè)歷時(shí)洪量的聯(lián)合分布。熊立華等[13]用Clayton Copula求得了同一河流上兩站年最大洪水的聯(lián)合分布。李亦凡等[14]通過構(gòu)造Khoudraji copula求得了干旱歷時(shí)與干旱峰值的聯(lián)合分布。李天元等[15]以Gumbel Copula求得了宜昌站最大日流量、7天洪量和30天洪量的三維聯(lián)合分布。

現(xiàn)有利用Copula函數(shù)構(gòu)建洪水變量聯(lián)合分布的研究成果，一般為洪峰流量之間的聯(lián)合、歷時(shí)洪量之間的聯(lián)合、洪峰和歷時(shí)洪量之間的聯(lián)合等；同時(shí)考慮了變量非一致性的成果不多[16]，并且較少將洪水發(fā)生時(shí)間作為研究變量[17- 18]。雅礱江干流中下游河段是我國十三大水電基地之一。近些年來，下游河段上以錦屏一級為龍頭水電站的五級梯級電站的陸續(xù)建成，人為改變了天然洪水的時(shí)空分布過程，使其不再滿足一致性。目前，研究雅礱江流域洪水的文獻(xiàn)不多[19- 20]，且缺乏對該流域在變化環(huán)境下的洪水進(jìn)行非一致性分析的相關(guān)文獻(xiàn)。

針對上述問題，本文將洪水的非一致性與其發(fā)生時(shí)間聯(lián)系起來，以雅礱江下游錦屏一級水電站洪水為研究對象，采用混合分布法分析洪水的非一致性，并通過Copula函數(shù)求取洪水與其發(fā)生日期的聯(lián)合分布(分析流程如圖1)，以便能夠掌握錦屏一級水電站洪水與其發(fā)生時(shí)間的變化規(guī)律、二者的組合情況及一定量級洪水發(fā)生在不同時(shí)間或某一時(shí)間發(fā)生不同量級洪水的可能性。文中數(shù)據(jù)處理采用Matlab R2014a、EViews 8等軟件。

圖1 錦屏一級水電站非一致性洪水分析流程

1 非一致性洪水研究方法

1.1 洪水變量的混合分布

設(shè)年最大洪峰流量(或年最大時(shí)段洪量)為隨機(jī)變量X，研究序列年數(shù)為N，在序列發(fā)生變異的情況下，可用混合分布函數(shù)描述其變化特征。即

(1)

(2)

1.2 時(shí)間變量的Von Mises分布

設(shè)年最大洪峰流量(或年最大時(shí)段洪量)發(fā)生日期為隨機(jī)變量Y，汛期天數(shù)為n。Di(1≤i≤n)為從汛期開始第i天的日期，通過適當(dāng)?shù)娜亲儞Q，將日期轉(zhuǎn)換為方向數(shù)據(jù)y后用Von Mises分布函數(shù)F(y)表示為[21- 22]

(3)

式中，00；I0(k)為第一類變型Bessel函數(shù)零階值。

日期與方向數(shù)據(jù)的轉(zhuǎn)換及參數(shù)θ0的計(jì)算步驟如下：①在一單位圓上選定0°方向位置，從0°方向開始逆時(shí)針將日期按時(shí)間順序依次均勻分布在圓周上，則D1，D2，…，Dn對應(yīng)的圓心角弧度(即方向數(shù)據(jù))分別為2π/n，4π/n，…，2π；②確定第j(1≤j≤N)年的年最大洪峰流量(或年最大時(shí)段洪量)發(fā)生日期在單位圓周上對應(yīng)點(diǎn)的位置及角度θj；③通過式(4)求arctan(b/a)，并根據(jù)a和b的正負(fù)，得到角度θ0，對應(yīng)的日期為年最大洪峰流量(或年最大時(shí)段洪量)發(fā)生的平均日期。即

(4)

1.3 分布的檢驗(yàn)

以P-P圖法進(jìn)行邊緣分布的擬合效果評價(jià)。P-P圖法是將變量的經(jīng)驗(yàn)累積頻率與擬采用的分布函數(shù)得到的理論累積頻率點(diǎn)繪在一張圖中，觀察點(diǎn)據(jù)在45°方向的分布情況，并計(jì)算回歸系數(shù)R2。點(diǎn)據(jù)越靠近45°方向，R2越接近于1，表明擬采用的分布函數(shù)的擬合效果越好。

2 基于Copula函數(shù)理論的兩變量聯(lián)合分布

2.1 Copula函數(shù)優(yōu)選及參數(shù)估計(jì)

單參數(shù)對稱Archimedean Copula函數(shù)有確定的解析形式，在水文變量頻率分析中應(yīng)用廣泛，兩變量X和Y的聯(lián)合分布函數(shù)(X≤x和Y≤y同時(shí)發(fā)生的概率)表示為

F(x，y)=C(u1，u2)=
φ-1[φ(u1;θ)+φ(u2;θ)]

(5)

式中，C(u1，u2)為Copula求得的聯(lián)合分布函數(shù)，u1和u2分別為隨機(jī)變量X和Y的分布函數(shù)F(x)和F(y)；φ(·)為生成函數(shù)，φ-1(·)為φ(·)的逆函數(shù)；θ為參數(shù)。本文選取4種Archimedean Copula函數(shù)進(jìn)行優(yōu)選，包括Gumbel-Hougaard(GH)、Clayton、Frank和Ali-Mikhail-Haq(AMH)，并采用秩相關(guān)系數(shù)法[23]估計(jì)Copula函數(shù)參數(shù)θ。

2.2 Copula函數(shù)檢驗(yàn)

兩變量聯(lián)合分布函數(shù)的檢驗(yàn)也采用K-S檢驗(yàn)的方法，統(tǒng)計(jì)量如式(6)，檢驗(yàn)過程可參考1.3節(jié)。即

(6)

式中，mj是二維觀測值為(xj，yj)時(shí)，同時(shí)滿足x≤xj，y≤yj的二維觀測值的個(gè)數(shù)；Cj為觀測值(xj，yj)的Copula值。

通過均方根誤差(RMSE)法與AIC法[10]對聯(lián)合分布的擬合效果進(jìn)行評價(jià)。兩種方法的指標(biāo)值越小，說明分布的擬合效果越好。具體為

(7)

(8)

式中，d為Copula參數(shù)個(gè)數(shù)；Pc(·)為Copula聯(lián)合分布計(jì)算值；Pe(·)為聯(lián)合分布經(jīng)驗(yàn)值，Pe(j)=(mj-0.44)/(N+0.12)。

2.3 聯(lián)合超越概率分布及條件超越概率分布

在洪水事件中，通常關(guān)注的是指標(biāo)變量大于某一閾值的概率，為區(qū)別于概率論中變量概率分布函數(shù)的含義，本文對于兩變量X和Y，作如下定義：

聯(lián)合超越概率分布函數(shù)(X≥x和Y≥y同時(shí)發(fā)生的概率)為

F′(x，y)=1-u1-u2+C(u1，u2)

(9)

條件超越概率分布函數(shù)(給定X≥x時(shí)，Y≥y發(fā)生的概率或給定Y≥y時(shí)，X≥x發(fā)生的概率)為

F′Y|X(y|x)=[1-u1-u2+
C(u1，u2)]/(1-u1)

(10)

F′X|Y(x|y)=[1-u1-
u2+C(u1，u2)]/(1-u2)

(11)

根據(jù)2.2節(jié)優(yōu)選出的Copula函數(shù)可求得上述各函數(shù)值。

3 實(shí)例分析

3.1 研究區(qū)概況及數(shù)據(jù)資料

雅礱江流域水能資源極其豐富，干流下游規(guī)劃有錦屏一級、錦屏二級、官地、二灘和桐子林5座水電站。其中，由錦屏一級、錦屏二級、官地構(gòu)成的“錦官電源組”為國家電力調(diào)度控制中心直調(diào)電站，錦屏一級為控制性電站，是一座以發(fā)電為主，并兼有分擔(dān)長江中下游地區(qū)防洪任務(wù)的年調(diào)節(jié)電站。

錦屏一級水電站水庫正常蓄水位為1 880 m，正常蓄水位下水庫庫容為77.6億m3，屬于大(1)型水利樞紐工程。考慮到大型水庫的調(diào)節(jié)性能高，其防洪調(diào)度主要通過洪量進(jìn)行控制。因此，本文依據(jù)錦屏一級水電站1972年～2010年入庫洪水的年最大時(shí)段洪量來分析洪水的非一致性。經(jīng)數(shù)據(jù)資料審查，滿足可靠性和代表性要求。

3.2 年最大時(shí)段洪量非一致性分析

3.2.1變異點(diǎn)的確定及顯著性檢驗(yàn)

采用累積距平法[24]對錦屏一級水電站年最大1日洪量序列和年最大3日洪量序列的變異情況進(jìn)行定性分析(見圖2)。

爸跟我說：“你媽走了，我怎么辦？”爸的頭發(fā)一下子全白了。我感覺他不是以前那個(gè)老爸了。終于，有一天，他也病了，也是腎結(jié)石。我們一家都有腎結(jié)石，可能是水的原因。老秦給他配了草藥，說老人年紀(jì)大了，最好保守治療，用草藥給化了?？墒遣菟帥]有見效，爸疼得喘不上氣。老秦說，不能再拖了，趕快住院做手術(shù)。

圖2 年最大時(shí)段洪量累積距平過程

根據(jù)宋松柏等人的建議[25]，考慮將序列分成兩段；因此，在圖2中年最大1日洪量序列、年最大3日洪量序列的累積距平過程曲線上，均取最小極值點(diǎn)對應(yīng)的年份1989年作為可能變異點(diǎn)來進(jìn)行研究。

采用滑動(dòng)t檢驗(yàn)法[26]對1989年變異顯著性進(jìn)行定量分析。在1989年前后分別獲取步長為10a的子序列，檢驗(yàn)結(jié)果見表1。

表1 年最大時(shí)段洪量滑動(dòng)t法檢驗(yàn)結(jié)果

注：顯著性水平α為0.1時(shí)，臨界值為1.734 1。

由表1可知，年最大1日洪量序列和年最大3日洪量序列在1989年均發(fā)生顯著變異，兩個(gè)時(shí)段洪量序列均存在非一致性情況。下面以年最大1日洪量為例，說明非一致洪水頻率與當(dāng)前設(shè)計(jì)頻率的差異。

3.2.2子序列及平穩(wěn)性檢驗(yàn)

將年最大1日洪量序列以變異年份1989年劃分為1972年～1989年和1990年～2010年兩個(gè)子序列。檢驗(yàn)序列是否平穩(wěn)的方法有時(shí)序圖法、樣本自相關(guān)函數(shù)法及單位根法等。其中，單位根法中的ADF檢驗(yàn)[27]是目前應(yīng)用較多且最為有效的。本文采用ADF檢驗(yàn)法對子序列進(jìn)行平穩(wěn)性檢驗(yàn)(見表2)，顯著性水平取0.05。

表2 子序列ADF檢驗(yàn)結(jié)果

由表2可知，通過ADF檢驗(yàn)得到兩個(gè)子序列不平穩(wěn)的概率均小于顯著性水平，可認(rèn)為兩個(gè)子序列是平穩(wěn)和一致的。

3.3 非一致年最大1日洪量及其發(fā)生日期的分布

將年最大1日洪量作為洪水變量X，設(shè)1972年～1989年和1990年～2010年兩子序列的分布函數(shù)分別為F1(x)和F2(x)。選取一年中的6月1日～9月30日為汛期，年最大1日洪量發(fā)生在汛期的日期為時(shí)間變量Y。兩子序列分布參數(shù)、當(dāng)前設(shè)計(jì)分布參數(shù)及年最大1日洪量發(fā)生日期分布的參數(shù)計(jì)算結(jié)果見表3。取顯著性水平α為0.05，年最大1日洪量與其發(fā)生日期的K-S檢驗(yàn)結(jié)果見表4。

將兩變量各自分布的經(jīng)驗(yàn)頻率和由式(1)～式(3)得到的對應(yīng)理論頻率點(diǎn)繪于圖3中。由圖3a和圖3b可知，各點(diǎn)在45°方向分布均勻，回歸系數(shù)R2均在0.95以上，理論分布均擬合較好。

另外，本文將考慮變異的混合分布曲線、當(dāng)前設(shè)計(jì)分布曲線、經(jīng)驗(yàn)頻率點(diǎn)據(jù)置于同一圖中(見圖4)。由圖4中曲線擬合情況可知，混合分布曲線與經(jīng)驗(yàn)頻率點(diǎn)據(jù)在整體上擬合較好，特別是頻率較小的中上部，而當(dāng)前設(shè)計(jì)分布曲線只對頻率大于20%的點(diǎn)據(jù)擬合較好，說明混合分布模型較當(dāng)前設(shè)計(jì)分布合理。

3.4 Copula函數(shù)優(yōu)選及兩變量聯(lián)合分布

經(jīng)SPSS軟件分析得到年最大1日洪量與發(fā)生日期兩變量的Kendall秩相關(guān)系數(shù)τ為0.118，可知二者為正相關(guān)關(guān)系，各Copula函數(shù)參數(shù)、統(tǒng)計(jì)量及擬合指標(biāo)值計(jì)算結(jié)果見表5。

由表5可知：4種Copula函數(shù)的K-S檢驗(yàn)統(tǒng)計(jì)量Z值都小于臨界值0.218，均通過檢驗(yàn)；但Clayton Copula的統(tǒng)計(jì)量值最小，檢驗(yàn)結(jié)果最為顯著。從擬合指標(biāo)值來看，RMSE值和AIC值最小的也是Clayton Copula。Clayton Copula也適合構(gòu)造具有正相關(guān)關(guān)系的變量之間的聯(lián)合分布[17]。

圖3 非一致年最大1日洪量及發(fā)生日期示意

圖4 不同分布擬合效果

Copula參數(shù)θK-S統(tǒng)計(jì)量Z擬合指標(biāo)RMSEAICGH1.1340.1130.043-242.554Clayton0.2680.1120.040-248.860Frank1.0730.1150.041-245.782AMH0.6480.1260.044-240.907

綜合上述考慮，本文以Clayton Copula函數(shù)來構(gòu)建非一致年最大1日洪量與其發(fā)生日期的聯(lián)合分布，并求得聯(lián)合超越概率分布，其三維曲面圖和典型概率等值線圖如圖5、6所示。

由圖5，一方面可得到超過某一量級洪水發(fā)生在不同日期后的概率及發(fā)生在某時(shí)間段的概率。如，W≥2.5億m3的洪水在6月4日后發(fā)生的概率為96.4%，在7月1日后發(fā)生的概率為83.2%；則發(fā)生在6月4日～7月1日之間的概率就是13.2%。另一方面，也可得到某一日期后發(fā)生超過一定量級洪水的概率及發(fā)生某一量級段洪水的概率，如，8月10日后發(fā)生W≥2.8億m3洪水的概率為42.8%，發(fā)生W≥4.25億m3洪水的概率為25.8%，則發(fā)生2.8～4.25億m3洪水的概率為17%。

圖5 聯(lián)合超越概率分布

圖6 聯(lián)合超越概率分布等值線

由圖6中等值線可得到給定某一概率時(shí)，洪水與其發(fā)生日期的不同組合情況。另外，等值線的疏密程度說明了洪水發(fā)生的可能性大小，等值線越密，發(fā)生的可能性越大。由圖可知，縱坐標(biāo)軸上7月20日～8月20日的等值線較密集，對應(yīng)橫坐標(biāo)軸上洪水的等值線分布也較密，說明該電站較容易在7月下旬到8月下旬這個(gè)時(shí)間段內(nèi)發(fā)生3.5億～5億m3大小的洪水。

3.5 條件超越概率密度

根據(jù)式(10)～(11)可得到Clayton copula函數(shù)的條件超越概率密度

(12)

(13)

圖7a、7b中各繪制了3條非一致年最大1日洪量和發(fā)生日期的條件超越概率密度曲線，其中圖7a為給定洪水量級X的情況下，發(fā)生在日期Y后的條件概率密度曲線，圖7b為給定日期Y的情況下，發(fā)生超過X量級洪水的條件概率密度曲線。

圖7 條件超越概率密度

圖7a中不同量級洪水發(fā)生的概率密度均隨著日期的向后推移呈先增大后減小的變化趨勢；曲線最高點(diǎn)處的概率密度反映了年最大1日洪量發(fā)生在對應(yīng)橫坐標(biāo)日期附近的概率也大。如，4億m3量級洪水的條件超越概率密度最大值為0.251，說明該量級洪水發(fā)生在8月9日左右的可能性最大。根據(jù)各曲線較大概率密度所在范圍，得出錦屏一級水電站的主汛期在7月19日～8月30日左右，這與實(shí)際情況相一致。前汛期小量級洪水發(fā)生的概率大于大量級洪水發(fā)生的概率，主汛期和后汛期的情況與前汛期相反。同理，圖7b中不同日期發(fā)生洪水的概率密度均隨著洪水量級的增加呈先增大后減小的變化趨勢；曲線最高點(diǎn)的概率密度反映了在該日期發(fā)生對應(yīng)橫坐標(biāo)量級左右洪水的概率最大。如，9月1日最有可能發(fā)生4.5億m3量級左右的洪水。

4 結(jié) 論

本文以雅礱江下游錦屏一級水電站的洪水為研究對象，分析了其在變化環(huán)境下的非一致性，得出以下主要結(jié)論：

(1)利用累積距平法與滑動(dòng)t檢驗(yàn)法對年最大1日洪量和年最大3日洪量時(shí)間序列的變異點(diǎn)的定性分析和定量判斷結(jié)果表明，兩個(gè)序列均在1989年發(fā)生顯著變異。

(2)年最大1日洪量的P-Ⅲ混合分布與當(dāng)前設(shè)計(jì)分布相比，更符合經(jīng)驗(yàn)頻率點(diǎn)據(jù)的分布，隨著頻率的減小，同頻率的P-Ⅲ混合分布洪量值要小于當(dāng)前設(shè)計(jì)分布的洪量值。

(3)通過K-S檢驗(yàn)與RMSE、AIC擬合評價(jià)指標(biāo)對4種Archimedean Copula進(jìn)行優(yōu)選，得出Clayton Copula函數(shù)的檢驗(yàn)結(jié)果最為顯著，擬合效果最好，適合用來求解非一致年最大1日洪量與其發(fā)生日期兩變量的聯(lián)合分布。

(4)由Clayton Copula得到了非一致年最大1日洪量與其發(fā)生日期的聯(lián)合分布，并求得了聯(lián)合超越概率分布、條件超越概率密度等，從不同角度挖掘了二者的變化規(guī)律。

本文研究可為錦屏一級水電站設(shè)計(jì)洪水的修正、汛期洪水分期、防洪調(diào)度等提供一定的參考依據(jù)。