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

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

（北京師范大學(xué)水科學(xué)研究院，北京100875）

0 引 言

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

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

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況與數(shù)據(jù)來(lái)源

黑龍江省位于中國(guó)東北部，地理范圍為東經(jīng)121°11'～135°05'，北緯43°26'～53°33'，東西跨越14 個(gè)經(jīng)度，根據(jù)黑龍江省耕地和水體的空間分布情況，本文選取黑龍江省富裕和海倫2 個(gè)代表性站點(diǎn)進(jìn)行研究（圖1）。黑龍江省屬于寒溫帶與溫帶大陸性季風(fēng)氣候，全省的主要?dú)夂蛱卣鳛椋捍杭镜蜏馗珊?，夏季溫?zé)岫嘤辏锛疽诐吃缢?，冬季寒冷漫長(zhǎng)，無(wú)霜期短，氣候地域性差異大。降水的季風(fēng)性特征明顯，具體表現(xiàn)為夏季降水充沛，冬季干燥少雨，年降水量一般處于400～650 mm 之間[17]?？臻g上，中部山區(qū)降水較多，西部和北部較少。全省多年平均氣溫-5～5℃。年太陽(yáng)輻射總量在4.4×109～5.0×109J/m3之間[18]。全年約5 000～7 000 h 的風(fēng)速不小于3 m/s，風(fēng)能密度為200 W/m3[19]。該地區(qū)的耕作制度為一熟制，主要作物有春玉米、大豆、春小麥等，其中玉米產(chǎn)量占全國(guó)產(chǎn)量的比例為23%，本文以春玉米生育期為代表進(jìn)行研究，即5-9月[20]。

本文的日尺度降水、最高氣溫、最低氣溫、相對(duì)濕度、日照時(shí)數(shù)和平均風(fēng)速數(shù)據(jù)庫(kù)由中國(guó)氣象局提供。數(shù)據(jù)處理主要參照全球觀測(cè)系統(tǒng)世界氣象組織指南（World Meteorological Organization’s guide to the Global Observing System）和中國(guó)氣象局天氣觀測(cè)技術(shù)規(guī)范標(biāo)準(zhǔn)，數(shù)據(jù)下載網(wǎng)站參照http：//cdc.cma.gov.cn。其中日照時(shí)數(shù)通過(guò)公式（參考FAO）轉(zhuǎn)化為日太陽(yáng)輻射量為計(jì)算參考作物騰發(fā)量所用。作物生育期數(shù)據(jù)來(lái)自于中國(guó)氣象局的中國(guó)農(nóng)作物產(chǎn)量資料旬值數(shù)據(jù)集。本文研究所需數(shù)據(jù)的時(shí)間范圍見表1。

表1 研究數(shù)據(jù)描述Tab.1 Description of study data

1.2 研究方法

1.2.1 參考作物騰發(fā)量

1998年FAO推薦的Penman-monteith 公式，所使用的數(shù)據(jù)包括測(cè)站的高程、緯度、風(fēng)速測(cè)量、高度、日最高氣溫、日最低氣溫、日平均氣溫、日平均風(fēng)速、日平均相對(duì)濕度和日照時(shí)數(shù)等。其具體形式如下：

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

1.2.2 二維Copula 函數(shù)

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

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

表2 Copula函數(shù)類型Tab.2 Copula function types

1.2.3 擬合檢驗(yàn)和擬合優(yōu)度評(píng)價(jià)指標(biāo)

本文采用Kolmogorov-Simirnov(K－S)檢驗(yàn)對(duì)Copula 函數(shù)進(jìn)行擬合檢驗(yàn)，采用離差平方和(OLS)最小準(zhǔn)則對(duì)Copula 函數(shù)進(jìn)行擬合優(yōu)度評(píng)價(jià)，其中K－S檢驗(yàn)統(tǒng)計(jì)量D和離差平方和OLS的定義如下：

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

式中：Pi和Pei分別為聯(lián)分布的理論頻率和經(jīng)驗(yàn)頻率。

式中：Num表示聯(lián)合觀測(cè)值對(duì)小于等于值對(duì)(xi,yi)的個(gè)數(shù)；N為聯(lián)合觀測(cè)值的總數(shù)。

1.2.4 聯(lián)合分布概率與重現(xiàn)期

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

豐豐型：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)。

考慮到本文主要研究的是黑龍江省春玉米生育期豐枯遭遇風(fēng)險(xiǎn)研究，所以重點(diǎn)關(guān)注P-ET0枯豐型的兩種聯(lián)合分布概率及其重現(xiàn)期，即：

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

2 結(jié)果與分析

2.1 P、ET0邊緣分布與豐枯劃分結(jié)果

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

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

2.2 Copula函數(shù)評(píng)價(jià)結(jié)果

表3 單變量邊緣分布參數(shù)估計(jì)擬合檢驗(yàn)Tab.3 Univariate marginal distribution parameter estimation fit test

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

富裕和海倫春玉米生育期降雨量和參考作物騰發(fā)量的Kendall 秩相關(guān)系數(shù)τ分別為-0.391、-0.511，且對(duì)應(yīng)的PKendall小于0.05，說(shuō)明春玉米生育期降雨量和參考作物騰發(fā)量之間具有顯著的負(fù)相關(guān)性，由τ與θ的關(guān)系（表2）可知，只有Frank Copula 函數(shù)可作為春玉米生育期降雨量和參考作物騰發(fā)量的聯(lián)結(jié)函數(shù)（表5），其計(jì)算公式為：

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

表5 二維Copula函數(shù)參數(shù)值及擬合優(yōu)度檢驗(yàn)Tab.5 Two-dimensional Copula function parameter value and goodness-of-fit test

對(duì)Frank Copula 函數(shù)進(jìn)行K-S 檢驗(yàn)，富裕和海倫的檢驗(yàn)統(tǒng)計(jì)量D分別為0.104、0.092，均小于D0，OLS分別為0.039、0.041，且Frank Copula函數(shù)得出的理論頻率能夠較好地與經(jīng)驗(yàn)頻率擬合，兩個(gè)站點(diǎn)經(jīng)驗(yàn)頻率與理論頻率的R2均在0.9 以上（圖3），說(shuō)明由Frank Copula 函數(shù)建立的聯(lián)合分布模型對(duì)春玉米生育期降雨量和參考作物騰發(fā)量的聯(lián)合分布情況擬合較好，聯(lián)合分布理論頻率分布模型比較合理。兩個(gè)站點(diǎn)春玉米生育期降雨量與參考作物騰發(fā)聯(lián)合分布模型分別為：

（1）富裕：

（2）海倫：

2.3 P-ET0豐枯遭遇聯(lián)合分布頻率分析

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

綜合富裕和海倫的P-ET0枯豐型聯(lián)合分布概率，由式（7）計(jì)算得出黑龍江省春玉米生育期P枯、ET0豐這兩個(gè)事件至少有一個(gè)事件發(fā)生的聯(lián)合重現(xiàn)期為2.50年，P枯、ET0豐這兩個(gè)事件同時(shí)發(fā)生的同現(xiàn)重現(xiàn)期為3.33年。說(shuō)明在天然來(lái)水條件下黑龍江省供水和需水處于不協(xié)調(diào)狀況的概率較高，事件發(fā)生的間隔時(shí)間較短，不能滿足作物需水的風(fēng)險(xiǎn)較大。在農(nóng)業(yè)生產(chǎn)中，需要進(jìn)行灌溉活動(dòng)、科學(xué)調(diào)配水資源。

3 結(jié) 語(yǔ)

隨著糧食產(chǎn)量的減少、水資源的供給不足以及氣候變化等問(wèn)題的出現(xiàn)，分析作物生育期豐枯遭遇特征，對(duì)指導(dǎo)農(nóng)業(yè)防旱抗旱有重要作用。本文以黑龍江省富裕和海倫1980-2017年的春玉米生育期（5-9月）降雨量和參考作物騰發(fā)量為研究對(duì)象，運(yùn)用Copula函數(shù)構(gòu)建P-ET0豐枯遭遇聯(lián)合分布模型，并對(duì)其進(jìn)行分析得出以下結(jié)論：

（1）黑龍江省春玉米生育期的降雨量與參考作物騰發(fā)量均服從Gamma分布，由Frank Copula函數(shù)建立的聯(lián)合分布模型對(duì)春玉米生育期降雨量和參考作物騰發(fā)量的聯(lián)合分布情況擬合較好，其經(jīng)驗(yàn)頻率與理論頻率的R2大于0.9。

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

（3）針對(duì)水資源供需不協(xié)調(diào)、P-ET0枯豐遭遇風(fēng)險(xiǎn)較大且重現(xiàn)期間隔短的情況，可通過(guò)降雨量和ET0的不同遭遇組合情況對(duì)其干旱風(fēng)險(xiǎn)進(jìn)行定性描述和定量評(píng)估，從而為制定和調(diào)整灌溉規(guī)劃、科學(xué)調(diào)配水資源提供一定參考依據(jù)。