程 旭，馬細霞，肖 遙，2，王倩麗

（1.鄭州大學水利科學與工程學院，鄭州450001；2.中水東北勘測設計研究有限責任公司，長春130021）

0 引言

設計洪水是水利工程規(guī)劃設計、防洪策略和水資源管理的重要水文數(shù)據，由流量資料推求設計洪水是設計洪水計算的主要方法［1，2］。傳統(tǒng)上，通過擬合各歷史洪水要素觀測值概率分布函數(shù)以估計其不同重現(xiàn)期下的單變量設計洪水值，以其峰值、洪量和持續(xù)時間這3個內部存在物理聯(lián)系而并不相互獨立的隨機變量分析洪水過程，對于一個完整的洪水事件，僅依靠單變量概率分布分析存在一定的局限性，不能完整描述洪水過程［3，4］。因此要進行更完整的設計洪水頻率分析，就需要對各洪水要素進行多變量聯(lián)合統(tǒng)計分析［5，6］。設計洪水地區(qū)組成是流域開發(fā)方案設計以及上下游、干支流工程聯(lián)合調洪運行的重要依據［7，8］，目前，常用的地區(qū)組成方法有同頻率法和典型年法，該類方法主觀地考慮“對防洪不利”［9-11］。雖然便于使用和規(guī)劃，但更多是主觀臆斷而缺少客觀的判斷，不能充分反映各分區(qū)之間的空間相關性，也無法定量地描述各分區(qū)與設計斷面之間的相關程度［12］。Copula 函數(shù)是將變量聯(lián)合累積分布函數(shù)同變量邊緣累積分布函數(shù)連接起來的函數(shù)，能夠構建各洪水要素之間的聯(lián)合分布，也能夠構建各分區(qū)洪水之間的聯(lián)合分布［13-15］。

為求得同時考慮洪水要素相關性和空間相關性的設計洪水結果，本文以黃河中游河南段洪水地區(qū)組成復雜的花園口水文站上游為研究對象，選取12天洪量、洪峰流量兩個洪水要素，基于Copula 函數(shù)建立兩個洪水要素以及不同斷面之間洪水要素的聯(lián)合分布，并與遺傳算法相結合構建Copula-GA（Copula-Genetic Algorithms）算法，推求最有可能的設計洪水地區(qū)組成，以期為防洪預案制定提供參考依據。

1 研究區(qū)概況及數(shù)據來源

花園口水文站位于黃河中游，控制流域面積73 萬km2，占黃河流域總面積的92%［16］?；▓@口斷面上游干流河南段建有三門峽和小浪底兩座大型水庫，流域水系圖如圖1所示，研究區(qū)水系縱橫，呈不對稱分布，水患災害頻發(fā)，各量級洪水以上大洪水和上下同時大洪水為主占95.2%，即以全流域性洪水為主［17］。各防洪斷面洪水組成成分的確定對于水庫間防洪優(yōu)化聯(lián)調工作的開展起著至關重要的作用［18］?；▓@口斷面的洪水概率分布直接受到干支流水庫調節(jié)作用的影響［19］。所以，應與黃河中游各水系實際情況相結合，確定花園口洪水地區(qū)組成方案。本文將以花園口站過境洪量為控制標準，探求花園口以上，小浪底和相應區(qū)間洪水的相關關系及其地區(qū)組成。

圖1 研究區(qū)水系圖Fig.1 Water system diagram of the study area

本文基本資料為1950-2015年黃河干流河南段的花園口、小浪底以及小-花區(qū)間各站點的實測流量。1958年7月份洪水是花園口水文站記載以來的最大洪水，花園口站洪峰流量高達22 300 m3/s，該場洪水預見期短，對下游威脅較大，具有一定的代表性，對該場實測洪水過程進行插值，對洪峰所在時段取實測洪峰值，得到分辨率為6 h的洪水典型洪水過程。

2 研究方法

本文基于Copula 函數(shù)特性以及研究區(qū)洪水情況，建立花園口站洪量和洪峰流量的聯(lián)合分布，計算相應重現(xiàn)期的時段洪量和洪峰流量，并推求設計洪水過程線；采用Copula 函數(shù)構造各干支流洪水要素之間的聯(lián)合分布，并與遺傳算法結合構建Copula-GA算法，推算各干支流設計洪水的最可能地區(qū)組成。

2.1 兩變量重現(xiàn)期

選取時段洪量、洪峰流量兩個變量作為洪水特征要素，應用一組隨機變量（Qi，Wi）描述相互獨立的洪水事件系列｛I1，I2，...，In｝中的每一場洪水事件Ii，假設其聯(lián)合分布為F（q，w）。在設計洪水要素推求過程中，洪峰和洪量的表示方式為：

將以上單變量事件整合為聯(lián)合變量事件，表達方式為：

式中：Q和W分別表示洪水事件的洪峰流量和洪量；洪水事件Iq，w表示Q，W至少一個超過設定閾值。

若Q和W為年最大值，其單變量重現(xiàn)期為Tq、Tw［20］：

此時，Eq，w事件的重現(xiàn)期為Tq，w：

FQ（q）和FW（w）均大于或等于F（q，w），由（4）和（5）可得：

假定采用單變量方法得到的T年一遇設計洪峰qT和洪量wT，即變量Q和W分別取T年一遇的設計值，根據式（6）可得：

從上式可知，若要求Tq，w（qT，wT），則必定q和w分別大于或等于qT和wT。

2.2 建立洪峰流量和洪量的聯(lián)合分布

選取洪峰流量以及洪量兩個洪水要素，使用Copula 函數(shù)建立兩者的聯(lián)合分布，本文以GH-Copula函數(shù)為例［21，22］：

式中：FQ（q）和FW（w）分別為洪峰流量Q和洪量W的邊緣分布；θ為函數(shù)參數(shù)，變量之間的相關性越強，其值越大，θ的參數(shù)估計為：

式中：τ為Kendall 秩相關系數(shù)；（xi-xj）為觀測點據；sign（*）為符號函數(shù)，當（xi-xj）（yi-yj）＞0 時，sign（*）=1；當（xi-xj）（yi-yj）＜0 時，sign（*）=-1；當（xi-xj）（yi-yj）=0時，sign（*）=0。

對于特定重現(xiàn)期Tq，w，將式（8）帶入式（5）得［23］：

式中：Tu，v=Tq，w?？梢詽M足式（13）的q，w的無數(shù)種組合，在平面上繪制q，w關系圖，會出現(xiàn)相應重現(xiàn)期下的等值線，借助于式（8）～（10），式（13）可以表示為：

式中：v=L（u）為相應重現(xiàn)期下等值線的表達式。

根據式（14），可以得到重現(xiàn)期為Tq，w洪峰流量q與洪量w的等值線。該等值線上所有點的重現(xiàn)期均為Tq，w，當Tq，w=T時，該等值線上點（q，w）分別大于相應的qT和wT值。在設計洪水的洪峰流量頻率u與洪量頻率v相等時，有唯一解，即等值線v=L（u）與直線v=u的唯一交點，其解析表達式為：

式中：qTO和wTO為兩變量重現(xiàn)期Tq，w（Tq，w=Tu，v）對應的Q和W的設計值。

2.3 推求設計洪水過程線

根據所得到的設計洪峰流量qTO、設計洪量wTO以及典型洪水過程，采用變倍比放大法獲得設計洪水過程［24］：

式中：DF（t）、TF（t）分別為設計洪水和典型洪水在t時刻的流量；QTD和WTD分別為典型洪水過程的洪峰流量和歷時為TD的洪量。

2.4 設計洪水地區(qū)組成及遺傳算法的優(yōu)化

2.4.1 洪水地區(qū)組成的數(shù)學描述

設計洪水地區(qū)組成的核心內容是在匯流洪量一定的情況下，確定洪水來源的各分區(qū)洪水所占比例［25，26］。根據流域水系分布情況，本文重點研究花園口水文站、小浪底水庫及小-花區(qū)間的設計洪水地區(qū)組成。如圖2所示，A為小浪底水庫，X為小浪底入庫洪水；B為小浪底-花園口區(qū)間，Y小浪底-花園口區(qū)間洪水；Z為花園口水文站斷面，其中Z=X+Y。

圖2 研究區(qū)內洪水地區(qū)組成示意圖Fig.2 Schematic diagram of the flood region composition of study areas

本文各分區(qū)獨立選樣并采用P-III型分布作為頻率曲線，其概率密度函數(shù)為［27］：

式中：α、β、a0分別為P-III 分布的形狀、尺度以及位置參數(shù)，α＞0，β＞0；Γ(α)為α的伽馬函數(shù)。

α、β、a0與總體的統(tǒng)計參數(shù)Xˉ、Cv、Cs的關系為：

2.4.2 同頻率地區(qū)組成法

根據防洪要求，上游分區(qū)的洪量與下游設計斷面洪量重現(xiàn)期（設計頻率）相同。對于本研究，同頻率法可以有以下的兩種組合形式，分別為：當花園口斷面發(fā)生重現(xiàn)期為T（頻率為P）的洪水zp時，小浪底水庫也發(fā)生重現(xiàn)期為T的洪水xp，小-花區(qū)間發(fā)生相應洪水y，有y=zp-xp；當花園口斷面發(fā)生重現(xiàn)期為T（頻率為P）的洪水zp時，小-花區(qū)間也發(fā)生重現(xiàn)期為T的洪水yp，小浪底水庫發(fā)生相應洪水x，即x=zp-yp。

2.4.3 最可能洪水地區(qū)組成

針對傳統(tǒng)洪水地區(qū)組成算法缺乏對空間相關性的考慮，本文利用GH-Copula 函數(shù)構建研究區(qū)內小浪底洪水X與小-花區(qū)間洪水Y之間的聯(lián)合分布，其中，各個斷面獨立選樣并計算邊緣分布，以此推求在花園口站發(fā)生一定量級的洪水時，小浪底和小-花區(qū)間最可能的設計洪水過程。

利用GH-Copula 函數(shù)建立小浪底水庫洪水X和小-花區(qū)間洪水Y的聯(lián)合分布如下式：

式中：FX（x）和FY（y）分別為小浪底水庫洪水與小-花區(qū)間洪水的P-III型分布函數(shù)。

由式（19）可得聯(lián)合密度函數(shù)如下：

當下游花園口斷面發(fā)生設計洪水zp時，尋求上游小浪底和小-花區(qū)間洪水過程的最可能組合，即在滿足x+y=zp的條件下，求解概率密度函數(shù)f（x，y）的最大值點。將y=zp-xp代入f（x，y）中，可式（19）轉化為單變量x的函數(shù)，如下式：

上式對x求導可得：

令df/dx，整理得：

式中：αx、βx、a0.x和αy、βy、a0，y分別為小浪底與小花區(qū)間洪水要素的P-III分布FX（x）和FY（y）對應參數(shù)。

對上述的方程求解，可以得到在給定的花園口站設計洪水值zp時，小-花區(qū)間洪水與小浪底的最可能的洪水地區(qū)組合（x，y）。

2.4.4 構建Copula-GA算法

對于已經構建的基于Copula 洪水地區(qū)組成模型，求得最有可能的組合情況，即在整個變量序列里求得Copula 密度函數(shù)的極值所對應的設計洪量值。但是由于Copula 函數(shù)特性，該模型求解過程繁瑣，不能有效和準確的推求其全局最有可能的極值組合。因此，本文將遺傳算法（Genetic Algorithms）對與構建的基于Copula 洪水地區(qū)組成模型相耦合，以探尋既客觀又高效的地區(qū)組成結果［28-30］。

此處所構建的Copula-GA（Copula-Genetic Algorithms）模型，以式（23）作為遺傳算法的目標函數(shù)，即適應度函數(shù)，計算個體適應度：

Copula-GA模型的運行過程可概述如下：

（1）建立Copula洪水地區(qū)組成模型；

（2）確定目標函數(shù)f（x），即率密度函數(shù)式（23）；

（3）初始化種群；

（4）隨機生成初始種群；

（5）利用式（23）計算個體適應度，并保留適應度值較高的染色體；

（6）對種群交叉、變異得到下一代群體；

（7）重復（5）～（6），滿足終止條件后輸出適應度最大的染色體。

3 結果分析

3.1 花園口站設計洪水

本文的研究區(qū)屬于半濕潤區(qū)，暴雨強度大，一次暴雨歷時2～3 d，最長達5 d，洪水過程歷時長，可達10～12 d，因此本文構建年最大12 d 洪量W12與年最大洪峰流量Q的聯(lián)合分布。Q與W12的邊緣分布函數(shù)分別為FQ（q）和FW（w12）的P-III 型分布。使用線性矩法計算FQ（q）、FW（w12）的參數(shù)，成果如表1所示。

表1 花園口站洪水統(tǒng)計特征值成果Tab.1 The results of flood statistics characteristic value of Huayuankou Hydrological Station

由花園口站實測流量資料，計算得花園口站洪峰流量和12天洪量的Kendall 秩相關系數(shù)τ= 0.776，對應θ≈4.464。使用聯(lián)合變量同頻率法以及單變量同頻率法對花園口站設計洪峰流量以及12天洪量進行推求，成果見圖3、4以及表2所示。

圖3 花園口站洪峰流量和洪量聯(lián)合頻率分布圖Fig.3 Joint probability diagram of flood peak and flood volume at Huayuankou Station

由表2可知，采用聯(lián)合變量同頻率法，計算得到的各個重現(xiàn)期的設計洪峰流量Q和設計洪量W12大于通過單變量同頻率法計算獲得的設計洪峰流量Q和設計洪量W12。即使傳統(tǒng)單變量法能夠單獨控制變量的重現(xiàn)期，但是缺乏對變量之間相關性的考慮，然而聯(lián)合變量重現(xiàn)期所得的設計值不會低于設計標準。鑒于對防洪尤為不利的情況，取聯(lián)合變量重現(xiàn)期作為本次推求花園口站的設計洪水特征值，將會加大花園口洪水威脅程度，而由此設計結果優(yōu)化得到的調度方案會提高花園口站整體抵御洪水的能力，從而降低花園口站遭受洪災的風險，按變倍比放大獲得花園口站設計洪水過程線，結果如圖5所示。

圖5 花園口站設計洪水過程線Fig.5 Design flood process of Huayuankou Hydrological Station

表2 花園口水文站設計洪峰和洪量設計成果Tab.2 Design results of flood peak and flood volume of Huayuankou Hydrological Station

3.2 設計洪水地區(qū)組成

小浪底壩址、小-花區(qū)間采用1955-2015年連續(xù)還原天然狀態(tài)下洪水系列資料，通過矩法估計獲得小浪底水庫以及小-花區(qū)間洪水統(tǒng)計參數(shù)，見表3。

表3 小浪底和小-花區(qū)間洪水統(tǒng)計參數(shù)Tab.3 Results of flood statistics in Xiaolangdi and Xiao-hua interval

小浪底壩址和小-花區(qū)間及花園口的12天洪量邊緣分布分別采用P-III型分布，參數(shù)估計結果如表4所示。

表4 花園口、小浪底以及小-花區(qū)間P-III型分布參數(shù)估計結果Tab.4 Estimation results of P-III distribution parameters in Huayuankou Hydrological Station，Xiaolangdi and Xiao-hua interval

根據現(xiàn)有資料系列，分別采用同頻率地區(qū)組成法和基于GH-Copula函數(shù)和遺傳算法構建的Copula-GA 算法對小浪底與小-花區(qū)間洪水的地區(qū)組成進行計算。

（1）同頻率地區(qū)組成法通過假設小浪底洪水重現(xiàn)期與花園口相同推求小浪底設計洪水特征值，再利用干支流關系計算小-花區(qū)間洪水特征值，計算結果如表5所示。

表5 小浪底與花園口同頻率地區(qū)組成計算結果Tab.5 Flood region composition results of Xiaolangdi and Huayuankou with the same frequency

（2）基于Copula 函數(shù)建立小浪底和小-花區(qū)間聯(lián)合分布，并推求其Copula 密度函數(shù)，利用干支流洪量組合關系，將在特定重現(xiàn)期下的Copula 密度函數(shù)轉變?yōu)閱巫兞棵芏群瘮?shù)。然后采用GA 算法尋優(yōu)求解，以推求其密度函數(shù)最大極值的可能組合，并推求小浪底與小-花區(qū)間12 天洪量的對應值，得到花園口站設計洪水最可能地區(qū)組成，成果見圖6以及表6。

圖4 花園口站洪峰和洪量聯(lián)合重現(xiàn)期概率圖Fig.4 Probability Diagram of Joint Return Period of Flood Peak and Flood Volume at Huayuankou Hydrological Station

圖6 小浪底和小-花區(qū)間聯(lián)合分布函數(shù)圖及小浪底洪量概率密度圖（P=2%）Fig.6 Xiaolangdi and Xiao-Hua interval joint distribution function diagram and Xiaolangdi flood probability density diagram（P=2%）

通過表6 可以看出，隨著花園口洪量的增大，小浪底與小-花區(qū)間洪量量級都相應的隨之增大，并且小浪底洪量所占的比例明顯加大。由此可以得出：小浪底的洪水對花園口洪水的發(fā)生起到了決定性作用，如果僅發(fā)生小-花區(qū)間為主的洪水，花園口一般不會出現(xiàn)災害性的洪水；只有發(fā)生全流域型洪水時，即小-花區(qū)間與小浪底同時發(fā)生洪水時花園口才會出現(xiàn)災害性大洪水，這與黃河中下游洪水成因及來源相印證［17，31］。

表6 花園口站設計洪水最可能地區(qū)組成成果Tab.6 The most likely flood region composition results of the design flood at Huayuankou Station

由花園口站設計洪水最可能地區(qū)組成成果可知，將小浪底在各頻率地區(qū)組成所分配的洪量值，由小浪底洪量頻率曲線推求出其相對應的頻率值；然后再根據洪峰、洪量同頻率的這一假設，分別推求出與小浪底各洪量頻率對應相同的洪峰值，成果見表7。

表7 小浪底壩址設計洪水洪峰、洪量成果Tab.7 Design flood peak and flood volume results of Xiaolangdi dam site

對同頻率地區(qū)組成法和Copula-GA 法計算的小浪底的洪峰流量和洪量結果分析可知：當花園口發(fā)生的洪水重現(xiàn)期為200年（P=1%）、100年（P=1%）以及50年（P=2%）時，小浪底最可能發(fā)生的洪水重現(xiàn)期約為52 a（P=1.92%）、29 a（P=3.44%）和18 a（P=5.76%）；通過假設小浪底洪水重現(xiàn)期與花園口相同計算得到的小浪底設計洪水，相較Copula-GA 法計算得到的最可能的設計洪水地區(qū)組成，前者所推求的小浪底洪水特征值較大，進而導致小-花區(qū)間洪水洪水特征值較小，不能達到相應的設計頻率。

3.3 小浪底及小-花區(qū)間設計洪水過程

根據求得的小浪底對應花園口各設計頻率的洪峰、洪量設計值，放大所選取的典型洪水，推求小浪底52年一遇、29年一遇、18年一遇設計洪水過程，成果見圖7。

圖7 小浪底設計洪水過程線Fig.7 Xiaolangdi design flood process

根據小浪底設計洪水過程，由馬斯京根方法將小浪底洪水演進至花園口斷面。通過花園口各重現(xiàn)期設計洪水過程線逐時刻減去小浪底演進至花園口過程線，即可得到小-花區(qū)間設計洪水過程線。

4 結論

本文使用Copula 函數(shù)構建了花園口站洪峰流量與洪量的聯(lián)合分布，并計算得到了特定重現(xiàn)期的洪峰流量以及洪量，充分地考慮了洪峰流量與洪量之間的相關性，通過對比分析現(xiàn)有設計成果發(fā)現(xiàn)，兼顧洪峰和洪量相關性計算得到的各重現(xiàn)期下的洪水峰值流量和時段洪量設計值略高于現(xiàn)有的設計成果值［13］，更全面的考慮洪水要素之間的相互作用，也將更適用于受洪峰和洪量共同控制的區(qū)域。

在計算花園口站各重現(xiàn)期下的設計洪水地區(qū)組成時，針對傳統(tǒng)的設計洪水地區(qū)組成計算方法中存在問題，本文采用Copula 函數(shù)構建了小浪底與小-花區(qū)間12 天洪量的聯(lián)合分布，并與遺傳算法相耦合構建了Copula-GA 模型，以快速準確的推求Copula 密度函數(shù)的最大值組合，即設計洪水最可能地區(qū)組成，得到的結果更具備客觀性；當下游防洪斷面將流經特定重現(xiàn)期洪水時，通過本研究最可能的地區(qū)組成方法所推得上游各斷面洪水重現(xiàn)期，相比傳統(tǒng)同頻率方法，也更符合研究區(qū)洪水的實際情況。

洪水歷時作為另一個洪水特征量，并與洪峰流量、洪量存在一定的相關性［21］。本文僅考慮了洪峰流量、洪量的相關性，建立了洪峰流量、洪量的聯(lián)合分布，在后續(xù)研究中還應建立洪峰流量、洪量以及洪水歷時間的聯(lián)合分布，推求不同重現(xiàn)期下的洪峰流量、洪量、洪水歷時組合?！?/p>