基于Copula函數(shù)的鄱陽(yáng)湖都昌站枯水多變量頻率分析

，，

(武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

1 研究背景

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

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

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

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

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

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

2.1 枯水識(shí)別與枯水特征變量提取

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

圖1 湖泊枯水事件游程圖

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

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

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

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

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

Sklar定理：令H為一個(gè)n-維分布函數(shù)，其邊緣分布為F1,F2,…,Fn，則存在一個(gè)n-Copula函數(shù)C，使得對(duì)任意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)簡(jiǎn)單，計(jì)算簡(jiǎn)便，可以構(gòu)造出多種形式多樣、適應(yīng)性強(qiáng)的多變量聯(lián)合分布函數(shù)，在水文分析計(jì)算中得到了廣泛的應(yīng)用[14,16-18]，根據(jù)構(gòu)造方式的不同可以分為對(duì)稱型和非對(duì)稱型2種。

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

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通過(guò)極大似然法進(jìn)行估計(jì)。

2.2.3 Copula函數(shù)的選擇

擬合優(yōu)度評(píng)價(jià)指標(biāo)是選擇分布線型的一個(gè)重要標(biāo)準(zhǔn)[14]。采用離差平方和最小準(zhǔn)則(OLS)來(lái)評(píng)價(jià)不同Copula函數(shù)的擬合效果，選擇最合適的Copula函數(shù)，其表達(dá)式為

(5)

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

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

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

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

(6)

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

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

(7)

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

(8)

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

(9)

3 結(jié)果與分析

3.1 邊緣分布的建立

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

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

圖2 理論邊緣分布擬合

3.2 相關(guān)性分析

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

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

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

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

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

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

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

圖3 理論分布與經(jīng)驗(yàn)點(diǎn)據(jù)擬合

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

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

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

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

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

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

4 結(jié) 語(yǔ)

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

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

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

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

參考文獻(xiàn)：

[1] 殷福才, 王在高, 梁 虹. 枯水研究進(jìn)展[J]. 水科學(xué)進(jìn)展, 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] 馮國(guó)章, 王雙銀. 河流枯水流量特征研究[J]. 自然資源學(xué)報(bào), 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é)報(bào) (工學(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ì)算[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] 譚國(guó)良, 郭生練, 王 俊, 等. 鄱陽(yáng)湖生態(tài)經(jīng)濟(jì)區(qū)水文水資源演變規(guī)律研究[M]. 北京: 中國(guó)水利水電出版社, 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年鄱陽(yáng)湖枯水變化分析[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é)進(jìn)展,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ù)在多變量水文分析計(jì)算中的應(yīng)用及研究進(jìn)展[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é)報(bào), 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.