基于Copula的汾河上游水文干旱頻率分析

任　璐，趙雪花

基于Copula的汾河上游水文干旱頻率分析

任璐，趙雪花

（太原理工大學(xué)水利科學(xué)與工程學(xué)院，山西太原 030024）

運用游程理論對汾河上游4個水文站的月徑流資料進行干旱識別；選擇4種分布函數(shù)擬合干旱特征變量，通過Kolmogorov-Smirnov檢驗法優(yōu)選單變量邊緣分布；利用Copula函數(shù)建立干旱特征變量的二維聯(lián)合分布，計算重現(xiàn)期。結(jié)果表明，汾河上游干旱特征變量的最佳邊緣分布是對數(shù)正態(tài)分布；汾河上游大多數(shù)干旱重現(xiàn)期小于5年，且干旱歷時與烈度的單變量重現(xiàn)期與二維重現(xiàn)期較為接近；干旱歷時與烈度峰值、烈度與烈度峰值的單變量重現(xiàn)期與二維重現(xiàn)期差距較大。

水文干旱；頻率分析；Copula函數(shù)；月徑流；汾河上游

0　引言

干旱是一種復(fù)雜、多因素并具有全球影響力的自然災(zāi)害［1］。它經(jīng)常發(fā)生卻甚少為人了解。我國屬于東亞季風(fēng)氣候，年季間季風(fēng)的不穩(wěn)定性是造成我國干旱頻發(fā)的一個主要原因［2］。氣候變化和人類活動的共同作用，導(dǎo)致干旱的發(fā)生頻率不斷上升、造成損失不斷加大，嚴(yán)重制約了我國社會經(jīng)濟的發(fā)展。緩解現(xiàn)今嚴(yán)峻的干旱形勢迫在眉睫，因此對于干旱事件的頻率分析是極其必要的。

最初的多變量干旱頻率的研究建立在一系列假定的基礎(chǔ)上。即，各變量的邊緣分布必須服從同種分布［3-6］或者聯(lián)合分布必須采用正態(tài)分布和經(jīng)過數(shù)學(xué)轉(zhuǎn)換得到的正態(tài)分布［6-8］。然而，實際情況往往不能滿足以上假定，難以應(yīng)用到干旱頻率分析中。隨后出現(xiàn)的非參數(shù)研究方法［9-10］雖然可避免以上假定，但是其數(shù)學(xué)推導(dǎo)及參數(shù)計算都太過復(fù)雜，在實際應(yīng)用中難以推廣。相比以上幾種方法，Copula函數(shù)法既可以構(gòu)建服從不同邊緣分布的特征變量的聯(lián)合分布函數(shù)［11］，又便于計算，更符合干旱頻率分析的實際情況；因此可用于水文干旱頻率分析中。由于水文現(xiàn)象具有很明顯的區(qū)域性，不同地區(qū)水文情況不同，干旱特征變量的邊緣分布類型也不盡相同。部分研究者直接擬定單變量邊緣分布類型，可能會給后續(xù)研究造成不必要的誤差。前人對于水文干旱頻率的研究多集中于廣東、新疆等省份或者黑河、鄱陽湖等流域，鮮有對汾河上游的水文干旱頻率的研究。因此，本文先選擇4種水文分析中常見的分布對汾河上游的干旱特征變量進行擬合，再從中優(yōu)選出適合各特征變量的邊緣分布，在此基礎(chǔ)上利用Copula函數(shù)構(gòu)建二維干旱變量的聯(lián)合分布函數(shù)，并計算從干旱序列中提取出的代表數(shù)對應(yīng)的單變量重現(xiàn)期、二維聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期，以此來分析汾河上游的水文干旱現(xiàn)象。

本文選取汾河上游的4個典型水文站：上靜游（1955年 ～2012年）、汾河水庫（1958年 ～2000年）、寨上（1956年 ～2000年）、蘭村（1951年 ～2012年）的月徑流資料，對汾河上游的水文干旱情況進行分析，以為該區(qū)域今后的農(nóng)業(yè)發(fā)展、工業(yè)生產(chǎn)以及防旱抗旱工作提供科學(xué)依據(jù)。

1　研究方法

1.1干旱事件的定義與識別

利用游程理論［12］提取干旱特征變量：選取月徑流距平百分比為-10%作為截取水平R1，當(dāng)徑流序列Ri（t=1，2，…，N）在一個或多個時段內(nèi)連續(xù)小于R1時，出現(xiàn)負游程，即認為發(fā)生干旱；負游程的長度作為干旱歷時D（月）；負游程的陰影部分面積作為干旱烈度S（萬m3）；負游程的極值作為干旱烈度峰值M（萬m3）。

實際分析中常出現(xiàn)這樣一種情況：一些歷時很短的非干旱事件夾雜在一場嚴(yán)重干旱事件之間，導(dǎo)致一場嚴(yán)重干旱被分隔為多場一般干旱，不能反映干旱的實際情況，因此需要按照一定的準(zhǔn)則和方法對這些干旱進行合并［13］。取R0為徑流距平百分比為0時所對應(yīng)的徑流值，R2為徑流距平百分比為 -50%時所對應(yīng)的徑流值。如圖1所示，R0、R1、R2將徑流量劃分為4部分，則按以下方法判斷是否發(fā)生干旱：①當(dāng)徑流序列中出現(xiàn)單月徑流量R＜R2時，認為這個月為一場干旱事件（如圖中AB段）；②當(dāng)徑流序列中出現(xiàn)單月徑流量R2＜R＜R1時，并不認為這個月是一次干旱事件（如圖中IJ段）；③兩次干旱事件之間只間隔一個月且該月的月徑流量R＜R0時，認為這兩個干旱事件有關(guān)聯(lián)，應(yīng)合并為一場干旱事件（如圖中CD段和EF段應(yīng)合并為一場干旱事件，干旱歷時D2=DCD+DDE+DEF，干旱烈度S= S2+S2’）；④當(dāng)兩次干旱事件之間只間隔一個月且該月的月徑流量R＞R0時，認為其仍為兩場干旱事件（如圖中EF段和GH段不需要合并為一場干旱）。

圖1　游程理論識別干旱變量示意

1.2干旱特征變量的邊緣分布

本文選取對數(shù)正態(tài)分布（兩參數(shù)，即2P，下同）、指數(shù)分布（2P）、伽馬分布（2P）和廣義帕累托分布（3P）4種分布對汾河上游的4個水文站的干旱歷時、烈度和烈度峰值進行擬合，具體表達式如下：

式中，α為尺度參數(shù)；ε為位置參數(shù)；k為形狀參數(shù)；erf為誤差函數(shù)，即為伽馬函數(shù)

采用矩法進行參數(shù)估計，利用 Kolmogorovmirnov檢驗法（K-S法）［14］進行擬合優(yōu)度檢驗，得到最優(yōu)的單變量邊緣分布函數(shù)。

1.3二維干旱特征變量模型

1.3.1Copula函數(shù)定義

Copula可以連接一維的邊緣分布，形成在［0，1］上的多元分布的函數(shù)［14］。設(shè)X，Y為連續(xù)的隨機變量，其邊緣分布函數(shù)為FX和FY，F(xiàn)（x，y）為變量X和Y的聯(lián)合分布函數(shù)，則存在唯一的Copula函數(shù)C，使得

式中，θ為待定參數(shù)。

1.3.2二維干旱特征變量的函數(shù)選擇、參數(shù)估計及擬合度檢驗

本文選用Clayton Copula、Gumbel-Hougaard（GH）Copula、Ali-Mikhail-Haq（AMH）Copula和 Nelson NO.2 Copula 4種Copula函數(shù)擬合二維干旱特征變量的聯(lián)合分布。利用 Kendall秩相關(guān)系數(shù) τ和Archimedean Copula函數(shù)的參數(shù) θ的解析關(guān)系［14，15］，對參數(shù)θ進行估計。利用均方根誤差（Root Mean Square Error，RMSE）準(zhǔn)則和赤池信息量準(zhǔn)則（Akaike Information Criterion，AIC）［14］進行擬合優(yōu)度檢驗，依據(jù)是RMSE和AIC值最小，表達式：

式中，n為樣本容量；m為Copula參數(shù)估算數(shù)；Pc為Copula多元聯(lián)合分布計算值；P0為多元聯(lián)合分布經(jīng)驗值。

1.4干旱重現(xiàn)期計算

1.4.1單變量重現(xiàn)期

若干旱事件系列（N年）發(fā)生n次，則發(fā)生一次干旱事件的年數(shù)（重現(xiàn)期）［10］可記為

式中，F(xiàn)（x）為干旱特征變量的邊緣分布，N為研究站點的徑流資料長度，n為干旱事件次數(shù)。

可推導(dǎo)出干旱歷時、烈度和烈度峰值對應(yīng)的單變量重現(xiàn)期，見文獻［16］。

1.4.2二維干旱變量重現(xiàn)期

二維干旱變量重現(xiàn)期包括聯(lián)合重現(xiàn)期與同現(xiàn)重現(xiàn)期。聯(lián)合重現(xiàn)期Tα與同現(xiàn)重現(xiàn)期T0的計算見文獻［14］。

2　實例分析

2.1單變量的邊緣分布

根據(jù)游程理論對汾河上游的上靜游、汾河水庫、寨上、蘭村4個水文站的月徑流資料進行干旱識別，結(jié)果見表1。

由表1可看出，在調(diào)查年限內(nèi)汾河上游的4個水文站都曾發(fā)生過多場干旱，并且都發(fā)生過歷時超過24個月或烈度超過1億m3或峰值超過300萬m3的嚴(yán)重干旱，越靠近中下游，發(fā)生過的最嚴(yán)重的干旱烈度越大，峰值越高。

采用對數(shù)正態(tài)分布、指數(shù)分布、伽馬分布和廣義帕累托分布分別對汾河上游的4個水文站的干旱歷時、烈度和烈度峰值進行擬合，計算公式見式（1）～式（4）。采用矩法進行參數(shù)估計，利用K-S法進行擬合優(yōu)度檢驗，依據(jù)是K-S統(tǒng)計量的值小于臨界值，結(jié)果見表2。

表1　4個水文站的基本信息及干旱識別結(jié)果

表2　干旱特征變量的4種概率分布的K-S檢驗結(jié)果

由表2可看出，①對于汾河上游4個站的干旱特征變量的擬合，對數(shù)正態(tài)分布全部通過K-S檢驗，擬合效果最佳；指數(shù)分布有50%通過K-S檢驗，廣義帕累托分布有33.3%通過K-S檢驗，擬合效果次之，伽馬分布全部未通過K-S檢驗，擬合效果最差。②上靜游的干旱歷時、烈度和烈度峰值只能由對數(shù)正態(tài)分布擬合。汾河水庫的干旱歷時和烈度可由對數(shù)正態(tài)分布、指數(shù)分布和廣義帕累托分布擬合，其中對數(shù)正態(tài)分布對二者的擬合效果最好；烈度峰值只可由對數(shù)正態(tài)分布擬合。寨上的干旱歷時和烈度可由對數(shù)正態(tài)分布、指數(shù)分布和廣義帕累托分布擬合，其中歷時由對數(shù)正態(tài)分布擬合效果最好；烈度由廣義帕累托分布擬合效果最好；烈度峰值只可由對數(shù)正態(tài)分布擬合。蘭村的干旱歷時和烈度可由對數(shù)正態(tài)分布和指數(shù)分布擬合，指數(shù)分布對二者的擬合效果最好；烈度峰值只可由對數(shù)正態(tài)分布擬合。

利用矩法對各站的干旱特征變量邊緣分布進行參數(shù)估計，結(jié)果見表3。

表3　干旱特征變量邊緣分布及其參數(shù)

2.2干旱特征變量的二維聯(lián)合分布

選取Clayton Copula、GH Copula、AMH Copula 和Nelson No.2 Copula 4種Copula函數(shù)建立干旱特征變量的二維聯(lián)合分布。根據(jù)Kendall秩相關(guān)系數(shù)τ 和Archimedean Copula函數(shù)的參數(shù)θ的解析關(guān)系確定參數(shù)θ，再根據(jù)式（6）進行擬合優(yōu)度評價，依據(jù)是RMSE與AIC值越小，擬合效果越好（見表4）。

由表4可看出，4個站的干旱歷時、烈度與烈度峰值兩兩之間的Kendall秩相關(guān)系數(shù)τ分別介于0.815～0.889、0.583～0.663、0.687～0.729之間，表明干旱特征變量間具有明顯的相關(guān)關(guān)系，可用Copula函數(shù)構(gòu)建其聯(lián)合分布。上靜游站，Clayton Copula對于干旱歷時、烈度和烈度峰值兩兩之間的擬合的效果都最好；汾河水庫站，GH Copula對于干旱歷時和烈度之間擬合效果最好，Clayton Copula對于干旱歷時和烈度峰值、烈度和烈度峰值之間的擬合的效果最好；寨上站，Nelson NO.2 Copula對于干旱歷時和烈度之間擬合效果最好，Clayton Copula對于干旱歷時和烈度峰值之間擬合效果最好，GH Copula對于烈度和烈度峰值之間的擬合的效果最好；蘭村站，GH Copula對于干旱歷時、烈度和烈度峰值兩兩之間的擬合的效果都最好。由于篇幅所限，以寨上站為例，將優(yōu)選的Copula函數(shù)的經(jīng)驗頻率和理論頻率對比如圖2。

表4　二維干旱特征變量的參數(shù)及擬合優(yōu)度評價

由圖2可看出，橫坐標(biāo)表示經(jīng)驗頻率，縱坐標(biāo)表示理論頻率的點據(jù)均勻分布在橫縱坐標(biāo)軸的對角線上及兩側(cè)，表明經(jīng)驗頻率與理論頻率近似相等，說明Copula函數(shù)對于寨上的干旱特征變量的二維分布擬合效果較好。

圖2　寨上站的干旱歷時與干旱烈度、干旱歷時與烈度峰值、干旱烈度與烈度峰值擬合

2.3重現(xiàn)期計算

一組數(shù)據(jù)中最具代表性且最能體現(xiàn)數(shù)據(jù)特征的5個數(shù)是中位數(shù)Me，下四分位數(shù)Q1，上四分位數(shù)Q3，最小值Min和最大值max。因此，選取各站提取的各干旱特征變量序列中的這5個代表性數(shù)據(jù)，分別計算其單變量重現(xiàn)期，再計算二維Copula函數(shù)值，求得對應(yīng)的聯(lián)合重現(xiàn)期（Tα）與同現(xiàn)重現(xiàn)期（T0），結(jié)果見表5。

表5　基于最優(yōu)Copula函數(shù)的水文干旱特征變量的重現(xiàn)期

從表5中的單變量重現(xiàn)期可看出：①4個站的干旱特征變量序列的最小值min、下四分位數(shù)Q1、中位數(shù)me和上四分位數(shù)Q3對應(yīng)的重現(xiàn)期分別為0.74～1.12年、0.9～1.49年、1.37～2.86年和2.18～5.10年。說明這4個站在調(diào)查年份中所發(fā)生的干旱大多數(shù)重現(xiàn)期小于5年。②4個站的干旱特征變量序列的最大值max對應(yīng)的重現(xiàn)期變化范圍很大，其中烈度峰值最大值對應(yīng)的重現(xiàn)期較小，從4.71～17.15 a，干旱歷時與烈度最大值對應(yīng)的重現(xiàn)期較大，分別為 73.73～1045.4 a和 72.95～990.09 a。說明所發(fā)生過的干旱事件中，烈度峰值對應(yīng)重現(xiàn)期比干旱歷時與烈度對應(yīng)的重現(xiàn)期小很多。③從干旱歷時與烈度看，蘭村站在調(diào)查年分內(nèi)干旱最嚴(yán)重，已發(fā)生過超越或接近千年一遇的大旱。

從表5的三種重現(xiàn)期對比可看出：①4個站的二維聯(lián)合重現(xiàn)期小于等于其對應(yīng)的單變量重現(xiàn)期；單變量重現(xiàn)期小于或等于與其對應(yīng)的二維同現(xiàn)重現(xiàn)期。②由于4個站干旱歷時與烈度的相關(guān)性很強，二者的三種重現(xiàn)期都較為接近；干旱歷時與烈度峰值、烈度與烈度峰值的相關(guān)性較弱，它們的三種重現(xiàn)期差距較大。這給今后汾河上游的水文干旱的分析預(yù)測提供了一定依據(jù)。

3　結(jié)論

（1）對數(shù)正態(tài)分布對于汾河上游的干旱特征變量的邊緣分布擬合效果最佳，指數(shù)分布和廣義帕累托分布次之，伽馬分布不適合汾河上游的干旱特征變量的邊緣分布擬合。

（2）汾河上游的4個水文站的單變量重現(xiàn)期介于對應(yīng)的二維聯(lián)合重現(xiàn)期與同現(xiàn)重現(xiàn)期之間。

（3）汾河上游所發(fā)生的干旱大多數(shù)重現(xiàn)期較短，且所發(fā)生的極少數(shù)嚴(yán)重干旱的特點是歷時長、烈度大、峰值并不太高。

（4）汾河上游干旱歷時與烈度的三種重現(xiàn)期都較為接近；干旱歷時與烈度峰值、烈度與烈度峰值的三種重現(xiàn)期差距較大。

［1］SHIAU J T，MODARRES R.Copula-based drought severity-durationfrequency analysis in Iran［J］.Meteorological applications，2009，16 （4）：481-489.

［2］張強，潘學(xué)標(biāo)，馬柱國，等.干旱［M］.北京：氣象出版社，2009：30.

［3］BACCHI B，BECCIU G，KOTTEGODA N T.Bivariate exponential model applied to intensities and durations of extreme rainfall［J］. Journal of Hydrology，1994，55（1/2）：225-236.

［4］Yue S，Ouarda T B M J，Bobee B.A review of bivariate gamma distributions for hydrological application［J］.Journal of Hydrology，2001，246（1）：1-18.

［5］SHIAU J T.Return period of bivariate distributed hydrological events ［J］.Stochastic Environmental Research and Risk Assessment，2003，17（1/2）：42-57.

［6］YUE S.Applying bivariate normal distribution to flood frequency analysis［J］.Water International，1999，124（3）：248-254.

［7］BOX G，COX D.An analysis of transformations［J］.Journal of the Royal Statistical Society，1964，26（2）：211-252.

［8］戴昌軍，梁忠民.多維聯(lián)合分布計算方法及其在水文中的應(yīng)用［J］.水利學(xué)報，2006，37（2）：160-165.

［9］王文圣，丁晶.基于核估計的多變量非參數(shù)隨機模型初步研究［J］.水利學(xué)報，2003，34（2）：9-14.

［10］KIM T W，VALDES J B，YOO C.Nonparametric approach for estimating return periods of droughts in arid regions［J］.Journal of Hydrologic Engineering，2003，8（5）：237-246.

［11］劉和昌，梁忠民，姚軼，等.基于Copula函數(shù)的水文變量條件組合分析［J］.水力發(fā)電，2014，40（5）：13-16.

［12］YEVJEVICH V M.An objective approach to definitions and investigations of continental hydrologic droughts［M］.Colorado State University，1967.

［13］TALLAKSEN L M，MADSEN H，CLAUSEN B.On the definition and modelling of streamflow drought duration and deficit volume［J］. Hydrological Sciences Journal，1997，42（1）：15-33.

［14］宋松柏，蔡煥杰，金菊良，等.Copulas函數(shù)及其在水文中的應(yīng)用［M］.北京：科學(xué)出版社，2012.

［15］NELSON R B.An Introduction to Copulas［M］.New York:Springer，2006.

［16］SHIAU J T，SHEN H W.Recurrence analysis of hydrologic droughts of differing severity［J］.Journal of Water Resources Planning and Management，2001，127（1）：30-40.

（責(zé)任編輯陳萍）

錦屏一級水電站左岸基礎(chǔ)處理工程通過驗收

錦屏一級水電站左岸基礎(chǔ)處理工程順利通過完工驗收。電站左岸基礎(chǔ)處理工程的成功建設(shè)，較好地解決了復(fù)雜地質(zhì)條件制約特高拱壩建設(shè)的關(guān)鍵技術(shù)難題，標(biāo)志著我院復(fù)雜地基處理技術(shù)水平處于世界領(lǐng)先水平。

錦屏一級水電站拱壩高305m，是世界最高拱壩?；A(chǔ)工程地質(zhì)條件極為復(fù)雜，左岸抗力體發(fā)育有斷層、煌斑巖脈、深部裂縫及低波速巖帶、層間擠壓錯動帶等，這些軟弱結(jié)構(gòu)面規(guī)模大、性狀差，對壩體受力狀態(tài)和變形穩(wěn)定產(chǎn)生較大影響。如此復(fù)雜的地質(zhì)條件，極大地增加了拱壩設(shè)計和基礎(chǔ)處理設(shè)計難度，其設(shè)計技術(shù)與施工超出了現(xiàn)行設(shè)計規(guī)范和已有工程經(jīng)驗。

左岸基礎(chǔ)處理工程采取了防滲帷幕灌漿、壩基排水及抗力體排水、左岸混凝土墊座、左岸抗力體固結(jié)灌漿、斷層及煌斑巖脈混凝土網(wǎng)格置換、基礎(chǔ)洞室回填、軟弱巖帶水泥-化學(xué)復(fù)合灌漿等綜合加固處理措施，處理措施的規(guī)模與難度在水電界少有。其經(jīng)驗和成果對促進復(fù)雜地基條件下高拱壩的建設(shè)，具有較高的借鑒價值和指導(dǎo)意義。

（中國電建集團成都勘測設(shè)計研究院有限公司）

Copula-based Analysis of Hydrologic Drought Frequency in the Upper Reaches of Fenhe River

REN Lu，ZHAO Xuehua
（College of Water Resources Science and Engineering，Taiyuan University of Technology，Taiyuan 030024，Shanxi，China）

Monthly runoff data taking from four hydrometric stations in the upper reaches of Fenhe River are applied to identify the droughts by using Run Theory.The drought characteristic variables are fitted by four distribution functions，and the Kolmogorov-Smirnov method is employed to evaluate the goodness-of-fit of univariate marginal distributions.The Copula functions are employed to construct two-dimensional joint distribution models of drought characteristic variables，and then their return periods are determined.The results show that:（a）the lognormal distribution is the candidate marginal distribution function with highest goodness-of-fit in fitting drought characteristic variables in the upper reaches of Fenhe River；（b）the return periods of most of the droughts are less than 5 years，and the univariate return periods and the two-dimensional return periods of drought duration and severity are close to each other；and（c）there are large differences between the univariate return periods and the two-dimensional return periods of drought duration and peak intensity，as well as the drought severity and peak intensity.

hydrologic drought；frequency analysis；Copula function；monthly runoff；upper reaches of Fenhe River

P333.3

A

0559-9342（2016）02-0011-06

2015-03-09

國家自然科學(xué)基金資助項目（40901018）；山西省科技攻關(guān)項目（20140313023-4）；太原理工大學(xué)校青年團隊啟動項目（2013T039）

任璐（1991—），女，山西運城人，碩士研究生，研究方向為水資源系統(tǒng)工程；趙雪花（通訊作者）.