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

任　璐，趙雪花

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

任璐，趙雪花

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

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

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

0　引言

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

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

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

1　研究方法

1.1干旱事件的定義與識別

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

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

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

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

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

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

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

1.3二維干旱特征變量模型

1.3.1Copula函數定義

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

式中，θ為待定參數。

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

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

式中，n為樣本容量；m為Copula參數估算數；Pc為Copula多元聯合分布計算值；P0為多元聯合分布經驗值。

1.4干旱重現期計算

1.4.1單變量重現期

若干旱事件系列（N年）發生n次，則發生一次干旱事件的年數（重現期）［10］可記為

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

可推導出干旱歷時、烈度和烈度峰值對應的單變量重現期，見文獻［16］。

1.4.2二維干旱變量重現期

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

2　實例分析

2.1單變量的邊緣分布

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

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

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

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

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

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

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

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

2.2干旱特征變量的二維聯合分布

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

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

表4　二維干旱特征變量的參數及擬合優度評價

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

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

2.3重現期計算

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

表5　基于最優Copula函數的水文干旱特征變量的重現期

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

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

3　結論

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

（2）汾河上游的4個水文站的單變量重現期介于對應的二維聯合重現期與同現重現期之間。

（3）汾河上游所發生的干旱大多數重現期較短，且所發生的極少數嚴重干旱的特點是歷時長、烈度大、峰值并不太高。

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

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

［2］張強，潘學標，馬柱國，等.干旱［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］戴昌軍，梁忠民.多維聯合分布計算方法及其在水文中的應用［J］.水利學報，2006，37（2）：160-165.

［9］王文圣，丁晶.基于核估計的多變量非參數隨機模型初步研究［J］.水利學報，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函數的水文變量條件組合分析［J］.水力發電，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函數及其在水文中的應用［M］.北京：科學出版社，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.

（責任編輯陳萍）

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

錦屏一級水電站左岸基礎處理工程順利通過完工驗收。電站左岸基礎處理工程的成功建設，較好地解決了復雜地質條件制約特高拱壩建設的關鍵技術難題，標志著我院復雜地基處理技術水平處于世界領先水平。

錦屏一級水電站拱壩高305m，是世界最高拱壩。基礎工程地質條件極為復雜，左岸抗力體發育有斷層、煌斑巖脈、深部裂縫及低波速巖帶、層間擠壓錯動帶等，這些軟弱結構面規模大、性狀差，對壩體受力狀態和變形穩定產生較大影響。如此復雜的地質條件，極大地增加了拱壩設計和基礎處理設計難度，其設計技術與施工超出了現行設計規范和已有工程經驗。

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

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

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

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

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