，，

(武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072)

1 研究背景

枯水是天然的河流、湖泊每年都會發(fā)生的季節(jié)性現(xiàn)象，是水文過程中不可或缺的重要組成部分。但由于早期枯水對人類的影響較小，其研究水平遠不及洪水。隨著社會經(jīng)濟的發(fā)展，水資源供需矛盾日益突出，尤其在枯水季節(jié)，由于河流徑流量減少，加之人類過度開發(fā)利用枯水資源帶來了水質(zhì)污染加重、河流斷流、地面沉降等一系列問題，人們越發(fā)重視枯水問題[1]。對于干旱缺水地區(qū)，枯水的主要基本特征，更是供水工程規(guī)劃設(shè)計、水資源優(yōu)化配置和合理保護的重要依據(jù)[2]。

目前，大部分的枯水研究集中于河流枯水徑流的研究[1-6]。但對于處在有明顯干、濕季節(jié)的地區(qū)的天然湖泊，由于降雨的季節(jié)性差異，直接影響了地表和地下徑流的形成，因而湖泊與之對應(yīng)分為豐水期和枯水期，其水文特征迥異[7]，因而很有必要對其枯水特征開展研究。湖泊枯水的研究成果較為有限，如梁明川[7]對天然湖泊進行枯水期水文分析, 建立水均衡方程，然后利用水文、氣象資料計算流域地下水資源。譚國良等[8]針對近年來鄱陽湖頻繁出現(xiàn)的干旱現(xiàn)象，采用長系列枯水期資料分析了鄱陽湖區(qū)枯水特征及其成因?，F(xiàn)有工作均未全面揭示湖泊枯水特征變量的概率統(tǒng)計特征。

枯水并不等同于干旱，連續(xù)的季節(jié)性枯水事件并不一定會形成干旱，但由于它們在發(fā)生成因和演變規(guī)律上存在一定的聯(lián)系且相互影響，因而在研究方法上可以相互借鑒。本文基于干旱分析中廣泛采用的游程理論識別湖泊枯水，提取湖泊枯水特征變量，然后采用Copula函數(shù)進行湖泊枯水特征變量的多變量頻率分析，為鄱陽湖區(qū)枯水季節(jié)的水資源管理提供理論參考。

2 研究區(qū)域與方法

鄱陽湖納贛、撫、信、饒、修等5大河及湖區(qū)區(qū)間來水，經(jīng)湖盆調(diào)蓄后由湖口注入長江。流域水網(wǎng)比較密集，鄱陽湖同時受到長江干流和五河的影響，由于江、湖、河枯期時間基本重疊，在一定程度上延長了湖區(qū)枯水持續(xù)時間，加劇了枯水程度，對于水資源管理極為不利[8]。近年來，鄱陽湖流域正面臨著嚴重的干旱問題，給鄱陽湖流域和長江中下游地區(qū)造成了重大的經(jīng)濟損失。因此，正確認識鄱陽湖湖區(qū)枯水規(guī)律，是指導(dǎo)抗旱、保護生態(tài)、開發(fā)利用水資源和保證流域社會經(jīng)濟協(xié)調(diào)發(fā)展的重要準備工作。

本文選取都昌水文站作為鄱陽湖區(qū)代表水文站。都昌站地處鄱陽湖水域中間，其水位能較準確地代表鄱陽湖的水位數(shù)據(jù)，能反映出鄱陽湖的整體水情[9]。本文采用的數(shù)據(jù)為都昌水文站1958—2007年長系列日平均水位資料，數(shù)據(jù)來自江西省水文局。

2.1 枯水識別與枯水特征變量提取

一般而言，枯水期的起訖時間和歷時，主要取決于水體的補給狀況的確定，但實際操作上存在一定困難。本文采用在干旱分析中廣泛采用的游程理論[10]，進行湖泊枯水的識別，并提取枯水特征變量。游程理論一個重要前提是閾值的選取，本文作者認為，在一個完整的水文年內(nèi)(枯水分析常用每年7月至次年6月[9])，當湖區(qū)代表站點日平均水位低于事先確定的閾值時，湖泊即進入枯水期。閾值可根據(jù)水位數(shù)據(jù)的分位數(shù)來確定，本文采用文獻[9]中建議的鄱陽湖都昌站枯水劃分標準，即都昌水文站水位12.8 m作為鄱陽湖進入枯水期的水位閾值。枯水事件的游程分析圖如圖1所示。

圖1 湖泊枯水事件游程圖

在基于日資料進行枯水統(tǒng)計的過程中，可能出現(xiàn)水位短暫地低于閾值，而后迅速恢復(fù)至閾值以上且量級較小的枯水過程，使得統(tǒng)計變得復(fù)雜，如圖1(b)區(qū)中箭頭所示。綜合參考陸桂華等[11]，Madsen和Rosbjerg[12]提出的方法，若某一段枯水過程的最大枯水強度小于造成明顯取用水困難的臨界強度(本文取1.0 m，對應(yīng)水位11.8 m[9])，歷時小于0.1E(D)時，該枯水過程忽略不計，E(D)指枯水歷時的期望。

此外，當2個不能忽略的枯水過程中間還摻雜著若干可以忽略的枯水過程，如圖1(d)區(qū)所示，此時枯水歷時D4=t8-t7，相應(yīng)得到平均枯水強度。

2.2 基于Copula函數(shù)的枯水多變量聯(lián)合分布

2.2.1 Copula函數(shù)的性質(zhì)

一個完整的水文事件往往包含多個特征變量，且特征變量之間通常存在相關(guān)性，Copula理論正是通過構(gòu)造聯(lián)合分布來描述這種相關(guān)結(jié)構(gòu)的[14]。Copula是定義域為[0,1]均勻分布的多維聯(lián)合分布函數(shù)，它可以將多個隨機變量的邊緣分布連接起來得到它們的聯(lián)合分布。

Sklar定理：令H為一個n-維分布函數(shù)，其邊緣分布為F1,F2,…,Fn，則存在一個n-Copula函數(shù)C，使得對任意x∈Rn [15]，即

H(x1,x2,…,xn)=

Cθ(F1(x1),F2(x2),…,Fn(xn)) 。

(1)

式中θ為Copula函數(shù)的參數(shù)。

2.2.2 Archimedean Copula函數(shù)

Archimedean Copula函數(shù)結(jié)構(gòu)簡單，計算簡便，可以構(gòu)造出多種形式多樣、適應(yīng)性強的多變量聯(lián)合分布函數(shù)，在水文分析計算中得到了廣泛的應(yīng)用[14,16-18]，根據(jù)構(gòu)造方式的不同可以分為對稱型和非對稱型2種。

基于前述方法識別枯水事件，分別提取枯水歷時D、平均枯水強度AveI和最大枯水強度MaxI作為湖泊枯水特征變量，分別確定其邊緣分布函數(shù)，然后通過三維非對稱型Copula函數(shù)構(gòu)造其聯(lián)合分布。本文選用在水文領(lǐng)域較為常用的3種Archimedean Copula函數(shù)，即Gumbel-Hougaard，Clayton，F(xiàn)rank Copula函數(shù)，3種Copula函數(shù)的三維非對稱形式如下。

Gumbel-Hougaard Copula函數(shù)：

C(u1,u2,u3)=exp{-[((-lnu1)θ2+

(-lnu2)θ2)θ1/θ2+(-lnu3)θ1]1/θ1},

1≤θ1≤θ2。

(2)

Clayton Copula函數(shù)：

C(u1,u2,u3)=

[(u1-θ2+u2-θ2-1)θ1/θ2+u3-θ1-1]-1/θ1,

1≤θ1≤θ2。

(3)

Frank Copula函數(shù)：

C(u1,u2,u3)=

(4)

式中：u1,u2,u3分別代表邊緣分布函數(shù);參數(shù)θ1和θ2通過極大似然法進行估計。

2.2.3 Copula函數(shù)的選擇

擬合優(yōu)度評價指標是選擇分布線型的一個重要標準[14]。采用離差平方和最小準則(OLS)來評價不同Copula函數(shù)的擬合效果，選擇最合適的Copula函數(shù)，其表達式為

(5)

式中：n為樣本容量;Pei,Pi分別為經(jīng)驗頻率和理論頻率。

2.3 枯水事件的重現(xiàn)期分析

在水文分析中，確定某一稀遇事件或特大災(zāi)害發(fā)生的概率，具有十分重要的現(xiàn)實意義。工程中常用重現(xiàn)期來度量風(fēng)險的大小。

枯水事件的單變量重現(xiàn)期為：

(6)

式中：TD,TAvel,TMaxI分別為枯水歷時、平均枯水強度、最大枯水強度的單變量重現(xiàn)期；FD(d),FAvel(avei),FMaxI(maxi)分別為枯水歷時、平均枯水強度、最大枯水強度的邊緣分布函數(shù)。

一場完整的枯水事件包括枯水歷時、平均枯水強度、最大枯水強度等多個變量，僅考慮某一個的重現(xiàn)期很可能高估或低估枯水事件出現(xiàn)的風(fēng)險，因而必須分析多個枯水變量之間的相關(guān)性及其出現(xiàn)的聯(lián)合重現(xiàn)期。

(7)

在給定變量Xk≤xk的前提下，事件Xi≤xi,Xj≤xj出現(xiàn)的條件概率為

(8)

對應(yīng)的條件重現(xiàn)期的表達式為

(9)

3 結(jié)果與分析

3.1 邊緣分布的建立

首先對提取的枯水特征變量進行“三性”檢驗，各變量均符合獨立同分布假設(shè)。采用指數(shù)分布、Gamma分布、P-III分布、Weibull分布、Logistic分布和廣義極值(GEV)分布分別擬合枯水歷時、平均枯水強度和最大枯水強度，采用線性矩法估計分布函數(shù)的參數(shù)，對分布擬合效果進行K-S檢驗，并以前述離差平方和(OLS)公式(5)最小為準則比較確定最適宜擬合各變量的線型。

經(jīng)計算，選用Weibull分布擬合干旱歷時，Logistic分布擬合平均枯水強度，GEV分布擬合最大枯水程度。經(jīng)驗頻率和理論分布的擬合結(jié)果如圖2所示。擬合效果均通過了95%置信水平的K-S檢驗。從圖2中可以看出，各種理論分布的擬合效果均較理想。

圖2 理論邊緣分布擬合

3.2 相關(guān)性分析

分別采用Pearson相關(guān)系數(shù)和Kendall秩相關(guān)系數(shù)計算3個變量的兩兩相關(guān)性，結(jié)果如表1所示。表中對角線上方元素為Kendall秩相關(guān)系數(shù)，下方元素為Pearson相關(guān)系數(shù)。從表中可知，各變量之間呈現(xiàn)明顯的正相關(guān)性，平均枯水強度和最大枯水強度之間的相關(guān)性最高，枯水歷時和最大枯水強度次之，枯水歷時和平均枯水強度最弱。

3.3 聯(lián)合分布的建立

應(yīng)用極大似然法分別估計3種Copula函數(shù)的參數(shù)，構(gòu)建3個變量聯(lián)合分布，并計算各自理論頻率與經(jīng)驗頻率的離差平方和，相關(guān)結(jié)果如表2所示。

表1 枯水變量Pearson和Kendall秩相關(guān)系數(shù)

表2 Copula函數(shù)參數(shù)估計及OLS

表3 聯(lián)合概率(重現(xiàn)期)和條件概率結(jié)果統(tǒng)計

從表2中可以看出Clayton Copula的擬合效果優(yōu)于其他2個連接函數(shù)，其OLS為0.034 7。

圖3 理論分布與經(jīng)驗點據(jù)擬合

根據(jù)Genest和Rivest[19]提出的一種比較直觀地評價Copula 函數(shù)的方法，分別計算理論概率估計值Kc(t)和經(jīng)驗頻率估計值Ke(t)，然后點繪Kc(t)和Ke(t)關(guān)系圖，如果兩者相關(guān)系數(shù)較高且圖上的點均勻地落在45°線附近，那么表明Kc(t)和Ke(t) Copula函數(shù)擬合得較好。點繪Clayton Copula函數(shù)計算得到的理論分布與經(jīng)驗頻率關(guān)系圖如圖3所示。從圖3可以看出，理論和經(jīng)驗分布的相關(guān)系數(shù)為0.990 3，且點據(jù)均勻地落在45°線附近，可知構(gòu)造的Copula函數(shù)是合理的。

3.4 枯水概率和重現(xiàn)期分析

根據(jù)建立的聯(lián)合分布，可以推求枯水事件的聯(lián)合重現(xiàn)期和給定其中一個變量時的條件重現(xiàn)期。單變量的重現(xiàn)期T由邊緣分布確定。根據(jù)聯(lián)合概率，由式(7)計算聯(lián)合重現(xiàn)期TJ，結(jié)果均列于表3。由表3可知：3個變量的聯(lián)合重現(xiàn)期隨著單變量的重現(xiàn)期增加而增加，單變量的重現(xiàn)期要高于3個變量的聯(lián)合重現(xiàn)期。以邊緣分布函數(shù)F(x)等于0.99時為例，此時單變量重現(xiàn)期為100 a，而聯(lián)合重現(xiàn)期為34.47 a，二者差異明顯。

圖4 枯水事件的條件重現(xiàn)期

采用式(9)，計算分別給定枯水歷時、平均枯水強度和最大枯水強度的單變量重現(xiàn)期為5 a和50 a時，事件E={Xi>xiorXj>xj|Xk≤xk}的條件重現(xiàn)期，結(jié)果如圖4所示。由圖4可知：給定枯水歷時的單變量重現(xiàn)期時計算得到的條件重現(xiàn)期低于給定另外2個變量時的條件重現(xiàn)期。在給定枯水事件任一單變量的重現(xiàn)期時，計算得到的條件重現(xiàn)期隨著另2個變量的單變量重現(xiàn)期的增大而增大。隨著給定的單變量重現(xiàn)期的增大，另2個變量量級相同時的枯水事件條件重現(xiàn)期逐漸減小。

以上結(jié)果均表明構(gòu)造的聯(lián)合分布能提供更多的湖泊枯水信息，能更全面地反映枯水事件變量間的相關(guān)性特征。

4 結(jié) 語

研究湖泊枯水規(guī)律，對于湖泊枯水季節(jié)的管理具有重要意義。本文借鑒干旱分析中應(yīng)用廣泛的游程理論識別湖泊枯水，從鄱陽湖湖區(qū)都昌水文站長系列日平均水位資料中提取枯水歷時、平均枯水強度和最大枯水強度作為湖泊枯水特征變量，確定其邊緣分布函數(shù)，采用Archimedean Copula函數(shù)構(gòu)造了變量的聯(lián)合分布，依據(jù)離差平方和最小準則選用最優(yōu)的Copula函數(shù)，分析枯水事件出現(xiàn)的概率，計算了枯水事件的聯(lián)合重現(xiàn)期和條件重現(xiàn)期，結(jié)果表明：

(1) 建議選擇Weibull分布、Logistic分布和GEV分布，分別作為邊緣分布函數(shù)擬合鄱陽湖都昌站干旱歷時、平均枯水強度和最大枯水程度。

(2) 依據(jù)離差平方和最小準則優(yōu)選的Clayton Copula函數(shù)，適于用來擬合鄱陽湖都昌站枯水多變量聯(lián)合分布。

(3) 聯(lián)合分布考慮了枯水事件的多個變量，能提供更多的湖泊枯水特征信息，更能全面揭示湖泊枯水的概率統(tǒng)計特征。

參考文獻：

[1] 殷福才, 王在高, 梁 虹. 枯水研究進展[J]. 水科學(xué)進展, 2004, 15(2): 249-254. (YIN Fu-cai, WANG Zai-gao, LIANG Hong. Advances in Low Flow Research [J]. Advances in Water Science, 2004, 15(2): 249-254. (in Chinese))

[2] 馮國章, 王雙銀. 河流枯水流量特征研究[J]. 自然資源學(xué)報, 1995, 10(2): 127-135. (FENG Guo-zhang, WANG Shuang-yin. A Study of the Characteristics of the Low Flows of Rivers [J]. Journal of Nature Resources, 1995, 10(2): 127-135. (in Chinese))

[3] SMAKHTIN V U. Low Flow Hydrology: A Review [J]. Journal of Hydrology, 2001, 240(3): 147-186.

[4] GOTTSCHALK L, YU K, LEBLOIS E,etal. Statistics of Low Flow: Theoretical Derivation of the Distribution of Minimum Streamflow Series [J]. Journal of Hydrology, 2013, 481: 204-219.

[5] 周 芬, 郭生練, 熊立華, 等. 枯水頻率分析線型的比較研究[J]. 水文, 2006, 26(1): 28-33. (ZHOU Fen, GUO Sheng-lian, XIONG Li-hua,etal. Comparative Study on Distribution Functions for Low Flow [J]. Journal of China Hydrology, 2006, 26(1): 28-33. (in Chinese))

[6] 方 彬, 郭生練, 肖 義, 等. 基于低定量取樣的枯水頻率分析[J]. 武漢大學(xué)學(xué)報 (工學(xué)版),2006,39(2):1-4. (FANG Bin, GUO Sheng-lian, XIAO Yi,etal. Low Flow Frequency Analysis Based on LUT Sampling [J]. Engineering Journal of Wuhan University,2006,39(2):1-4. (in Chinese))

[7] 梁明川. 天然湖泊枯水期水文分析及流域地下水資源計算[J]. 水文地質(zhì)工程地質(zhì), 1984, (4): 58-59. (LIANG Ming-chuan. Hydrology Analysis and Groundwater Resource Calculation of Nature Lake during Dry Season[J]. Hydrogeology & Engineering Geology, 1984, (4): 58-59. (in Chinese))

[8] 譚國良, 郭生練, 王 俊, 等. 鄱陽湖生態(tài)經(jīng)濟區(qū)水文水資源演變規(guī)律研究[M]. 北京: 中國水利水電出版社, 2013. (TAN Guo-liang, GUO Sheng-lian, WANG Jun,etal. Hydrology and Water Resource Variation Law of Poyang Lake Basin [M]. Beijing: China Water Power Press, 2013. (in Chinese))

[9] 閔 騫, 占臘生. 1952—2011年鄱陽湖枯水變化分析[J]. 湖泊科學(xué), 2012, (5): 675-678. (MIN Qian, ZHAN La-sheng. Characteristics of Low-water Level Changes in Lake Poyang during 1952-2011 [J]. Journal of Lake Sciences, 2012, (5): 675-678. (in Chinese))

[10] YEVJEVICH V. An Objective Approach to Definitions and Investigations of Continental Hydrologic Droughts [M]. Fort Collins: Colorado State University, 1967.

[11] 陸桂華, 閆桂霞, 吳志勇, 等. 基于copula函數(shù)的區(qū)域干旱分析方法[J]. 水科學(xué)進展,2010,21(2):188-193. (LU Gui-hua, YAN Gui-xia, WU Zhi-yong,etal. Regional Drought Analysis Approach Based on Copula Function [J]. Advances in Water Science,2010,21(2): 188-193. (in Chinese))

[12] MADSEN H, ROSBJERG D. On the Modeling of Extreme Droughts [J]. IAHS Publications-Series of Proceedings and Reports-Intern Assoc Hydrological Sciences, 1995, 231: 377-386.

[13] TALLAKSEN L M, MADSEN H, CLAUSEN B. On the Definition and Modeling of Streamflow Drought Duration and Deficit Volume [J]. Hydrological Sciences Journal, 1997, 42(1): 15-33.

[14] 郭生練, 閆寶偉, 肖 義, 等. Copula 函數(shù)在多變量水文分析計算中的應(yīng)用及研究進展[J]. 水文, 2008, 28(3): 1-7. (GUO Sheng-lian, YAN Bao-wei, XIAO Yi,etal. Multivariate Hydrological Analysis and Estimation [J]. Journal of China Hydrology,2008,28(3):1-7. (in Chinese))

[15] NELSEN R B. An Introduction to Copulas [M]. USA: Springer, 2006.

[16] 閆寶偉, 郭生練, 肖 義, 等. 基于兩變量聯(lián)合分布的干旱特征分析[J]. 干旱區(qū)研究,2007,24(4):537-542. (YAN Bao-wei, GUO Sheng-lian, XIAO Yi,etal. Analysis on Drought Characteristics Based on Bivariate Joint Distribution [J]. Arid Zone Research, 2007, 24(4): 537-542. (in Chinese))

[17] 閆寶偉, 郭生練, 肖 義. 南水北調(diào)中線水源區(qū)與受水區(qū)降水豐枯遭遇研究[J]. 水利學(xué)報, 2007,38(10):1178-1185. (YAN Bao-wei, GUO Sheng-lian, XIAO Yi,etal. Synchronous-asynchronous Encounter Probability of Rich-poor Precipitation between Water Source Area and Water Receiving Areas in the Middle Route of South-to-North Water Transfer Project [J]. Journal of Hydraulic Engineering, 2007,38(10):1178-1185. (in Chinese))

[18] CHEN L, SINGH V P, GUO S,etal. Drought Analysis Using Copulas [J]. Journal of Hydrologic Engineering, 2013, 18: 797-808.

[19] GENEST C, RIVEST L P. Statistical Inference Procedures for Bivariate Archimedean Copulas [J]. Journal of the American Statistical Association, 1993, 88(423): 1034-1043.