顧世祥，趙 眾，陳 晶，陳金明,3，張劉東

基于高維Copula函數(shù)的逐日潛在蒸散量及氣象干旱預(yù)測(cè)

顧世祥1,2，趙 眾1，陳 晶1,2，陳金明1,2,3，張劉東1

（1. 云南農(nóng)業(yè)大學(xué)城鄉(xiāng)水安全與節(jié)水減排高校重點(diǎn)實(shí)驗(yàn)室，昆明 650201；2. 云南省水利水電勘測(cè)設(shè)計(jì)研究院，昆明 650021；3. 武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072）

嘗試引入高維Copula函數(shù)對(duì)影響參考作物蒸散量ET0的氣象因素進(jìn)行聯(lián)合分布構(gòu)建，揭示不同變量間的相關(guān)結(jié)構(gòu)，建立多元?dú)庀笠蛩貙?duì)ET0的聯(lián)合分布模型，對(duì)逐日ET0及短期干旱等級(jí)進(jìn)行預(yù)測(cè)，并將枯季1—4月份的多維Copula聯(lián)合分布預(yù)測(cè)模型的系統(tǒng)性偏差構(gòu)造成修正函數(shù)，代回ET0預(yù)報(bào)模型以改善預(yù)報(bào)效果，利用洱海流域內(nèi)大理站1954—2018年逐日氣象觀測(cè)數(shù)據(jù)，以FAO Penman-Monteith方程為標(biāo)準(zhǔn)值對(duì)比分析。結(jié)果表明：1）平均氣溫（）和最高氣溫（max）2個(gè)氣象因子組合時(shí)，二維Normal Copula模型對(duì)逐日ET0預(yù)測(cè)的精度最高，疊加上修正函數(shù)項(xiàng)之后，相對(duì)誤差小于10%、15%、20%、25%的樣本比例分別提高到71.6%、84.4%、91.4%、96.5%，全年符合指數(shù)IA變化范圍為0.98～0.99，平均偏差ME為0.17～0.30，均方根誤差RMSE為0.54～0.64，Nash-Sutcliffe效率系數(shù)為0.90～0.98；2）將逐日ET0預(yù)測(cè)方法應(yīng)用于逐日氣象干旱預(yù)測(cè)評(píng)估（以逐日SPEI指數(shù)為例），逐日SPEI指數(shù)預(yù)測(cè)值與標(biāo)準(zhǔn)值的相關(guān)系數(shù)為0.95～0.99，平均偏差ME為-0.10～0.35，均方根誤差RMSE為0.20～0.30，符合指數(shù)IA為0.97～0.98，Nash-Sutcliffe效率系數(shù)NSE為0.91～0.97，在降水量多的季節(jié)，Copula函數(shù)模型預(yù)測(cè)ET0的精度更高一些，且逐日SPEI預(yù)測(cè)的誤差參數(shù)都優(yōu)于逐日ET0的預(yù)測(cè)結(jié)果。

干旱；蒸散量；SPEI指數(shù)；高維Copula函數(shù)；預(yù)測(cè)；洱海流域

0 引 言

潛在蒸散量ET0是影響作物需水氣象因素分項(xiàng)的綜合反映[1]，是氣象因子之間復(fù)雜的非線性關(guān)系描述。ET0除了傳統(tǒng)意義上作為灌區(qū)和田間尺度的灌溉用水調(diào)度管理、農(nóng)業(yè)水土資源平衡及水資源優(yōu)化配置的輸入項(xiàng)外，還是揭示全球及區(qū)域農(nóng)業(yè)氣候變化、干旱災(zāi)害、生態(tài)環(huán)境監(jiān)測(cè)等的重要指標(biāo)[2-4]。逐日ET0計(jì)算模型主要包括Penman Monteith方程[1,4]、Priestley-Taylor模型[5]、Hargreaves公式[6]、Irmark-Allen法[7]等，或者是基于地面能量平衡的區(qū)域ET0估算、地表蒸散發(fā)遙感反演[8-9]。楊永剛等[10]利用Arc GIS空間插值、敏感性分析和貢獻(xiàn)率等對(duì)中國(guó)糧食主產(chǎn)區(qū)265個(gè)站點(diǎn)1961—2013年53a氣象數(shù)據(jù)及ET0進(jìn)行分析。Traore等[11]采用氣象因子基于人工智能網(wǎng)絡(luò)預(yù)測(cè)ET0，發(fā)現(xiàn)最高氣溫是很重要的因素，模型預(yù)測(cè)的準(zhǔn)確性取決于太陽(yáng)凈輻射（Rs）的精度和準(zhǔn)確的天氣預(yù)報(bào)信息。Yoo等[12]在中東、北非、西歐等地區(qū)現(xiàn)有11~21個(gè)區(qū)域模型基礎(chǔ)上，構(gòu)建了以輻射因子作為修正參數(shù)的大區(qū)域ET0計(jì)算模型。王振龍等[13]采用三基點(diǎn)溫度分析了逐日蒸散的動(dòng)態(tài)變化，發(fā)現(xiàn)通過(guò)溫度模擬冬小麥和夏玉米作物系數(shù)變化的擬合度較高。Islam等[14]對(duì)于孟加拉國(guó)的水稻蒸散量進(jìn)行研究，發(fā)現(xiàn)水稻蒸散量增加的主要原因是最高氣溫升高的速率大于太陽(yáng)凈輻射（Rs）升高率。以上成果均是基于單變量或雙變量對(duì)ET0預(yù)測(cè)分析，鮮有全面考慮平均氣溫（）、最低氣溫（min）、最高氣溫（max）、日照時(shí)數(shù)（）、風(fēng)速（）、相對(duì)濕度（RH）和降雨量開(kāi)展對(duì)ET0的聯(lián)合分布研究。高維Copula函數(shù)能有效地描述氣象因素特征變量間的相依性，構(gòu)造任意邊緣分布的聯(lián)合分布函數(shù)，較好地刻畫(huà)變量間的相關(guān)結(jié)構(gòu)，被越來(lái)越多的應(yīng)用到多變量事件分析中[15-17]，如嘗試用概率密度匹配估算潛在蒸散量[18]。Copula函數(shù)理論方法在水文水資源學(xué)領(lǐng)域的降水過(guò)程預(yù)報(bào)、洪水頻率分析、干旱評(píng)估、多站徑流模擬等的應(yīng)用已較廣泛[19-20]。ET0變化是多個(gè)氣象因素共同作用的結(jié)果，可探索通過(guò)maxmin、、、RH等構(gòu)建多元影響因素聯(lián)合分布，研究其用于ET0實(shí)時(shí)預(yù)報(bào)的可行性。地表濕潤(rùn)指數(shù)（Surface Wetness Index，SWI）、標(biāo)準(zhǔn)化降水指數(shù)（Standardized Precipitation Index，SPI）、標(biāo)準(zhǔn)化降水蒸散指數(shù)（Standardized Precipitation Evapotranspiration Index，SPEI）是目前常用的干旱評(píng)估指數(shù)[21-23]，在國(guó)家、區(qū)域、流域等不同空間尺度下的干旱強(qiáng)度分布格局、長(zhǎng)期變化趨勢(shì)及動(dòng)態(tài)監(jiān)測(cè)評(píng)估等已有廣泛應(yīng)用[24-26]。

洱海位于云南省大理白族自治州，屬瀾滄江—湄公河水系，流域面積2 785 km2，年均氣溫15.1 ℃，年降水1 057 mm，水資源總量10.7億m3。流域干濕季明顯，濃郁的蒼山洱海自然風(fēng)情和南詔古國(guó)文明遺跡，具有“東方日內(nèi)瓦”美譽(yù)[27]。全球氣候變化及低緯度高原的區(qū)域響應(yīng)，洱海地區(qū)遭遇了2010—2015年的連續(xù)干旱災(zāi)害，農(nóng)業(yè)用水及農(nóng)田面源加大、湖水位下降、局部區(qū)域藍(lán)藻爆發(fā)，洱海流域水生態(tài)保護(hù)治理已成為全國(guó)關(guān)注焦點(diǎn)之一。利用短期氣象資料研究洱海地區(qū)的逐日ET0及SPEI氣象干旱指數(shù)預(yù)測(cè)評(píng)估方法模型，對(duì)流域內(nèi)灌溉用水的精細(xì)化調(diào)度管理、農(nóng)業(yè)節(jié)水減排及洱海入湖面源控制等具有重要意義。

1 材料與方法

構(gòu)建maxminRH等6種氣象因素的邊緣分布函數(shù)，建立六維Copula模型，并對(duì)其進(jìn)行優(yōu)選。篩選出合適的Copula模型之后，在對(duì)其分別構(gòu)建二維和三維的Copula模型。將氣象因素代入模型之后，計(jì)算出不同氣象因素之間的聯(lián)合分布，記為預(yù)測(cè)ET0的分布概率，再代回其邊緣分布函數(shù)中計(jì)算出ET0即為預(yù)測(cè)值。為評(píng)估本文方法預(yù)測(cè)ET0的精度和效用，逐日ET0標(biāo)準(zhǔn)值采用聯(lián)合國(guó)糧農(nóng)組織（Food and Agriculture Organization of the United Nations，F(xiàn)AO）推薦的標(biāo)準(zhǔn)ET0估算方法（FAO-56 PM），即Penman-Monteith方程計(jì)算[1,4,10]。將上述方法運(yùn)用于日SPEI的計(jì)算模型中，可進(jìn)行短期干旱預(yù)報(bào)。分別將實(shí)際ET0和預(yù)測(cè)ET0運(yùn)用于日SPEI計(jì)算，并以實(shí)際ET0進(jìn)行SWI干旱指數(shù)計(jì)算與之對(duì)比，驗(yàn)證其準(zhǔn)確性。

Copula函數(shù)能把多維隨機(jī)變量1,…,X的聯(lián)合分布(1,…,x)與它們各自的邊緣分布1(1),…,F(x)相連接。本研究主要包括2種：正態(tài)Copula函數(shù)[28]和t-Copula函數(shù)[29]。在求Copula函數(shù)的參數(shù)時(shí)，定義不同的相關(guān)陣來(lái)反映數(shù)據(jù)資料的特點(diǎn)，對(duì)于maxminRH 6種氣象因素，選擇以下3種相關(guān)結(jié)構(gòu)：可交換相關(guān)結(jié)構(gòu)、Toeplitz結(jié)構(gòu)和無(wú)結(jié)構(gòu)相關(guān)。3種相關(guān)結(jié)構(gòu)形式分別為

采用極大似然數(shù)函數(shù)估計(jì)法對(duì)上述Copula函數(shù)進(jìn)行參數(shù)估計(jì)，運(yùn)用K-S檢驗(yàn)法進(jìn)行擬合優(yōu)度檢驗(yàn)，運(yùn)用AIC準(zhǔn)則和BIC準(zhǔn)則進(jìn)行擬合優(yōu)度評(píng)價(jià)，從而比選出最優(yōu)Copula函數(shù)[30]。表1為Copula模型的擬合優(yōu)度評(píng)價(jià)結(jié)果，模型選擇為無(wú)結(jié)構(gòu)相關(guān)的Normal Copula模型。

表1 Copula模型選擇結(jié)果

由于ET0是多個(gè)氣象因素相互作用的結(jié)果，Copula函數(shù)的計(jì)算為氣象因素的聯(lián)合分布，即Copula函數(shù)計(jì)算結(jié)果的聯(lián)合分布與ET0的概率分布結(jié)果相等同。則可將Copula函數(shù)的結(jié)果記為ET0的預(yù)測(cè)概率，將ET0的預(yù)測(cè)概率代入ET0的概率分布函數(shù)，還原得ET0的預(yù)測(cè)值。

為計(jì)算出ET0的預(yù)測(cè)值，運(yùn)用正態(tài)分布、Gamma分布、Lognormal分布、Weibull分布等分別計(jì)算ET0的分布，選擇出最優(yōu)分布，將ET0預(yù)測(cè)概率值代入分布函數(shù)，計(jì)算出ET0預(yù)測(cè)值。如表2所示，根據(jù)ET0的各個(gè)分布函數(shù)之?dāng)M合優(yōu)度檢驗(yàn)結(jié)果，基于最小值準(zhǔn)則，ET0的分布函數(shù)選擇為Gamma函數(shù)。

賈艷青等[26]改進(jìn)了SPEI，改進(jìn)的逐日SPEI計(jì)算過(guò)程與月SPEI類似。計(jì)算SPEI時(shí)應(yīng)先給定時(shí)間尺度，然后計(jì)算P-PE的累積序列，再采用廣義邏輯分布對(duì)累積序列進(jìn)行擬合，最后對(duì)累積概率密度進(jìn)行標(biāo)準(zhǔn)化處理得到日SPEI。干旱劃分等級(jí)參照國(guó)際上通用的SPEI指數(shù)干旱等級(jí)劃分標(biāo)準(zhǔn)[31]。SWI定義為年降水量與年潛在蒸散量的比值[25]，其干旱指數(shù)采用Ma等[32]劃分標(biāo)準(zhǔn)。

表2 逐日ET0的分布函數(shù)擬合優(yōu)度檢驗(yàn)

Note:K-S, Kolmogorov-Smirnov;AIC, Akaike Information Criterion; BIC, Bayesian Information Criterion.

使用4個(gè)常用的統(tǒng)計(jì)指數(shù)來(lái)評(píng)估模擬的效果，檢驗(yàn)Copula模型預(yù)報(bào)逐日ET0的精度，分別為平均偏差（Mean Deviation，ME）、均方根誤差（Root Mean Square Error，RMSE）、符合指數(shù)（Coincidence Index，IA）、Nash-Sutcliffe效率系數(shù)（Nash-Sutcliffe Efficiency Coefficient，NSE）[33]。符合指數(shù)（IA）在0和1之間的范圍內(nèi)，IA越大模擬效果越好。Nash-Sutcliffe效率系數(shù)（NSE）變化范圍從-∞到1，值越接近1說(shuō)明模擬值和實(shí)際值越接近。

2 數(shù)據(jù)資料來(lái)源

從洱海流域內(nèi)大理氣象站獲取1954—2018年逐日平均氣溫、最低氣溫、最高氣溫、日照時(shí)數(shù)、風(fēng)速、相對(duì)濕度和降雨量等氣象觀測(cè)資料；大理、洱源、賓川等洱海相關(guān)市縣的自然地理、社會(huì)經(jīng)濟(jì)、農(nóng)業(yè)水利等綜合統(tǒng)計(jì)年報(bào)（年鑒）和文獻(xiàn)資料獲得相關(guān)現(xiàn)狀信息。

3 結(jié)果與分析

3.1 建立Copula模型預(yù)測(cè)逐日ET0

選用有代表性的正態(tài)分布、Gamma分布、Lognormal分布、Weibull分布分別構(gòu)建maxminRH的邊緣分布函數(shù)，并從中選擇擬合效果最好的分布函數(shù)。

如表 3，邊緣分布函數(shù)擬合優(yōu)度檢驗(yàn)，選擇Weibull分布，max選擇Weibull分布，min選擇正態(tài)分布，選擇Weibull分布，選擇Gamma分布，RH選擇正態(tài)分布。

表3 邊緣分布函數(shù)的擬合優(yōu)度檢驗(yàn)

注：、min、max、、、RH分別表示平均氣溫、最低氣溫、最高氣溫、日照時(shí)數(shù)、風(fēng)速和相對(duì)濕度。

Note:,min,max,,and RH represent average temperature, minimum temperature, maximum temperature, sunshine duration, wind speed and relative humidity,respectively.

分別建立-min-max---max-minmax-max-max-min-min-min-等12種組合形式的二維Normal Copula模型。分別計(jì)算出Normal Copula聯(lián)合分布函數(shù)數(shù)值之后，將其結(jié)果標(biāo)記為預(yù)測(cè)ET0的概率值，代入ET0的邊緣分布函數(shù)Gamma函數(shù)中，求得ET0預(yù)測(cè)概率值。將逐日氣象數(shù)據(jù)分別導(dǎo)入上述對(duì)應(yīng)模型中，計(jì)算出對(duì)應(yīng)日期的ET0預(yù)測(cè)概率值。

3.2 預(yù)測(cè)結(jié)果與分析

為方便比較，選擇年降水頻率為特豐水年（=5%）1993年，豐水年（=25%）1987年，平水年（=50%）1996年，中等干旱年（=75%）1986年，特枯水年（=95%）1960年。將5個(gè)典型年的預(yù)測(cè)結(jié)果與其逐日ET0標(biāo)準(zhǔn)值進(jìn)行對(duì)比，并以此驗(yàn)證模型的準(zhǔn)確性。每個(gè)模型對(duì)比有1 827個(gè)樣本，在-max兩個(gè)氣象因素組合時(shí)預(yù)測(cè)精度最好，相對(duì)誤差（Relative Error，ERR）小于10%、15%、20%、25%、30%的樣本占比分別為59.7%、73.6%、86.0%、91.9%、96.4%。

分析發(fā)現(xiàn)，逐日ET0對(duì)于-max的組合相依性最好，其次為-組合。且只要有參與的Normal Copula模型，其計(jì)算結(jié)果誤差均較大。顯然，ET0對(duì)于-max的相依性最好，而其對(duì)于的相依性最差。

為繼續(xù)探索逐日ET0與其他氣象因素的相依性關(guān)系，再以-max分別與min和組合建立三維Normal Copula模型，為-max-min-max--max-和-max-等4組。逐日ET0預(yù)測(cè)值與標(biāo)準(zhǔn)值的相對(duì)誤差在-max-min組合時(shí)最小，在-max-組合時(shí)誤差最大，且在使用三維Normal Copula模型時(shí)-max-min3個(gè)氣象因子組合最優(yōu)，但仍不如-max組合的預(yù)測(cè)精度高。從各個(gè)氣象因子組合的二維、三維Copula模型進(jìn)行逐日ET0預(yù)測(cè)值與標(biāo)準(zhǔn)值比較的誤差曲線可以看出，1—4月份時(shí)都出現(xiàn)預(yù)測(cè)值整體性較標(biāo)準(zhǔn)值偏小且無(wú)法消除的現(xiàn)象，這可能與Copula函數(shù)本身對(duì)降水量敏感性高有關(guān)。為此，又另選出一組典型年：特豐水年（=7%）2000年、豐水年（=37%）2016年、平水年（=50%）2001年、中等干旱年（=77%）2003年、特枯水年（=97%）2012年，亦采用-max組合的Copula模型進(jìn)行逐日ET0預(yù)測(cè)，與前一組典型年對(duì)應(yīng)各個(gè)水文頻率年景的逐日ET0預(yù)測(cè)誤差曲線相比較，發(fā)現(xiàn)1—4月的變化趨勢(shì)基本一致，對(duì)應(yīng)到各旬統(tǒng)計(jì)其平均值，得到新的系列ET0,i=()，為旬序數(shù)，大致呈“M”型雙峰變化，峰值分別出現(xiàn)在1月中旬和3月中旬，且第一個(gè)峰值較小、持續(xù)時(shí)間短，第二個(gè)峰值持續(xù)了1個(gè)月左右，如圖1所示。將ET0,i作為Copula聯(lián)合分布模型預(yù)測(cè)逐日ET0的修正函數(shù)項(xiàng)，代回預(yù)報(bào)模型中以改善預(yù)測(cè)效果。將改進(jìn)后的各典型年逐日ET0預(yù)測(cè)的標(biāo)準(zhǔn)值與預(yù)測(cè)值比較如圖 2。由圖可知，從=5%特豐水年到=95%特枯水年，逐日ET0預(yù)測(cè)值與標(biāo)準(zhǔn)值的相關(guān)系數(shù)為0.903～0.946，相對(duì)誤差小于10%、15%、20%、25%的樣本數(shù)分別提高到71.6%、84.4%、91.4%、96.5%。

誤差分析結(jié)果見(jiàn)表4，修正后ME值全年變化范圍在0.17～0.30之間，1―6月在0.06～0.14之間變化，7—12月在0.20～0.48之間變化；RMSE值全年變化范圍在0.54～0.64之間，1―6月變化范圍在0.50～0.66之間，7―12月份變化范圍在0.55～0.74之間；全年平均IA變化范圍在0.98～0.99之間，1―6月份IA變化范圍在0.99～0.98之間，7―12月IA變化范圍在0.95～0.99之間；NSE值全年變化范圍在0.90～0.98之間，1―6月變化范圍在0.96～0.98之間，7―12月變化范圍在0.73～0.93之間，其預(yù)測(cè)值與標(biāo)準(zhǔn)值很接近。

圖1 修正函數(shù)變化過(guò)程示意圖

3.3 短期氣象干旱等級(jí)預(yù)測(cè)

不論是SPEI還是改進(jìn)的逐日SPEI，計(jì)算潛在蒸散量均采用是Penman Monteith方程[22,25-26]，該方程的計(jì)算結(jié)果理論最完備，但所需氣象因素較多，在一些氣象因素缺失、監(jiān)測(cè)不全及干旱風(fēng)險(xiǎn)管理實(shí)踐中，很難對(duì)其進(jìn)行準(zhǔn)確計(jì)算。本次引入Copula函數(shù)，對(duì)其氣象因素與ET0的相依性進(jìn)行分析，在缺失氣象因素時(shí)，只需要提供日最高氣溫和平均氣溫，即可對(duì)逐日ET0進(jìn)行預(yù)測(cè)計(jì)算。將Copula函數(shù)模型與逐日SPEI干旱指數(shù)相結(jié)合，既能保證干旱指數(shù)對(duì)氣象因素的充分考慮，又避免氣象因素缺失導(dǎo)致的計(jì)算失真問(wèn)題。本文選用-max二維Normal Copula模型對(duì)ET0的預(yù)測(cè)結(jié)果與逐日SPEI干旱指數(shù)結(jié)合進(jìn)行干旱預(yù)測(cè)。

表4 改進(jìn)T-Tmax二維Normal Copula模型逐日ET0預(yù)測(cè)的誤差分析

圖2 修正之后大理不同水文頻率下的ET0標(biāo)準(zhǔn)值與預(yù)測(cè)值

采用K-S檢驗(yàn)Log-logistic分布，檢驗(yàn)了不同水文頻率年份在30日尺度的逐日SPEI指數(shù)，水文頻率由豐到枯為=5%、25%、50%、75%、95%，不同典型年對(duì)應(yīng)的值為1.889、2.569、1.837、3.007、2.167，均滿足Log-logistic分布。圖3是采用預(yù)測(cè)ET0與降水量計(jì)算逐日SPEI干旱指數(shù)得出大理站1960年、1986年、1987年、1993年和1997年的逐日干旱預(yù)測(cè)情況，并與采用標(biāo)準(zhǔn)ET0與降水量計(jì)算的SPEI實(shí)際結(jié)果相比較。表5為加入ET0預(yù)測(cè)修正函數(shù)項(xiàng)之后的逐日SPEI預(yù)測(cè)誤差分析，逐日SPEI指數(shù)預(yù)測(cè)值與實(shí)際值相關(guān)系數(shù)為0.95～0.99，ME為-0.10～0.35，RMSE為0.20～0.30，IA為0.97～0.98，NSE為0.91～0.97。顯然，采用預(yù)測(cè)SPEI與實(shí)際SPEI對(duì)逐日干旱等級(jí)預(yù)測(cè)的趨勢(shì)相同，且逐日SPEI預(yù)測(cè)的誤差參數(shù)都優(yōu)于逐日ET0的預(yù)測(cè)，原因可能是在ET0計(jì)算模型中缺少的降水因子，被逐日SPEI預(yù)測(cè)模型補(bǔ)充進(jìn)來(lái)，綜合反映氣溫、降水等對(duì)干旱指數(shù)的作用，從而提高了逐日SPEI的預(yù)測(cè)精度。

圖3 修正后不同水文年下大理站逐日SPEI實(shí)際值與預(yù)測(cè)值

將逐日ET0預(yù)測(cè)方法應(yīng)用于逐日氣象干旱（日SPEI指數(shù)）預(yù)測(cè)評(píng)估可得，由豐水年至特枯水年，無(wú)旱及輕旱的天數(shù)比例由81.32%減至46%，中旱由10.68%增至27.87%，重旱及特旱由8.2%增至26.2%。在5個(gè)典型年中，采用逐日ET0標(biāo)準(zhǔn)值進(jìn)行干旱預(yù)測(cè)時(shí)，出現(xiàn)重旱及特旱的天數(shù)為302 d，頻次為16.53%，而采用修正預(yù)測(cè)逐日ET0進(jìn)行干旱預(yù)測(cè)時(shí)，出現(xiàn)重旱或特旱的天數(shù)為324 d，頻次為17.73%，相對(duì)偏差為1.2%。按季節(jié)分，洱海地區(qū)上半年1―6月無(wú)旱及輕旱頻次為36.51%，中旱頻次為30.37%，重旱及特旱頻次為33.11%，預(yù)測(cè)與實(shí)際偏差為2.1%；下半年7―12月無(wú)旱及輕旱頻次為89.73%，中旱頻次為9.07%，重旱及特旱頻次為1.2%，預(yù)測(cè)與實(shí)際偏差為1%。結(jié)合前文二維Normal Copula模型前半年1―6月份預(yù)測(cè)ET0與實(shí)際ET0的誤差比后半年7―12月份的大，可知降水量越高，該模型預(yù)測(cè)ET0及SPEI的精度會(huì)更高一些。這與其他方法一般在枯水（或干旱）年景時(shí)模擬逐日ET0預(yù)報(bào)效果更佳的規(guī)律相左。

表5 修正后逐日SPEI預(yù)測(cè)誤差分析

此外，采用SWI、SPEI等干旱指數(shù)，基于大理站1954—2018年逐月降水及潛在蒸散量分析了其65 a的干旱指數(shù)。根據(jù)K-S檢驗(yàn)結(jié)果，得到在1、3、12個(gè)月尺度下的指數(shù)分別為1.964、2.683、3.126，均滿足Log-logistic分布，這與王林等[34]通過(guò)干旱指數(shù)監(jiān)測(cè)中國(guó)干旱適用性的結(jié)果是一致的。近65年來(lái)洱海流域大理站全年無(wú)旱及輕旱的月份占比為44.4%，中旱月份占33.6%，重旱及特旱月份占22%。其中上半年1―6月份無(wú)旱及輕旱月份占比為8.9%，中旱月份占47.1%，重旱及特旱月份占44%；下半年7―12月份無(wú)旱及輕旱月份占80%，中旱月份占20%，重旱及特旱月份為0。即前半年為枯季，重旱及特旱主要出現(xiàn)在4―6月份，下半年雨季主要為輕旱及無(wú)旱。圖4和圖5為大理站近65年來(lái)各月的干旱等級(jí)及各干旱等級(jí)天數(shù)比例情況。在季節(jié)性干旱的氣候環(huán)境下，從枯水年至豐水年，重旱及特旱月份由5個(gè)月逐漸減少為0，輕旱及無(wú)旱月份由0逐漸增加至7個(gè)月，其結(jié)果與此前采用干濕指數(shù)Ia的研究結(jié)論一致[35]。圖6為其中1960、1986、1987、1993年和1997年等典型年逐日統(tǒng)計(jì)得到各月的SWI和SPEI變化情況，因2個(gè)干旱指標(biāo)的數(shù)量關(guān)系難以統(tǒng)一，以數(shù)據(jù)點(diǎn)落在相同干旱等級(jí)區(qū)域視為結(jié)果相同，則2種干旱指標(biāo)在各個(gè)典型年景的年內(nèi)變化趨勢(shì)基本一致，僅個(gè)別月份出現(xiàn)偏差，其余絕大部分時(shí)間2個(gè)指標(biāo)的評(píng)價(jià)結(jié)果完全相同，表明本文的方法合理可靠。

注：標(biāo)準(zhǔn)值實(shí)線，預(yù)測(cè)為虛線。

圖5 1954—2018年洱海流域大理站干旱等級(jí)

4 結(jié) 論

構(gòu)建了多維Normal Copula模型篩選比較，在二維模型之-max組合時(shí)，預(yù)測(cè)逐日潛在蒸散量ET0的精度最高，進(jìn)一步將-max與其他氣象因子構(gòu)建的三維Normal Copula模型，預(yù)測(cè)效果還變差了。引入修正函數(shù)項(xiàng)消除了Copula模型在1—4月逐日ET0預(yù)測(cè)值系統(tǒng)性偏小的現(xiàn)象，相對(duì)誤差小于10%、15%、20%、25%樣本數(shù)分別提高到71.6%、84.4%、91.4%、96.5%。全年平均的ME為0.17～0.30，RMSE為0.54～0.64，IA為0.98～0.99，NSE為0.90～0.98。

將改進(jìn)后的逐日ET0預(yù)測(cè)方法用于氣象干旱（日SPEI指數(shù)）預(yù)測(cè)評(píng)估，逐日SPEI指數(shù)預(yù)測(cè)值與實(shí)際值相關(guān)系數(shù)達(dá)0.95～0.99，ME為-0.10～0.35，RMSE為0.20～0.30，IA為0.97～0.98，NSE為0.91～0.97。由特豐水年至特枯水年，出現(xiàn)重旱或特旱的天數(shù)有324 d，頻次為17.73%，預(yù)測(cè)的相對(duì)偏差僅1.2%。洱海地區(qū)1-6月無(wú)旱及輕旱頻次為36.51%，中旱頻次為30.37%，重旱及特旱頻次為33.11%，預(yù)測(cè)偏差為2.1%；下半年7-12月無(wú)旱及輕旱頻次為89.73%，中旱頻次為9.07%，重旱及特旱頻次為1.2%，預(yù)測(cè)偏差為1.0%。上半年1-6月預(yù)測(cè)ET0與實(shí)際ET0的誤差比下半年7-12月偏大，即降水量豐富時(shí)段該模型預(yù)測(cè)ET0的精度越高，逐日SPEI預(yù)測(cè)的誤差參數(shù)都優(yōu)于同期ET0的預(yù)測(cè)，原因是ET0的數(shù)學(xué)模型中未含降水因子，而逐日SPEI預(yù)測(cè)模型補(bǔ)充進(jìn)來(lái)，綜合地反映了氣溫、降水等對(duì)干旱指數(shù)的作用。

[1] Allen R G, Pereira L S, Raes D, et al. Crop Evapotranspiration Guideline for Computing Crop Requirement[M]. Rome: FAO Irrigation and Drainage, 1988.

[2] Zhang Yongqiang, Pe?a-Arancibia J L, McVicar T R, et al. Multi-decadal trends in global terrestrial evapotranspiration and its components[J]. Nature: Scientific Report, 2016, 6, 19124. doi: 10.1038/srep19124.

[3] 張強(qiáng)，韓蘭英，王勝，等. 影響南方農(nóng)業(yè)干旱災(zāi)損率的氣候要素關(guān)鍵期特征[J]. 科學(xué)通報(bào)，2018，63(23)：2378-2392.Zhang Qiang, Han Lanying, Wang Sheng, et al. The affected characteristic of key period’s climate factor on the agricultural disaster loss caused by drought in the south China[J]. Chinese Science Bulletin, 2018, 63(23): 2378-2392.(in Chinese with English abstract)

[4] 焦醒，劉廣全，匡尚富，等. Penman-Monteith模型在森林植被蒸散研究中的應(yīng)用[J] .水利學(xué)報(bào)，2010，41(2)：245-252. Jiao Xing, Liu Guangquan, Kuang Shangfu, et al. Review on application of Penman-Monteith Equation to studying forest vegetation evapotranspiration[J]. Journal of Hydraulic Engineering, 2010, 41(2): 245-252. (in Chinese with English abstract)

[5] Liu Xiaoying, Lin Erda. Performance of the Priestley-Taylor equation in the semiarid climate of North China[J]. Agricultural Water Management, 2005, 71(1): 1-17.

[6] 賈悅，崔寧博，魏新平，等. 考慮輻射改進(jìn)Hargreaves模型計(jì)算川中丘陵區(qū)參考作物蒸散量[J]. 農(nóng)業(yè)工程學(xué)報(bào)， 2016， 32(21)： 152-160. Jia Yue, Cui Ningbo, Wei Xinping, et al. Modifying Hargreaves model considering radiation to calculate reference crop evapotranspiration in hilly area of central Sichuan Basin[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(21): 152-160. (in Chinese with English abstract)

[7] Irmak S, Irmak A, Allen R S, et al. Solar and net radiation-based equations to estimate reference evapotranspiration in humid climate[J]. ASCE Journal of Irrigation and Drainage Engineering, 2003, 129(5): 336-347.

[8] Wang Dakang, Zhan Yulin, Yu Tao, et al. Improving meteorological input for surface energy balance system utilizing mesoscale weather research and forecasting model for estimating daily actual evapotranspiration[J]. Water, 2020, 12(9): 1-15.

[9] 尹劍，歐照凡，付強(qiáng)，等. 區(qū)域尺度蒸散發(fā)遙感估算：反演與數(shù)據(jù)同化研究進(jìn)展[J]. 地理科學(xué)，2018，38(3)： 448-456. Yin Jian, Ou Zhaofan, Fu Qiang, et al. Review of current methodologies for regional evapotranspiration estimation: Inversion and data assimilation[J]. Scientia Geographica Sinica, 2018,38(3):448-456.

[10] 楊永剛，崔寧博，胡笑濤，等. 中國(guó)糧食主產(chǎn)區(qū)參考作物蒸散量演變特征與成因分析[J]. 中國(guó)農(nóng)業(yè)氣象，2018，39(4)：245-255. Yang Yonggang, Cui Ningbo, Hu Xiaotao, et al. Spatio-temporal variability and cause analysis of reference crop evapotranspiration in the main grain producing areas of China[J]. Chinese Journal of Agrometeorology, 2018, 39(4): 245-255. (in Chinese with English abstract)

[11] Traore S, Luo Y F, Fipps G. Deployment of artificial neural network for short-term forecasting of evapotranspiration using public weather forecast restricted messages[J]. Agricultural Water Management, 2016, 163: 363-379.

[12] Yoo B H, Kim J, Lee B W, et al. A surrogate weighted mean ensemble method to reduce the uncertainty at a regional scale for the calculation of potential evapotranspiration[J]. Nature: ScientificReports, 2020, 10. doi: 10.1038/s41598-020-57466-0.

[13] 王振龍，顧南，呂海深，等. 基于溫度效應(yīng)的作物系數(shù)及蒸散量計(jì)算方法[J] .水利學(xué)報(bào)， 2019， 50(2)：242-251. Wang Zhenlong, Gu Nan, Lü Haishen, et al. Calculation of crop coefficient and evapotranspiration based on temperature effect[J]. Journal of Hydraulic Engineering, 2019, 50(2): 242-251. (in Chinese with English abstract)

[14] Islam A R, Shen Shuanghe, Yang Shengbin, et al. Predicting design water requirement of winter paddy under climate change condition using frequency analysis in Bangladesh[J]. Agricultural Water Management, 2018: 58-70.

[15] 許紅師，練繼建，賓零陵，等. 臺(tái)風(fēng)災(zāi)害多元致災(zāi)因子聯(lián)合分布研究[J]. 地理科學(xué)， 2018，38(12)：2118-2124. Xu Hongshi, Lian Jijian, Bin Lingling, et al. Joint distribution of multiple typhoon hazard factors[J]. Scientia Geographica Sinica, 2018, 38(12): 2118-2124. (in Chinese with English abstract)

[16] Zhang Ye, Wang Kang, Liu Xinming, et al. The three-dimensional joint distributions of rainstorm factors based on Copula function: A case in Kuandian County, Liaoning Province[J]. Scientia Geographica Sinica, 2017, 37(4): 603-610.

[17] Li Ning, Gu Xiaotian, Liu Xueqin. Return period analysis based on joint distribution of three hazards in dust storm disaster[J]. Advances in Earth Science, 2013, 28(4): 490-496.

[18] 邱康俊，溫華洋，王根. 基于概率密度匹配的潛在蒸散量估算[J]. 江蘇農(nóng)業(yè)科學(xué)，2018，46(16)：204-208. Qiu Kangjun, Wen Huayang, Wang Gen. Estimation of potential evapotranspiration amount based on probability density matching[J]. Jiangsu Agricultural Sciences, 2018, 46(16): 204-208. (in Chinese with English abstract)

[19] Chen Lu, Guo Shenglian. Copulas and Its Application in Hydrology and Water Resources[M]. Singapore, 2019: 97-116.

[20] Scott C W, Scott S, William F, et al. Copula theory as a generalized framework for flow‐duration curve based streamflow estimates in ungaged and partially gaged catchments[J]. Water Resources Research, 2019, 55(11): 9378-9397.

[21] Liu Xianfeng, Zhu Xiufang, Pan Yaozhong, et al. Agricultural drought monitoring: Progress, challenges, and prospects[J]. Journal of Geography Science, 2016, 26(6): 750-767. doi: 10.1007/s11442-016-1297-9.

[22] 陳燕麗，蒙良莉，黃肖寒，等. 基于SPEI 的廣西甘蔗生育期干旱時(shí)空演變特征分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019， 35(14)：149-158. Chen Yanli, Meng Liangli, Huang Xiaohan, et al. Spatial and temporal evolution characteristics of drought in Guangxi during sugarcane growth period based on SPEI[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(14): 149-158. (in Chinese with English abstract)

[23] Vicente-Serrano S M, Beguería S, López-Moreno J I. A multi-scalar drought index sensitive to global warming: The standardized precipitation evapotranspiration index[J]. Journal of Climate, 2010, 23(7): 1696-1718.

[24] Yang Qing, Li Mingxing, Zheng Ziyan，et al. Regional applicability of seven meteorological drought indices in China[J]. Science China(Series D: Earth Sciences), 2017, 60: 745–760.

[25] Wang Fei, Wang Zongmin, Yang Haibo, et al. Study of the temporal and spatial patterns of drought in the Yellow River basin based on SPEI[J]. Science China (Series D: Earth Sciences), 2018, 61: 1098–1111.

[26] 賈艷青，張勃. 基于日SPEI的近55 a西南地區(qū)極端干旱事件時(shí)空演變特征[J].地理科學(xué)，2018，38(3)：474-483. Jia Yanqing, Zhang Bo. Spatial-temporal variability characteristics of extreme drought events based on daily SPEI in the southwest China in recent 55 years[J]. Scientia Geographica Sinica, 2018, 38(3): 474-483. (in Chinese with English abstract)

[27] 中科院南京地理與湖泊研究所編. 中國(guó)湖泊調(diào)查報(bào)告[M]. 北京：科學(xué)出版社，2019.

[28] Fouque J P, Zhou Xianwen. Perturbed Gaussian copula[J]. Advances in Econometrics, 2006, 22:103-121.

[29] Luo X, Shevchenko P V. Bayesian model choice of grouped t-Copula[J]. Methodology & Computing in Applied Probability, 2012, 14(4): 1097-1119.

[30] Genest C, Rémillard B, Beaudoin D. Goodness-of-fit tests for copulas: A review and a power study[J]. Insurance Mathematics & Economics, 2009, 44(2): 199-213.

[31] 溫家興，張?chǎng)?，王云，? 多時(shí)間尺度干旱對(duì)青海省東部農(nóng)業(yè)區(qū)小麥的影響[J]. 灌溉排水學(xué)報(bào)，2016，35(4)：92-97.Wen Jiaxing, Zhang Xin, Wang Yun, et al. Effects of drought in multi-time scale on wheat crop in eastern agricultural region of qinghai province[J]. Journal of Irrigation and Drainage, 2016, 35(4): 92-97. (in Chinese with English abstract)

[32] Ma Zhuguo, Fu Congbin. Global aridification in the second half of the 20th century and its relationship to large-scale climate background[J]. Science in China (Series D: Earth Sciences), 2007, 50(5): 776-788.

[33] Liang Hao, Hu Kelin, Batchelor W D, et al. An integrated soil-crop system model for water and nitrogen management in North China[J]. Nature: Scientific Report，2016，6，25755.

[34] 王林，陳文. 標(biāo)準(zhǔn)化降水蒸散指數(shù)在中國(guó)干旱監(jiān)測(cè)的適用性分析[J]. 高原氣象， 2014， 33(2)：423-431. Wang Lin, Chen Wen. Applicability analysis of standardized precipitation evapotranspiration index in drought monitoring in China[J]. Plateau Meteorology, 2014, 33(2): 423-431. (in Chinese with English abstract)

[35] Gu Shixiang, He Daming, Cui Yuamlai, et al. Temporal and spatial changes of agricultural water requirements in the Lancang River Basin[J]. Journal of Geographic Science，2012, 22(3): 441-450.

Daily reference evapotranspiration and meteorological drought forecast using high-dimensional Copula joint distribution model

Gu Shixiang1,2, Zhao Zhong1, Chen Jing1,2, Chen Jinming1,2,3, Zhang Liudong1

(1.,650021,; 2..,650021,; 3.,,430072,)

A high-dimensional copula function was introduced to construct the joint distribution of meteorological factors that affected by reference evapotranspiration (ET0). Specifically, an attempt was made to reveal the correlation structure between different variables in copula function, thereby to establish the joint distribution model of multiple meteorological factors on daily ET0prediction, and finally to analyze short-term drought level. Daily observation data were collected from Dali meteorological station in Erhai watershed from 1954 to 2018. FAO Penman Monteith equation was used to calculate the standard ET0value for the assessment of forecast precision. T-Tmaxtwo-dimensional normal copula function model was used to predict daily ET0after screening a variety of meteorological factor datasets. The systematic error appeared between January to April was necessary to be corrected, otherwise it can make the predicted value relatively smaller than the standard ET0. The empirical correction function with error curve was used for the daily ET0forecast model, to improve the prediction accuracy, thereby to realize the real-time prediction in irrigated region. The results show that: 1) When combining two meteorological factors of T-Tmax, the two-dimensional normal copula model can achieve the highest prediction accuracy for daily ET0, 71.6%, 84.4%, 91.4% and 96.5%, under the relative errors less than 10%, 15%, 20% and 25%, respectively. The annual compliance index IA range was 0.98- 0.99, the average deviation, ME, was 0.17-0.30, the root of mean square error, RMSE, was 0.54-0.64, and the Nash Sutcliffe efficiency coefficient was 0.90-0.98. 2) The daily ET0prediction method was applied to the prediction and evaluation of daily meteorological drought, taking the daily SPEI index as an example. The correlation coefficient between the prediction value of daily SPEI index and the actual value was 0.95- 0.99, ME was -0.10-0.35, RMSE was 0.20-0.30, IA was 0.97-0.98, NSE was 0.91-0.97, respectively. In the season with more precipitation, the accuracy of Copula function model was higher, and the error parameters of daily SPEI prediction were better, than that of daily ET0prediction. 3) From the extremely wet year to the extremely dry year, the proportion of humid and light drought days decreased from 81.3% to 46.0%, the proportion of medium drought days increased from 10.7% to 27.9%, and the proportion of heavy drought and extreme drought days increased from 8.2% to 26.2%. In the five typical years of annual precipitation frequency,= 5%, 25%, 50%, 75% and 95%, the relative deviation of heavy and extremely drought frequency was 1.5% between the predicted ET0and actual ET0, while reached 1.2% after corrected daily ET0prediction and actual ET0, to evaluate the daily meteorological drought level. 4) The results also revealed that the frequency of non-drought and light drought was 36.51%, the frequency of moderate drought was 30.37%, the frequency of severe drought and extreme drought was 33.11%, and the prediction deviation was 1.1% from January to June, whereas from July to December, the frequency of humid and light drought was 89.73%, the frequency of moderate drought was 9.07%, the frequency of severe drought and extreme drought was 1.2%, the prediction deviation was 1%, indicating the significant characteristics of seasonal drought.

drought; evapotranspiration; SPEI index; high-dimensional Copula function; forecast; Erhai Watershed

顧世祥，趙眾，陳晶，等. 基于高維Copula函數(shù)的逐日潛在蒸散量及氣象干旱預(yù)測(cè)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(9)：143-151.doi：10.11975/j.issn.1002-6819.2020.09.016 http://www.tcsae.org

Gu Shixiang, Zhao Zhong, Chen Jing, et al. Daily reference evapotranspiration and meteorological drought forecast using high-dimensional Copula joint distribution model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(9): 143-151. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.09.016 http://www.tcsae.org

2020-01-15

2020-03-23

云南省應(yīng)用基礎(chǔ)研究重點(diǎn)基金（2017FA022）；國(guó)家自然科學(xué)基金項(xiàng)目（51669035）；云南重點(diǎn)研發(fā)計(jì)劃（科技入滇專項(xiàng)）；國(guó)家瀾湄合作基金項(xiàng)目（2018-1177-02）；云南省創(chuàng)新團(tuán)隊(duì)建設(shè)專項(xiàng)（2018HC024）

顧世祥，博士，教授級(jí)高工，從事農(nóng)業(yè)節(jié)水灌溉理論技術(shù)研究。Email：gushxang@qq.com。

10.11975/j.issn.1002-6819.2020.09.016

S271

A

1002-6819(2020)-09-0143-09