應(yīng)用Copula函數(shù)的潘家口水庫設(shè)計(jì)洪水線推求

喀迪爾·麥麥提

（新疆英吉沙縣水利局，新疆 英吉沙 844500）

0 引言

洪水過程線隨機(jī)模擬是水利工程防洪設(shè)計(jì)的關(guān)鍵內(nèi)容，也是水利安全、洪水災(zāi)害控制的依據(jù)[1]。此前防洪設(shè)計(jì)中以傳統(tǒng)單變量倍比放大或同頻率設(shè)計(jì)方案為主，雖然其操作流程簡(jiǎn)易、各特征量重現(xiàn)周期易于滿足防洪標(biāo)準(zhǔn)，但仍需要指出的是存在明顯局限性：①單純考慮洪水特征量，忽視了期間非線性、相關(guān)性關(guān)系[2]；②將防洪安全標(biāo)準(zhǔn)與設(shè)計(jì)洪水過程線的特征量的頻率等同考慮，避開了洪水過程聯(lián)動(dòng)性等水文規(guī)律[3]。隨著計(jì)算科學(xué)的發(fā)展，國(guó)內(nèi)外學(xué)者提出了自回歸[4]、多元統(tǒng)計(jì)[5]、典型解集模型[6]、神經(jīng)網(wǎng)絡(luò)[7]等模型應(yīng)用于洪水過程隨機(jī)模擬，然而未能有效識(shí)別變量間共線性、依存關(guān)系，且難以實(shí)現(xiàn)多維變量聯(lián)合分布分析。Copula函數(shù)能夠聯(lián)立多變量間復(fù)雜關(guān)系，其邊緣分布靈活多樣，在水文研究中取得良好應(yīng)用。高超等[8]運(yùn)用用三維Copula函數(shù)以構(gòu)建洪水特征量聯(lián)合分布模型，并探究了洪水過程線形狀和三維對(duì)稱型Archimedean Copula函數(shù)公式對(duì)模擬精度的影響；劉和昌等[9]構(gòu)建了基于Copula函數(shù)的洪峰、洪量之間的函數(shù)關(guān)系，由此計(jì)算了峰量的邊緣分布頻率；李天元等[10]以隔延河水庫為例，對(duì)比分析了單變量同頻與多變量同頻洪水過程線設(shè)計(jì)的差異，結(jié)果表明，基于Copula函數(shù)構(gòu)造的聯(lián)合分布模型由于傳統(tǒng)同倍比方法。鑒于此，本文以二維Copula函數(shù)為基礎(chǔ)理論建立潘家口水庫峰量的聯(lián)合分布，并推求設(shè)計(jì)洪水，以期為水庫安全管理提供基礎(chǔ)信息。

1 原理與方法

1.1 Copula函數(shù)原理

Copula函數(shù)是基于多元分布理論而將聯(lián)合分布函數(shù)與邊緣分布函數(shù)連接構(gòu)建的一種相依性測(cè)度方法，對(duì)于傳統(tǒng)共線性問題不能準(zhǔn)確識(shí)別變量間相依信息，Copula能夠排除這些高度相關(guān)影響，以構(gòu)造二維離散分布模型進(jìn)行隨機(jī)模擬。Nelsen[11]將N元Copula函數(shù)定義域設(shè)定為[0～1]；其在洪水過程線模擬的應(yīng)用通常以1、2、3維度為主，本研究采樣雙變量同頻率設(shè)計(jì)方案，故針對(duì)二維Copula函數(shù)表達(dá)如下[12]：

式中：C 為 Copula 函數(shù)，U1(x1)，U2(x2)，…，Un(xn)在區(qū)間是連續(xù)分布的；則其相應(yīng)的聯(lián)合概率密度函數(shù)表示如下：

式中：u=FX(x),v=FY(x)分別代表分別為隨機(jī)變量X、Y的邊緣分布函數(shù)。洪水過程是洪峰、時(shí)段等特征量的函數(shù)，具有多元、非線性特征，因而二維變量聯(lián)合分布模擬多以Archimedean Copula函數(shù)函數(shù)簇應(yīng)用較多，且適應(yīng)性強(qiáng)，其能夠糾正共線性情形，辨識(shí)上尾相關(guān)性，其公式如下：

式中：u、v、r均為邊際分布函數(shù)，θ為Copula函數(shù)的參數(shù)。其密度函數(shù)為：

函數(shù)參數(shù)θ與Kendall秩相關(guān)系數(shù)τ關(guān)系表示如下：

1.2 邊緣分布頻率

洪峰、洪量是洪水過程的矢量之一，其隨時(shí)間變化符合PIII分布。已知洪峰Q、洪量W的分布函數(shù)為[13]：

式中：q0、w0分別表示洪峰、時(shí)段洪量給定值；a、β、a0依次為形狀、尺度、位置參數(shù)，其用于度量洪水特征。根據(jù)邊緣分布概率定理；則其二維離散型頻率計(jì)算公式如下：

式中：T為伽馬函數(shù)。

1.3 推求設(shè)計(jì)洪水過程線

以洪峰、洪量、洪水歷時(shí)等特征量構(gòu)建Copula聯(lián)合分布函數(shù)進(jìn)行聯(lián)合模擬，以描述洪水過程?；贑opula函數(shù)，在指定洪水歷時(shí)條件下，其分布如下[14]：

上式可轉(zhuǎn)化為：

假設(shè)Copula函數(shù)具有三變量重現(xiàn)期u、v、r，則在區(qū)間上存在無數(shù)個(gè)連續(xù)的u、v、r組合，然而實(shí)際洪水過程具有一定歷時(shí)，因此其組合個(gè)數(shù)的有限、可知的。遵循同頻倍比法，當(dāng)洪峰與時(shí)段洪量同頻時(shí)，即可滿足u=v=r條件，則洪水過程線聯(lián)合重現(xiàn)期可通過頻率組合C(u，v，r)來描述?；诖耍上冗\(yùn)用邊緣分布P-III函數(shù)推求洪峰與時(shí)段洪量的聯(lián)合設(shè)計(jì)值組合，在此基礎(chǔ)上進(jìn)行同倍比方法洪水過程，從而得到多遍了聯(lián)合分布模擬結(jié)果。

2 實(shí)例分析

2.1 潘家口水庫庫區(qū)概況

潘家口水庫位于灤河上游、唐山市燕山南麓，控制徑流面積8540 km2，是華北地區(qū)重要水利工程之一。庫區(qū)位于燕山迎風(fēng)坡降水中心地帶，年平均降水量介于400 mm～700 mm；區(qū)域水系發(fā)育充分、年內(nèi)分配不均，特別是7月～8月降水占全年降水的60%，為洪水多發(fā)季節(jié)。據(jù)多年暴雨資料測(cè)算顯示，區(qū)域洪水過程達(dá)3 d～4 d，洪峰歷時(shí)約30 h，落水在40 h～50 h之間，洪量占全年徑流量的20%～40%，洪水年際差異大。監(jiān)測(cè)期內(nèi)(1929年～2015年)最大洪峰流量出現(xiàn)在1962年，洪峰流量達(dá)18800 m3/s，最大六日洪量18.53億m3；1958年7月洪峰流量達(dá)9570 m3/s，最大六日洪量為16.52億m3；最近的一次洪水發(fā)生在2012年，洪峰流量達(dá)4280m3/s，相應(yīng)水位為226.69 m

2.2 潘家口水庫設(shè)計(jì)洪水參數(shù)選定

潘家口水庫兼具防洪、灌溉、飲水、養(yǎng)殖等多元功能，然而其設(shè)計(jì)初衷在于防患水患。其依據(jù)灤縣站、欒南站的1929年～1970年間42年序列監(jiān)測(cè)數(shù)據(jù)，并參考了1938年、1949年、1958年、1959年、1962年、1964年等年份的洪水資料。鑒于該項(xiàng)目的二期工程實(shí)施，本文將實(shí)測(cè)系列延伸至2015年；以3d、6d為基準(zhǔn)時(shí)段，運(yùn)用最大取樣法計(jì)算其設(shè)計(jì)流量，并以《潘家口水庫調(diào)度方式研究報(bào)告》為依據(jù)進(jìn)行復(fù)核，其設(shè)計(jì)結(jié)果見表1。計(jì)算了各項(xiàng)參數(shù)設(shè)計(jì)值與復(fù)核值之間的絕對(duì)誤差，結(jié)果表明其誤差介于0.8%～3.2%之間，說明設(shè)計(jì)洪水參數(shù)穩(wěn)定、可靠；另外各參數(shù)均值大部分略有偏小，只有極個(gè)別偏大，表明原設(shè)計(jì)成果相對(duì)保守、安全，基于多年來汛情特征和設(shè)計(jì)經(jīng)驗(yàn)，故仍以最初設(shè)計(jì)值為選定標(biāo)準(zhǔn)。

表1 潘家口水庫年最大設(shè)計(jì)洪水成果

2.3 基于Copula模型的洪峰流量聯(lián)合分布

洪水事件是典型的自然概率問題，最大似然估計(jì)可以精確求算一個(gè)樣本集的相關(guān)概率密度，因而被廣泛應(yīng)用于平穩(wěn)、非平穩(wěn)自然變量聯(lián)合分布研究中。鑒于此，參照文獻(xiàn)[11-12]中的方案，先應(yīng)用Gumbel-Hougaard Copula二維離散函數(shù)構(gòu)建3 d、6 d洪量的聯(lián)合分布模型，根據(jù)經(jīng)驗(yàn)頻率公式p=m(i)/(n+1)計(jì)算，并滿足n個(gè)樣本容量中有m(i)個(gè)聯(lián)合觀測(cè)樣本子集服從 x≤xi且 y≤yi；得到其 kendall秩相關(guān)系數(shù) τ分別為0.6321、0.5963；然后運(yùn)用非對(duì)稱 Gumbel-Hougaard Copula函數(shù)建立3 d、6 d洪量的三維聯(lián)合分布模型，基于最大似然法得到參數(shù)θ分別為0.2712、0.2518；雖然潘家口水庫上游建立了2個(gè)水庫工程，但其為微小型水庫、控制面積和庫容量均較小，因而可忽略其影響。潘家口水庫年最大洪峰與96 h洪量的邊緣分布(P-Ⅲ)頻率曲線見圖1，為便于直觀表達(dá)，統(tǒng)一采用升序排列；依圖可知其邊緣分布頻率符合指數(shù)模型 (yW=1760.9e-0.038x，R2=0.9966；yQ= 89.116e-0.038x，R2=0.9966)，擬合良好、精度可靠。

圖1 潘家口水庫年最大洪峰與96 h洪量頻率曲線

高超等[12]在類似于潘家口水庫群的設(shè)計(jì)洪水過程線研究中分別比較了 Clayton Copula、Gumbel Copula、Frank Copula 等三種函數(shù)對(duì)洪水特征量模擬效果的差異，多站數(shù)據(jù)均顯示以Frank Copula函數(shù)模擬精度最佳，基于此本文采用該函數(shù)計(jì)算潘家口水庫三變量間的理論聯(lián)合概率。圖2-a、2-b、2-c為依次為潘家口水庫的灤縣站、灤南站和二站綜合平均值的最大洪峰與96 h的洪量頻率分布曲線。從圖中可看出，各站和綜合評(píng)價(jià)資料的得出的洪水特征量的理論概率點(diǎn)和經(jīng)驗(yàn)概率值均密集分布于y=x直線附近，擬合函數(shù)為 (y=1.0043x+0.5142，R2=0.9966；y=1.0163x+0.3251，R2=0.998；y=1.0109x-0.9308，R2=0.9967)，均通過5%水平信度檢驗(yàn)，說明擬合結(jié)果較好，該模擬結(jié)果合理、可信。

圖2 潘家口年最大洪峰流量與168h 洪量的理論頻率與經(jīng)驗(yàn)頻率擬合圖

2.4 潘家口水庫協(xié)防調(diào)度規(guī)則

潘家口水庫目前執(zhí)行的是壩前各級(jí)控制水位與限制下泄流量相結(jié)合的調(diào)度方法，即當(dāng)洪水入庫時(shí)，以庫水位是否超過某一控制水位來判別洪水達(dá)到哪一級(jí)標(biāo)準(zhǔn)，在每一級(jí)控制水位之下則按相應(yīng)的防洪標(biāo)準(zhǔn)執(zhí)行某級(jí)限泄流量。潘家口水庫現(xiàn)狀洪水調(diào)度方式實(shí)施兩級(jí)控泄。具體調(diào)度原則為：

①上游入庫洪水達(dá)到50年一遇及以下標(biāo)準(zhǔn)洪水控泄10000 m3/s，保京山鐵路大橋安全。

②50年～500年一遇洪水控泄28000 m3/s，保潘家口電站安全。

③大于500年一遇洪水不控泄。

2.5 推求設(shè)計(jì)洪水過程線

在對(duì)潘家口水庫進(jìn)行補(bǔ)測(cè)時(shí)應(yīng)用階段用暴雨移置法、水汽入流指標(biāo)法及暴雨組合法，以1962年作為洪水典型年，設(shè)定聯(lián)合重現(xiàn)周期為1000a，另補(bǔ)充了50a的周期；于洪峰、洪量等值區(qū)間內(nèi)等距抽取200組設(shè)計(jì)值子集，經(jīng)通倍比放大后按照防洪調(diào)度規(guī)律演算其汛限水位，得到其限值依次為Zmax=229.6 m、Zmin=231.5 m，這一結(jié)果于1976年以來沿用至今?；诒?中給出的潘家口洪水特征量設(shè)計(jì)值(Q=40400 m3/s，W6h=36.7億m3)，經(jīng)同頻率放大推算了1955年設(shè)計(jì)洪水過程線，見圖3、表2。

可知，對(duì)1995年典型洪水過程推求，以單變量同頻率計(jì)算得到其最高水位為231.58 m，而雙變量同頻率得到其最高水位為231.69 m，這一結(jié)果表明兩種設(shè)計(jì)方法對(duì)汛限調(diào)度來講比較合理；但從防洪角度來看，單變量同頻設(shè)計(jì)相當(dāng)保守，不利于極端天氣條件下調(diào)洪協(xié)防。

表2 潘家口水庫千年一遇設(shè)計(jì)洪水調(diào)洪結(jié)果

鑒于1998年該流域發(fā)生特大洪水事件，另以該年為典型年，應(yīng)用雙變量同頻率并推求其洪水過程線，以驗(yàn)證該設(shè)計(jì)的合理性，結(jié)果見表2、圖4。由圖4可知1998年典型年潘家口水庫洪水過程歷時(shí)時(shí)間長(zhǎng)，最高洪峰流量高達(dá)21532 m3s-1，最高水位達(dá)230.304 m，洪峰歷時(shí)約30 h；入庫與出庫流量近似同步但略有滯后。由表2可知，1962年和1988年典型年中采用雙變量同頻設(shè)計(jì)的潘家口洪水特征參數(shù)值略高于單變量同頻設(shè)計(jì)方法，因此前者更利于庫區(qū)防洪；另外以雙變量的聯(lián)合重現(xiàn)期作為防洪標(biāo)準(zhǔn)，更符合洪水過程規(guī)律，因此該方案具有科學(xué)、合理性。

圖3 潘家口水庫千年、五十年一遇設(shè)計(jì)洪水過程線比較(1962典型年)

圖4 兩變量同頻率組合下潘家口水庫千年一遇設(shè)計(jì)洪水調(diào)洪過程線(1998年典型)

3 結(jié)論

區(qū)別于傳統(tǒng)單變量同頻設(shè)計(jì)洪水過程線獨(dú)立控制洪水各特征量重現(xiàn)周期，未能針對(duì)整個(gè)洪水過程聯(lián)合分布進(jìn)行模擬，本文利用Copula函數(shù)構(gòu)造了基于長(zhǎng)時(shí)歷史洪水資料的洪水特征量多變量聯(lián)合分布模型，該方案更符合防洪安全標(biāo)準(zhǔn)，適用于水文工程實(shí)際應(yīng)用。

洪水過程具有隨機(jī)非線性、多元化特征，Copula函數(shù)隨機(jī)模擬考慮了洪水特征量直接復(fù)雜相依關(guān)系?；陔p變量同頻率分析，經(jīng)水庫調(diào)洪演算后得到洪水歷時(shí)洪峰、洪量、最高水位，從而得到潘家口水庫洪水過程線。相較于單變量同頻模擬而言，以雙變量同頻推求的最高水位限值略大，并介于上下限制區(qū)間內(nèi)；由此可知，基于防洪安全角度，采用雙變量同頻組合設(shè)計(jì)能夠更好適應(yīng)于氣候變化背景下極端天氣引起的洪水不確定性。