基于Copula函數的區域作物生育階段干旱頻率分析

陳海濤，邱 林，李阿龍

(華北水利水電大學水利學院，鄭州 450011)

我國北方大部分區域，農業干旱問題突出，分析研究區域農業干旱特征，對于加強對干旱的預防和應急管理，促進區域水資源的合理配置具有重要意義[1,2]。

以往針對干旱事件的研究往往采用游程理論從長期實測資料中提取干旱指標(如干旱研究中常用的干旱歷時、干旱烈度、干旱強度等)，然后進行Copula擬合，從而分析研究區域的干旱分布狀況。Copula函數不受各單變量邊際分布類型的影響，可以方便地用于構建變量聯合分布模型，目前在水文事件多元頻率分析中得到了廣泛應用。閆寶偉等[3]應用Copula函數分別構造了最大洪水發生時間及發生時間與其量級之間的聯合分布；馮平等[4]采用Gumbel Copula構造降雨量和入境水量的聯合分布；劉曾美等[5]采用Copula函數構建了區間暴雨和外江洪水位的聯合分布模型；Genest和Favre[6]基于Copula函數建立洪峰洪量聯合分布模型。在干旱研究方面，Copula函數也有大量的應用并取得一定成果，周玉良等、陸桂華等……采用Copula進行區域水文干旱頻率、分布研究[7-11]；陳永勤等[12]、李記等[13]采用Copula函數進行了流域干旱頻率分析；于藝等[14]采用Copula函數進行多變量干旱特征分析；孫可可[15]等結合Copula函數建立抗旱能力指數與干旱頻率的對應關系，得到實際抗旱能力下的農業旱災損失風險曲線。

以往的方法在系統刻畫研究區域干旱特性，分析區域水文干旱、氣象干旱特性方面作用顯著。但由于作物不同生育階段生長需水量存在差異，而不同生育階段缺水對作物產量影響也不相同，而上述方法不能得出作物不同生育階段作物水分虧缺的概率分布狀況，因此，該方法在分析農業干旱方面稍顯不足。本文采用Archimedean Copula函數，結合干旱度指標對作物生育階段降雨量進行分期劃分，構建相鄰生育階段降雨量二維聯合分布函數、分析研究該地區的作物生育階段干旱特征，以期為區域干旱管理提供理論依據。

1 作物干旱指標與生育階段降雨量

作物生育階段水分的數學模型，包含供水時間和數量多少兩方面對作物產量的影響，稱為時間水分生產函數。該類模型將作物的連續生長過程劃分為若干個不同生育階段，認為在作物相同生育階段水分具有同效性，在作物不同生育階段才具有變化。M E Jensen(1968年)構建的Jensen模型給出了作物實際產量與作物相對蒸騰量之間的關系，而作物相對蒸騰量與生育階段相對土壤水分脅迫之間存在某種聯系。結合Jensen模型可定義出作物干旱度指標[16]，見下式。

(1)

WCk=w0,k+Pk+Gk+Ik

(2)

WNk=wk+TEc,k

(3)

w0=1 000H1ηθ1

(4)

(5)

k=(1,2,…,t)

式中：Dr為評價區域干旱度，0≤Dr≤1，Dr=1為作物無收成，水分虧缺程度達到最大；Dr=0為作物豐收，水分不虧缺；λk為作物k階段的敏感系數；WCk為作物第k階段的供水量，mm；w0,k為作物第k階段土壤初始貯水量供水量，mm；Pk為作物第k階段時的有效降雨量，mm；Gk為作物第k階段時的地下水補給量，mm；Ik為作物第k階段時的凈灌水量，mm；WNk為作物第k階段的需水量，mm；TEc,k為作物第k階段的蒸發蒸騰量，mm；wk為作物第k階段正常生長發育所允許的最小土壤貯水量界限值，mm；Hk為作物第k階段時的土壤計劃濕潤層深，m；θ1為作物第1階段時的土壤含水率(以占土壤空隙體積的%計)；t為作物生育階段。

作物干旱度指標反映了區域在特定水文、氣象、土壤及灌溉條件下由于缺水引起的作物減產情況。從以上公式可以看出，對于一個特定地區，作物敏感系數λk、需水量WNk是相對穩定的，作物干旱度指標主要與作物各階段供水量有關。供水量WCk包括了4項內容，其中，對于北方地下水埋深較大的地區，Gk可忽略不計；w0,k主要由前期降雨量決定(w0,k-1主要受更前期的降雨量影響，TEc,k-1相對穩定，Gk-1忽略不計)；階段凈灌水量Ik主要受人類活動影響，非自然因素；Pk為影響供水量的主要指標自然因素。通過上述分析可知，對于北方某一特定區域，在不考慮灌溉的情況下，各生育階段降雨量為該區域干旱度的主要影響因素，如果能夠對某種作物不同生育階段降雨量及不同階段之間的聯系做一系統刻畫，即可全面了解該區域的農業干旱狀況，并針對旱情為不同階段灌溉水量進行預估，從而為整個區域水資源調度提供依據。本文結合干旱度指標中的作物生育階段劃分，以實際調查為依據，統計某作物研究時段內各生育階段降雨量，從而獲取長系列以生育階段為劃分的降雨量序列。

2 基于生育階段降雨量的兩變量聯合分布及檢驗

2.1 Copula函數概述

當選用相鄰兩個生育階段降雨量來共同描述不同生育階段干旱頻率時，則需要度量兩者之間的聯系，即計算它們的聯合概率分布函數。Copula理論已成為當今實現這種相關性分析的有效、通用的方法，其中Archimedean Copula函數因形式簡單而被廣泛應用，常用的Archimedean Copula函數有Gumbel-Hougaard、Frank和ClaytonCopula等十余種，表1列出了這幾種常用Copula函數的表達式，θ為待定參數。

表1 4種常用的二維阿基米德Copula函數

2.2 各生育階段降雨量單變量分布及檢驗

根據sklar定理，單變量概率分布即為多變量聯合分布的邊緣分布函數，合理確定單變量分布對最終確定的Copula函數的精確性有重要意義。大量數據分析計算表明，即使同為生育期降雨量，他們最佳分布函數類型也可能是不一致的，因此，對于不同的生育階段，必須分別分析其分布函數。參考現有干旱水文干旱單變量分布研究成果，選定若干個備選分布函數，針對不同的分布函數對研究區各生育階段降雨數據進行適線，然后根據經驗點據與理論曲線的配合情況，選定各生育周期邊緣分布類型。本文備選的單變量分布函數見表3。

對于某個生育階段的降雨量，本文首先采用極大似然法(MLM)估計各備選分布函數的統計參數，然后進行分布函數的擬合度檢驗與評價，并最終選定該生育階段降雨量的分布類型。具體檢驗評價步驟如下：①給定置信水平α，本文取90%；②利用傳統的X2檢驗在給定置信水平條件下各備選分布函數是否通過檢驗；③若多個分布函數均通過檢驗，則進一步采用的是離差平方和最小準則(OLS)進行評價并最終選取該生育階段降雨量的分布函數。OLS的計算公式為：

(6)

式中：Pi為理論頻率；Pei為經驗頻率；n為樣本容量。

采用該方法可確定各生育階段降雨量的分布函數類型及其統計參數。

2.3 相依性度量

在用Copula函數描述變量間的相關性結構之前，還必須進行單變量相關性分析，以考察各變量之間的相關程度，確保它們是非獨立的。相關系數是衡量隨機變量之間相關性常用的指標，本文Spearman相關系數rn描述相鄰生育階段兩兩之間的相關關系：

(7)

本文采用t-檢驗檢驗兩變量間的其相關性，當計算出的rn大于臨界值rα時，則認為兩個變量具有相關關系。

2.4 基于Copula函數的分布及檢驗

本文采用分步法計算Copula函數的統計參數，即當各生育期降雨量分布確定以后，再對相鄰兩個生育周期降雨量進行Copula擬合。對于某個選定的Copula函數，仍采用MLM法進行估計，求解方程組(8)即可得到Copula函數對應參數估計值。

(8)

式中：θ為參數集；αi為邊緣分布函數參數；Fk為邊緣分布函數；c(·)則為Copula函數的密度函數。

對于多個Copula函數，本文采用AIC信息準則法對其進行優越性評價。AIC信息準則包括兩個部分：Copula函數擬合的偏差與Copula函數參數個數導致的不穩定性，AIC表達為：

(9)

AIC=nlnMSE+2m

(10)

(11)

式中：m為模型參數個數。

判斷依據為，AIC數值越小則說明擬合度越好。

3 二維組合概率和重現期

為了更好的論述不同生育周期聯合概率及其重現期，現對不同生育周期降雨量變量以JYLi(i=0,1,…,n)表示，i代表不同的生育階段，0階段為作物播種前的前期降雨階段，n為某作物劃分的生育階段總數。各生育階段對應的邊緣分布記為ui(i=0,1,…,n)。作為二維分布函數，概率與重現期公式，關鍵取決于所選取變量的不同，本文以第一、第二生育階段聯合分布為例說明各個情況下的聯合分布特性。本文選取如下3種概率/重現期從不同側面對作物生育階段干旱頻率進行分析。

3.1 同現概率與同現重現期

同現概率與同現重現期刻畫了連續兩個生育階段發生干旱狀況的概率分布。

JYL1與JYL2同現不超越概率：

P(JYL1≤jyl1and JYL2≤jyl2)=C(u1,u2)

(12)

JYL1與JYL2不超越某個歷時的同現重現期：

(13)

式中：E(L)and為JYL1與JYL2的同現歷時，即連續兩個JYL1≤jyl1and JYL2≤jyl2年份相距時長的數學期望，年。

3.2 聯合概率與聯合重現期

聯合概率與聯合重現期刻畫了連續兩個生育階段最少一個階段發生干旱的概率分布。

JYL1與JYL2聯合不超越概率：

P(JYL1≤jyl1or JYL2≤jyl2)=u1+u2-C(u1,u2)

(14)

JYL1與JYL2不超越某個歷時的聯合重現期：

(15)

式中：E(L)or為JYL1與JYL2的聯合歷時，即連續兩個JYL1≤jyl1or JYL2≤jyl2年份相距時長的數學期望，單位為年。

3.3 條件概率與重現期

條件概率與條件重現期刻畫了連續兩個生育階段在第一個發生某種干旱狀況的條件下第二生育階段發生干旱的概率分布。

給定JYL1小于等于某定值條件下JYL2的條件概率：

(16)

(17)

式中：E(L)2|1為給定JYL1小于等于某定值條件下JYL2的聯合歷時，即JYL1小于等于某定值條件下，連續兩個JYL2≤jyl2年份相距時長的數學期望，年。

4 應用實例

河南省渠村引黃灌區位于濮陽市西部，屬東亞季風區，灌區內農作物種植面積為22.91 萬hm2，主要作物為小麥、玉米和棉花，現有該灌區1961-2013共53 a逐月降雨資料。該地區玉米生育周期為[17]：第一階段，播種～苗期(06-05-06-25)，降雨量記為jyl1；第二階段，拔節～抽穗(06-25-07-11)降雨量記為jyl2；第三階段，抽穗～乳熟(07-12-08-12) 降雨量記為jyl3；第四階段，乳熟～收獲(08-13-09-02)，降雨量記為jyl4；以播種前15 d(05-20-06-05)作為播種前期以充分考慮前期降雨對土壤含水量的影響，降雨量記為jyl0。本文以濮陽渠村灌區玉米種植為例對Copula函數在農業干旱分析中的應用加以說明。根據上述生育階段的劃分，渠村灌區玉米不同生育階段降雨量統計特征值見表2。

4.1 邊緣分布的確定

針對不同的生育階段，參考水文統計常用的分布函數，利用3.2節所介紹方法進行單變量擬合檢驗。各生育階段參數估計見表3。

表2 各生育階段降雨量統計特征表

表3 各生育階段單變量參數估計成果表

X2檢驗結果為0表明通過檢驗，OLS越小，說明擬合度越好。從表3可以看出，jyl0服從gamma分布，jyl1、jyl2服從Lognormal分布，jyl3、jyl4服從Weibull分布。

4.2 相鄰生育階段降雨量相依性度量

利用3.3節所介紹方法進行變量相關性檢驗，檢驗結果見表4。

表4 相鄰生育階段降雨量相關系系數計算成果表

對相關系數進行顯著性檢驗，對于樣本容量n=53，置信水平α取90%時，臨界相關系數rα為0.2306，從表4可以看出，相鄰生育階段降雨量之間的相關系數均大于臨界相關系數，存在顯著的相關性。

4.3 聯合分布的確定

針對2.1節所描述的幾種Archimedean Copula函數，利用3.4節所述方法對相鄰兩個生育階段降雨量進行擬合并進行測評，其結果見表5。

表5 相鄰生育階段Copula函數參數估計與擬合度評價

表5中，“-”表示統計參數超出了定義域范圍，不參與測評。從表5可看出，相鄰兩個生育階段降雨量均以Gumbel-Houggard Copula為最佳擬合。

4.4 相鄰生育階段降雨量組合概率與重現期分析

根據4.2節確定的Gumbel-Houggard Copula函數參數，利用第3節所述方法計算相鄰生育階段降雨量的同現概率、聯合概率及條件概率，見圖1～圖3。并可計算出其相應的同現重現期、聯合重現期及條件重現期，由于篇幅有限，本文僅列出jyl0、jyl1兩個生育階段對應的成果，見表6～表8，表中“-”為實測數據區，實際應用中若數據落在該區間，則通過插值求得。

圖1 同現概率等值線圖(單位：%)

圖2 聯合概率等值線圖(單位：%)

圖3 條件概率等值線圖(單位：%)

表6 jyl0～ jyl1同現重現期計算成果表 年

表7 jyl0～ jyl1聯合重現期計算成果表 年

表8 jyl0～ jyl1條件重現期計算成果表 年

首先，圖1～圖3、表6～表8以玉米生長過程為例系統反應區域農業干旱狀況，便于對相鄰兩個區域農業干旱狀況進行對比。

其次，可對研究區域某一具體年份的當前農業干旱狀況做出評價，并為下階段灌溉制度提供依據。例如，假定該區域階段降雨量小于階段平均值即會對該階段作物生長產生干旱影響(表2可查出該區域階段降雨量平均值)，則從圖1～圖3可以查出P(JYL0≤25.92 and JYL1≤46.01)約為25%，P(JYL0≤25.92 and JYL1≤46.01)約為65%，P(JYL0≤25.92 and JYL1≤46.01)約為70%；從表6～表8可以查出相應的同現重現期為9年一遇，聯合重現期為2年一遇，條件重現期為4年一遇。這表明該區域玉米生長的前兩個生育階段均發生干旱的概率不是很大(25%)，但這兩個生育階段中至少一個階段發生干旱的概率較大(65%),因此，前兩個生育階段應做抗旱準備；尤其地，根據條件概率，若第一階段發生干旱，則第二階段發生干旱的概率達70%，即，若第一階段干旱，則有必要加強第二階段的抗旱工作，從而為灌區水資源配置提供預估。

5 結 語

本文以濮陽渠村灌區玉米種植為例，對作物生育階段降雨量進行分期劃分。首先，擬合各生育階段降雨量分布函數，然后，選用Archimedean Copula函數構造相鄰生育階段降雨量的聯合分布，計算出相鄰生育周期降雨量的組合概率與組合重現期。該方法分析了相鄰生育階段降雨量之間的多種組合情況，能夠全面地反映區域作物生育階段干旱的特征，可為區域內農業干旱分析和水資源配置提供指導。

[1] 于 藝，宋松柏，馬明衛，等．Archimedean 族 Copulas 函數在多變量干旱特征分析中的應用[J]．水文，2011,31(2):6-10．

[2] 孫可可，陳 進，金菊良，等．實際抗旱能力下的南方農業旱災損失風險曲線計算方法[J]．水利學報，2014,45(7):809-814．

[3] 唐 蘊，王 浩，嚴登華，等．嫩江流域近45年來徑流演變規律研究[J]．地理科學，2009，29(6):864-868．

[4] 金菊良，劉 麗，汪明武，等．基于三角模糊數隨機模擬的地下水環境系統綜合風險評價模型[J]．地理科學，2011，31(2): 143-147．

[5] 閆寶偉，郭生練，陳 璐，等．長江和清江洪水遭遇風險分析[J]．水利學報，2010，41(5) : 553-559．

[6] 馮 平，毛慧慧，王 勇．多變量情況下的水文頻率分析方法及其應用[J].水利學報，2009，40(1) :33-37．

[7] 劉曾美，陳子桑．區間暴雨和外江洪水位遭遇組合的風險[J]．水科學進展，2009，20(5):619-625．

[8] Genest C，Favre A C． Everything you always wanted to know about copula modeling but were afraid to ask [J]．Journal of Hydrologic Engineering，2007，12(4):347-367．

[9] 周玉良，袁瀟晨，金菊良，等．基于Copula的區域水文干旱頻率分析[J]. 地理科學，2011,31(11):1 383-1 388．

[10] 陸桂華，閆桂霞，吳志勇，等．基于copula 函數的區域干旱分析方法[J]．水科學進展，2010，2121(2):188-193.

[11] 徐春曉，袁瀟晨，金菊良，等．基于Copula的區域干旱空間分布特征分析[J]．資源科學，2011,33(12):2 301-2 313．

[12] Shiau J T．Fitting drought duration and severity with two-dimensional Copulas[J]．Water Resources Management，2006，20(5):795-815．

[13] 宋松柏，蔡煥杰，金菊良，等．Copulas函數及其在水文中的應用[M].北京：科學出版社，2012．

[14] 陳永勤，孫 鵬，張 強，等．基于Copula的鄱陽湖流域水文干旱頻率分析[J]．自然災害學報，2013,22(1):75-84．

[15] 李 計，李 毅，賀纏生．基于Copula函數的黑河流域干旱頻率分析[J]．西北農林科技大學學報(自然科學版)，2013，41(1) :213-220．

[16] 陳海濤，黃 鑫，邱 林，等．基于最大熵原理的區域農業干旱度概率分布模型[J]．水利學報，2013,44(2):221-226．

[17] 陳曉楠，段春青，劉昌明，等．基于兩層土壤計算模式的農業干旱風險評估模型[J]．農業工程學報，2009,25(9):51-55．