基于Copula函數的黑龍江省作物生育期豐枯遭遇研究

李紅芳，王會肖，趙茹欣，楊雅雪，郭嘉豪

（北京師范大學水科學研究院，北京100875）

0 引 言

隨著糧食產量的減少、水資源的供給不足以及氣候變化等問題的出現，分析研究區域農業干旱特征，對于加強對干旱的預防和應急管理，促進區域水資源的合理配置具有重要意義。降雨量和參考作物蒸騰量是灌溉系統中兩個具有相關關系的重要水文隨機變量，可表征灌區供需水變化關系[1-3]。因此十分有必要通過研究降雨量和參考作物騰發量的二維聯合概率分布模型，來揭示自然降水條件下的灌區水資源短缺風險[4]，近年來，多變量聯合分析成為水文計算領域的一個研究熱點，并證實比單變量分析能更好地描述水文事件的內在規律和分析各個特征屬性之間的相互關系[5]。目前，一些學者已將Copula 函數應用到水文水資源領域[6-9]，Copula 函數能夠獨立于隨機變量的邊緣分布來反映變量之間的相關性，不限定變量的邊緣分布，能夠通過邊緣分布和相關性結構兩部分構建多維聯合分布，并且在轉換過程中不會改變原始隨機變量的信息[10-12]。針對干旱事件的研究往往采用游程理論從長期實測資料中提取干旱指標（如干旱研究中常用的干旱歷時、干旱烈度、干旱強度等），然后進行Copula 擬合，從而分析研究區域的干旱分布狀況[13]。

以往的方法在系統刻畫研究區域干旱特性，分析區域水文干旱、氣象干旱特性方面作用顯著[14]。但由于作物生育期是作物生長需水的關鍵時期，且作物生育期缺水對作物產量影響非常嚴重，而上述方法不能得出作物生育期豐枯遭遇情況，因此，該方法在分析農業干旱方面稍顯不足。本文基于黑龍江省作物生育期降雨量和參考作物騰發量之間的天然相關性，收集研究區1980-2017年的氣象資料，根據Penman-monteith公式求得黑龍江省生育期參考作物騰發量，通過比較簡單易行且已在水文計算中得到成功應用的Copula 方法[15，16]構建作物生育期降雨量和參考作物騰發量聯合分布模型，利用Kolmogorov-Smirnov 檢驗和最小離差平方和（OLS）方法進行擬合檢驗；利用聯合分布模型對其豐枯遭遇頻率進行研究，以期為定量評價灌溉系統的干旱風險提供基礎依據和科學參考，從而促進糧食增產。

1 數據與方法

1.1 研究區概況與數據來源

黑龍江省位于中國東北部，地理范圍為東經121°11'～135°05'，北緯43°26'～53°33'，東西跨越14 個經度，根據黑龍江省耕地和水體的空間分布情況，本文選取黑龍江省富裕和海倫2 個代表性站點進行研究（圖1）。黑龍江省屬于寒溫帶與溫帶大陸性季風氣候，全省的主要氣候特征為：春季低溫干旱，夏季溫熱多雨，秋季易澇早霜，冬季寒冷漫長，無霜期短，氣候地域性差異大。降水的季風性特征明顯，具體表現為夏季降水充沛，冬季干燥少雨，年降水量一般處于400～650 mm 之間[17]。空間上，中部山區降水較多，西部和北部較少。全省多年平均氣溫-5～5℃。年太陽輻射總量在4.4×109～5.0×109J/m3之間[18]。全年約5 000～7 000 h 的風速不小于3 m/s，風能密度為200 W/m3[19]。該地區的耕作制度為一熟制，主要作物有春玉米、大豆、春小麥等，其中玉米產量占全國產量的比例為23%，本文以春玉米生育期為代表進行研究，即5-9月[20]。

本文的日尺度降水、最高氣溫、最低氣溫、相對濕度、日照時數和平均風速數據庫由中國氣象局提供。數據處理主要參照全球觀測系統世界氣象組織指南（World Meteorological Organization’s guide to the Global Observing System）和中國氣象局天氣觀測技術規范標準，數據下載網站參照http：//cdc.cma.gov.cn。其中日照時數通過公式（參考FAO）轉化為日太陽輻射量為計算參考作物騰發量所用。作物生育期數據來自于中國氣象局的中國農作物產量資料旬值數據集。本文研究所需數據的時間范圍見表1。

表1 研究數據描述Tab.1 Description of study data

1.2 研究方法

1.2.1 參考作物騰發量

1998年FAO推薦的Penman-monteith 公式，所使用的數據包括測站的高程、緯度、風速測量、高度、日最高氣溫、日最低氣溫、日平均氣溫、日平均風速、日平均相對濕度和日照時數等。其具體形式如下：

式中：ET0為參考作物騰發量，mm/d；?為飽和水汽壓—溫度曲線上的斜率，kPa/℃；Rn為植物冠層表面凈輻射，MJ/(m2·d)；G為土壤熱通量，MJ/(m2·d)，逐日計算G=0；γ為濕度計常數，kPa/℃；u2為2 m 高處的風速，m/s；es和ea分別為飽和水汽壓和實際水汽壓，kPa；T為日平均氣溫，℃。

1.2.2 二維Copula 函數

Copula 方法的理論基石是Sklar 定理：設X和Y分別為降雨量、參考作物蒸騰量，其聯合分布函數為F(x,y)，邊緣分布函數分別為Fx(x)和Fy(y)，則存在唯一聯合函數Copula 函數Cθ(u,v)，滿足：

Copula 理論聯合分布的建立具體包括以下幾個步驟：①確定春玉米生育期降雨量和參考作物騰發量的邊緣分布；②基于生育期降雨量和參考作物騰發量數據，計算Kendall 秩相關系數τ；③根據τ與θ的關系計算θ；④根據評價指標選取Copula 函數，建立聯合分布；⑤根據所建分布模型進行豐枯遭遇頻率分析。水文領域常用的Copula 函數中的3 種：Frank Copula、Clayton Copula 和Gumbel-Hougarrd Copula 函數[21，22]，本文利用這3 種函數構建二維干旱特征變量的聯合分布（表2）。

表2 Copula函數類型Tab.2 Copula function types

1.2.3 擬合檢驗和擬合優度評價指標

本文采用Kolmogorov-Simirnov(K－S)檢驗對Copula 函數進行擬合檢驗，采用離差平方和(OLS)最小準則對Copula 函數進行擬合優度評價，其中K－S檢驗統計量D和離差平方和OLS的定義如下：

式中：ck為聯合觀測值樣本(xk,yk)的Copula 值；mk為聯合觀測值樣本中滿足條件x≤xk且y≤yk的聯合觀測值的個數。

式中：Pi和Pei分別為聯分布的理論頻率和經驗頻率。

式中：Num表示聯合觀測值對小于等于值對(xi,yi)的個數；N為聯合觀測值的總數。

1.2.4 聯合分布概率與重現期

以pf＝37.5%，pk＝62.5%作為對富裕和尚志春玉米生育期降雨量與參考作物騰發量進行豐枯劃分的頻率[23-25]，運用已構建的二維Copula 聯合分布模型，對二者在年際間的豐枯遭遇情況進行分析，兩個站點的豐枯遭遇理論上可分為9 種情況[15]，其中，X表示降雨量序列，Y表示參考作物騰發量序列：

豐豐型：p1=p(X≥xpf,Y≥ypf)；豐平型：p2=p(X≥xpf,ypk＜Y＜ypf)；豐枯型：p3=p(X≥xpf,Y≤ypk)；

平豐型：p4=p(xpk＜X＜xpf,Y≥ypf)；平平型：p5=p(xpk＜X＜xpf,ypk＜Y＜ypf)；平枯型：p6=p(xpk＜X＜xpf,Y≤ypk)；

枯豐型：p7=p(X≤xpk,Y≥ypf)；枯平型：p8=p(X≤xpk,ypk＜Y＜ypf)；枯枯型：p9=p(X≤xpk,Y≤ypk)。

考慮到本文主要研究的是黑龍江省春玉米生育期豐枯遭遇風險研究，所以重點關注P-ET0枯豐型的兩種聯合分布概率及其重現期，即：

式中：P(X≤x or Y≥y)指X不超過某特定值、Y超過某特定值這兩個事件中至少有一個事件發生的聯合分布概率；P(X≤x,Y≥y)指X不超過某特定值且Y超過某特定值時事件發生的聯合分布概率。Tx,y為聯合重現期，表示兩事件至少有一個發生的重現期；T'x,y為同現重現期，表示兩事件同時發生的重現期。

2 結果與分析

2.1 P、ET0邊緣分布與豐枯劃分結果

圖2為富裕和海倫1980-2017年春玉米生育期降雨量與參考作物騰發量系列，整體上，富裕P-ET0差值波動幅度大于海倫。其中，富裕的參考作物騰發量均大于降雨量，僅1998年和2014年例外；海倫1980-2017年的參考作物騰發量波動比較穩定。

假設春玉米生育期降雨量和生育期參考作物騰發量服從Gamma 分布，Gamma 分布與P-III 型分布同屬于Γ 類分布，是P-III的一種特殊形式，在我國降水研究中適用性強[15]。取KS檢驗的顯著性水平α=0.05,n=38 時，對應的分位點D0=0.215 44，富裕的P和ET0對應的D分別為0.131、0.081，海倫的P和ET0對應的D分別為0.153、0.089，均小于D0，且對應的PK-S值均大于0.05，說明假設成立（表3），且由Gamma 分布可得兩個站點春玉米生育期P和ET0的豐枯劃分數值（表4）。

2.2 Copula函數評價結果

表3 單變量邊緣分布參數估計擬合檢驗Tab.3 Univariate marginal distribution parameter estimation fit test

表4 春玉米生育期降雨量和ET0的豐枯劃分值Tab.4 Precipitation and ET0 during the growth period of spring maize

富裕和海倫春玉米生育期降雨量和參考作物騰發量的Kendall 秩相關系數τ分別為-0.391、-0.511，且對應的PKendall小于0.05，說明春玉米生育期降雨量和參考作物騰發量之間具有顯著的負相關性，由τ與θ的關系（表2）可知，只有Frank Copula 函數可作為春玉米生育期降雨量和參考作物騰發量的聯結函數（表5），其計算公式為：

式中：u、v為降雨量和參考作物騰發量的邊際分布函數。

表5 二維Copula函數參數值及擬合優度檢驗Tab.5 Two-dimensional Copula function parameter value and goodness-of-fit test

對Frank Copula 函數進行K-S 檢驗，富裕和海倫的檢驗統計量D分別為0.104、0.092，均小于D0，OLS分別為0.039、0.041，且Frank Copula函數得出的理論頻率能夠較好地與經驗頻率擬合，兩個站點經驗頻率與理論頻率的R2均在0.9 以上（圖3），說明由Frank Copula 函數建立的聯合分布模型對春玉米生育期降雨量和參考作物騰發量的聯合分布情況擬合較好，聯合分布理論頻率分布模型比較合理。兩個站點春玉米生育期降雨量與參考作物騰發聯合分布模型分別為：

（1）富裕：

（2）海倫：

2.3 P-ET0豐枯遭遇聯合分布頻率分析

在9 種豐枯遭遇組合中，P-ET0枯豐型組合是最不利于灌溉調度的情況，圖4分別繪制了富裕和海倫春玉米生育期PET0聯合分布GX,Y(x,y)、G'X,Y(x,y)等值線，以及兩站點1980-2017年春玉米生育期降雨量和參考作物騰發量在其中的分布，GX,Y(x,y)表示P枯、ET0豐這兩個事件至少有一個事件發生的聯合分布概率；G'X,Y(x,y)表示P枯、ET0豐這兩個事件同時發生的聯合分布概率。整體上，富裕和海倫聯合概率值GX,Y(x,y)、G'X,Y(x,y)變化趨勢是相同的，都是隨著P的增大和ET0的減小而增大。局部看，富裕[圖4（a）]1980-2017年春玉米生育期P-ET0枯豐型中實測值(P≤351.71mmor ET0≥540.01mm)的聯合分布理論概率約為0.4，有13年的P-ET0在該概率范圍內，占總年份的34.21%，其中7年集中在2000-2010年；(P≤351.71mm,ET0≥540.01mm)的聯合分布概率約為0.3，有11年的P-ET0在該概率范圍內，占總年份的28.95%，其中6年集中在2000-2010年。海倫[見圖4（b）]1980-2017年春玉米生育期P-ET0枯豐型中實測值(P≤437.14 mmor ET0≥507.54 mm)的聯合分布理論概率約為0.4，有15年的P-ET0在該概率范圍內，占總年份的39.47%，其中1980-1989年和 2000-2010年分別占有5年；(P≤437.14 mm,ET0≥507.54 mm)的聯合分布概率約為0.3，有9年的P-ET0在該概率范圍內，占總年份的23.68%，其中有3年在1980-1989年，有4年在2000-2007年。

綜合富裕和海倫的P-ET0枯豐型聯合分布概率，由式（7）計算得出黑龍江省春玉米生育期P枯、ET0豐這兩個事件至少有一個事件發生的聯合重現期為2.50年，P枯、ET0豐這兩個事件同時發生的同現重現期為3.33年。說明在天然來水條件下黑龍江省供水和需水處于不協調狀況的概率較高，事件發生的間隔時間較短，不能滿足作物需水的風險較大。在農業生產中，需要進行灌溉活動、科學調配水資源。

3 結 語

隨著糧食產量的減少、水資源的供給不足以及氣候變化等問題的出現，分析作物生育期豐枯遭遇特征，對指導農業防旱抗旱有重要作用。本文以黑龍江省富裕和海倫1980-2017年的春玉米生育期（5-9月）降雨量和參考作物騰發量為研究對象，運用Copula函數構建P-ET0豐枯遭遇聯合分布模型，并對其進行分析得出以下結論：

（1）黑龍江省春玉米生育期的降雨量與參考作物騰發量均服從Gamma分布，由Frank Copula函數建立的聯合分布模型對春玉米生育期降雨量和參考作物騰發量的聯合分布情況擬合較好，其經驗頻率與理論頻率的R2大于0.9。

（2）黑龍江省的聯合概率值GX,Y(x,y)、G'X,Y(x,y)均隨著P的增大和ET0的減小而呈增大趨勢。在天然來水條件下，黑龍江省供水不能滿足作物需水的風險較大，且間隔時間較短，春玉米生育期P枯、ET0豐的聯合概率值GX,Y(x,y)、G'X,Y(x,y)分別為0.4、0.3，聯合重現期和同現重現期分別為年2.50、3.33年。

（3）針對水資源供需不協調、P-ET0枯豐遭遇風險較大且重現期間隔短的情況，可通過降雨量和ET0的不同遭遇組合情況對其干旱風險進行定性描述和定量評估，從而為制定和調整灌溉規劃、科學調配水資源提供一定參考依據。