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

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

0 引 言

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

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

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

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

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

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

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

1.2 研究方法

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

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

式中：ET0為參考作物騰發(fā)量，mm/d；?為飽和水汽壓—溫度曲線上的斜率，kPa/℃；Rn為植物冠層表面凈輻射，MJ/(m2·d)；G為土壤熱通量，MJ/(m2·d)，逐日計算G=0；γ為濕度計常數(shù)，kPa/℃；u2為2 m 高處的風速，m/s；es和ea分別為飽和水汽壓和實際水汽壓，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)合分布的建立具體包括以下幾個步驟：①確定春玉米生育期降雨量和參考作物騰發(fā)量的邊緣分布；②基于生育期降雨量和參考作物騰發(fā)量數(shù)據(jù)，計算Kendall 秩相關(guān)系數(shù)τ；③根據(jù)τ與θ的關(guān)系計算θ；④根據(jù)評價指標選取Copula 函數(shù)，建立聯(lián)合分布；⑤根據(jù)所建分布模型進行豐枯遭遇頻率分析。水文領(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ōu)度評價指標

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

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

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

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

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

以pf＝37.5%，pk＝62.5%作為對富裕和尚志春玉米生育期降雨量與參考作物騰發(fā)量進行豐枯劃分的頻率[23-25]，運用已構(gòu)建的二維Copula 聯(liá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)。

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

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

2 結(jié)果與分析

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

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

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

2.2 Copula函數(shù)評價結(jié)果

表3 單變量邊緣分布參數(shù)估計擬合檢驗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，且對應(yīng)的PKendall小于0.05，說明春玉米生育期降雨量和參考作物騰發(fā)量之間具有顯著的負相關(guān)性，由τ與θ的關(guān)系（表2）可知，只有Frank Copula 函數(shù)可作為春玉米生育期降雨量和參考作物騰發(fā)量的聯(lián)結(jié)函數(shù)（表5），其計算公式為：

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

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

對Frank Copula 函數(shù)進行K-S 檢驗，富裕和海倫的檢驗統(tǒng)計量D分別為0.104、0.092，均小于D0，OLS分別為0.039、0.041，且Frank Copula函數(shù)得出的理論頻率能夠較好地與經(jīng)驗頻率擬合，兩個站點經(jīng)驗頻率與理論頻率的R2均在0.9 以上（圖3），說明由Frank Copula 函數(shù)建立的聯(lián)合分布模型對春玉米生育期降雨量和參考作物騰發(fā)量的聯(lián)合分布情況擬合較好，聯(liá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)等值線，以及兩站點1980-2017年春玉米生育期降雨量和參考作物騰發(fā)量在其中的分布，GX,Y(x,y)表示P枯、ET0豐這兩個事件至少有一個事件發(fā)生的聯(lián)合分布概率；G'X,Y(x,y)表示P枯、ET0豐這兩個事件同時發(fā)生的聯(lián)合分布概率。整體上，富裕和海倫聯(lián)合概率值GX,Y(x,y)、G'X,Y(x,y)變化趨勢是相同的，都是隨著P的增大和ET0的減小而增大。局部看，富裕[圖4（a）]1980-2017年春玉米生育期P-ET0枯豐型中實測值(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枯豐型中實測值(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）計算得出黑龍江省春玉米生育期P枯、ET0豐這兩個事件至少有一個事件發(fā)生的聯(lián)合重現(xiàn)期為2.50年，P枯、ET0豐這兩個事件同時發(fā)生的同現(xiàn)重現(xiàn)期為3.33年。說明在天然來水條件下黑龍江省供水和需水處于不協(xié)調(diào)狀況的概率較高，事件發(fā)生的間隔時間較短，不能滿足作物需水的風險較大。在農(nóng)業(yè)生產(chǎn)中，需要進行灌溉活動、科學調(diào)配水資源。

3 結(jié) 語

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

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

（2）黑龍江省的聯(lián)合概率值GX,Y(x,y)、G'X,Y(x,y)均隨著P的增大和ET0的減小而呈增大趨勢。在天然來水條件下，黑龍江省供水不能滿足作物需水的風險較大，且間隔時間較短，春玉米生育期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）針對水資源供需不協(xié)調(diào)、P-ET0枯豐遭遇風險較大且重現(xiàn)期間隔短的情況，可通過降雨量和ET0的不同遭遇組合情況對其干旱風險進行定性描述和定量評估，從而為制定和調(diào)整灌溉規(guī)劃、科學調(diào)配水資源提供一定參考依據(jù)。