基于游程理論和Copula函數(shù)的遼寧省農(nóng)業(yè)氣象干旱特征研究

侯陳瑤，朱秀芳，肖名忠，肖國峰，陳昌為

(1. 環(huán)境演變與自然災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室，北京師范大學(xué) 地理科學(xué)學(xué)部，北京100875；2. 北京師范大學(xué) 地理科學(xué)學(xué)部遙感科學(xué)與工程研究院，北京 100875；3. 河海大學(xué) 水文水資源學(xué)院，江蘇 南京 210098)

據(jù)統(tǒng)計(jì)，近60年來中國年平均受旱面積2 087萬hm2，平均因旱損失糧食163億kg/年[1]，農(nóng)業(yè)氣象干旱災(zāi)情嚴(yán)峻。快速準(zhǔn)確識別農(nóng)業(yè)氣象干旱事件，分析其發(fā)生的頻率及重現(xiàn)期，有助于提高應(yīng)對農(nóng)業(yè)氣象干旱風(fēng)險的能力，減少農(nóng)業(yè)生產(chǎn)活動的損失。

在干旱頻率和重現(xiàn)期的相關(guān)研究中，最初多以單變量計(jì)算為基礎(chǔ)，而單變量重現(xiàn)期只反應(yīng)單一因素的影響，無法刻畫復(fù)雜的干旱特征變量。后來將多種因素進(jìn)行聯(lián)合分析的Copula方法受到干旱領(lǐng)域研究者的關(guān)注，其能較好的對干旱問題中的多特征變量進(jìn)行擬合，并對多變量進(jìn)行條件概率、重現(xiàn)期分析等[2]。例如李穎等在遼西地區(qū)基于SPI指數(shù)利用游程理論和Copula函數(shù)提取干旱歷時和干旱強(qiáng)度分析農(nóng)業(yè)氣象干旱發(fā)生的聯(lián)合概率和重現(xiàn)期[3]；馮星等利用多種Copula函數(shù)在渭河流域分別構(gòu)造干旱歷時、干旱烈度、干旱峰值之間的兩兩二維關(guān)系，分析各種干旱特征變量組合的Copula擬合效果及水文干旱的特點(diǎn)[4]；李阿龍?jiān)诤幽鲜「鶕?jù)游程理論和Copula多維聯(lián)合分布理論選取SPI為干旱指標(biāo)提取干旱特征變量對農(nóng)業(yè)干旱的危險性進(jìn)行分析[5]?，F(xiàn)有研究中對干旱事件的概率分析和對重現(xiàn)期計(jì)算分析多集中在站點(diǎn)尺度，較少有聯(lián)合概率和重現(xiàn)期的區(qū)域尺度分析，且選用的指標(biāo)多為降水據(jù)平百分率[6]、修正帕默爾干旱指數(shù)(SC-PDSI)[7]和標(biāo)準(zhǔn)化降水指數(shù)(SPI)[3]等。而利用由Vicente-Serrano提出的標(biāo)準(zhǔn)化降水蒸散指數(shù)(SPEI)[8]為農(nóng)業(yè)氣象干旱指標(biāo)較少。該指標(biāo)不僅考慮降水因素，還綜合考慮蒸散對干旱的綜合影響。其中3個月的SPEI-3指數(shù)代表季節(jié)尺度的水分平衡估計(jì)，可以有效表征農(nóng)耕區(qū)土壤墑情狀況[9-10]。

遼寧省是東北地區(qū)重要的糧食主產(chǎn)區(qū)，同時也是我國東北地區(qū)干旱的多發(fā)區(qū)，據(jù)《中國氣象災(zāi)害大典-遼寧卷》記載遼寧省干旱災(zāi)害頻發(fā)，嚴(yán)重影響社會經(jīng)濟(jì)的發(fā)展和人民生活水平的提高[11]。以往對于遼寧省的干旱研究多集中在利用干旱指標(biāo)描述遼寧省不同季節(jié)干旱時空分布特征[12-14]，干旱對于農(nóng)作物的影響[15-16]，或抗旱作物的發(fā)展[17-18]等方向，分析遼寧省不同場景下干旱發(fā)生的概率及重現(xiàn)期的研究少。開展干旱特征變量的聯(lián)合概率分布和重現(xiàn)期分析，有利于認(rèn)清遼寧省的旱災(zāi)頻率和聯(lián)合重現(xiàn)期分布規(guī)律，同時對其他區(qū)域干旱特征識別與分析有借鑒和指導(dǎo)意義。

鑒此，選擇遼寧省為研究區(qū)，基于遼寧省近30年的氣象站點(diǎn)逐月氣候數(shù)據(jù)，通過R語言計(jì)算3個月的SPEI-3指數(shù)作為農(nóng)業(yè)氣象干旱指標(biāo)，利用游程理論和Copula函數(shù)提取干旱特征變量并分別從站點(diǎn)尺度和區(qū)域尺度實(shí)現(xiàn)對遼寧省干旱事件特征變量的頻率特征以及重現(xiàn)期的分析。

1 數(shù)據(jù)來源

研究數(shù)據(jù)為遼寧省27個氣象臺站的逐月氣象數(shù)據(jù)序列，主要包括氣壓、氣溫、相對濕度、降水、蒸發(fā)、風(fēng)向風(fēng)速和日照時數(shù)等。其中，草河口、皮口和長海氣象臺站由于缺失數(shù)據(jù)被剔除，時間跨度為1991年1月到2016年5月，所選站點(diǎn)均經(jīng)過了嚴(yán)格的質(zhì)量檢查和控制，包括極值檢驗(yàn)和時間一致性檢驗(yàn)等，消除了非氣候因素造成的影響，以上數(shù)據(jù)取自中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://cdc.cma.gov.cn)。

2 研究方法

2.1 SPEI指數(shù)

標(biāo)準(zhǔn)化降水蒸散指數(shù)(SPEI)是基于月尺度的水平衡模型，對降水量與潛在蒸散量差值序列的累積概率值進(jìn)行正態(tài)標(biāo)準(zhǔn)化后構(gòu)建的指數(shù)。其分級標(biāo)準(zhǔn)是： (-∞, -2]為特旱，(-2.0, -1.5]重旱，(-1.5, -1.0]為中旱，(-1.0, -0.5]為輕旱，(-0.5, 0.5]為無旱[8]。已有學(xué)者基于SPEI-PM對東北地區(qū)干旱演變特征分析，結(jié)果表明SPEI-PM指數(shù)可以表征東北地區(qū)的干旱特征[19-20]?？紤]到SPEI最初主要是采用Thornthwaite模型估算潛在蒸散量[21]，在干旱區(qū)有較大的不確定性。本文采用聯(lián)合國糧農(nóng)組織推薦的Penman-Monteith(P-M)模型[22]計(jì)算潛在蒸散。

2.2 游程理論

游程理論也叫輪次理論，是分析時間序列的一種有效方法，Herbst等首先運(yùn)用該方法對干旱進(jìn)行識別[23]。運(yùn)用游程理論進(jìn)行干旱識別[24]，首先給定一個截斷水平即閾值來分離隨時間變化的干旱指標(biāo)序列，當(dāng)干旱指標(biāo)序列在一個或者多個時間內(nèi)出現(xiàn)連續(xù)大于閾值，則出現(xiàn)正游程，反之，出現(xiàn)負(fù)游程。本文以SPEI的分級標(biāo)準(zhǔn)為參考，選取-0.5作為截斷水平進(jìn)行閾值設(shè)定。

2.3 Copula函數(shù)和重現(xiàn)期

Copula函數(shù)是由Sklar和Nelsen分別提出和發(fā)展的聯(lián)合分布函數(shù)，能夠?qū)⒍鄠€一維邊緣分布函數(shù)在[0,1]區(qū)間連接在一起[25]。其中對稱型的Archimedean Copula函數(shù)構(gòu)造簡單、僅含一個參數(shù)、容易求解，已被廣泛的應(yīng)用于水文多變量頻率計(jì)算[26]。常用的對稱式Archimedean Copula主要包括Gumbel、Clayton和Frank Copula。本文選擇這三種函數(shù)分析的干旱特征變量之間的關(guān)系。

(1)

當(dāng)S?s的條件時，S的條件概率分布為：

(2)

與此同時，相應(yīng)的條件重現(xiàn)期表示為：

(3)

3 結(jié)果分析

本文選取在干旱事件特征研究中應(yīng)用廣泛的Gamma函數(shù)，對數(shù)正態(tài)分布，威布爾分布以及指數(shù)分布函數(shù)來擬合每個站點(diǎn)的干旱歷時和干旱強(qiáng)度的邊緣分布函數(shù)。采用全極大似然法估計(jì)對邊緣分布函數(shù)的參數(shù)估計(jì)，并使用K-S方法進(jìn)行擬合優(yōu)度檢驗(yàn)，在保證最多站點(diǎn)通過顯著性檢驗(yàn)的基礎(chǔ)上，選用通過顯著性檢驗(yàn)的各站點(diǎn)p值均值最高的分布。最終選擇采用p值均值為0.442的威布爾分布對干旱歷時進(jìn)行擬合，采用p值均值為0.793的對數(shù)正態(tài)分布對干旱強(qiáng)度進(jìn)行擬合。所有站點(diǎn)的Spearman相關(guān)系數(shù)均大于0.9以上，符合Copula函數(shù)在聯(lián)合分布研究中對兩組數(shù)據(jù)變量之間具有的相關(guān)性的要求。

基于干旱歷時和干旱強(qiáng)度的邊緣分布函數(shù)，建立了干旱歷時與干旱強(qiáng)度之間的三種Copula模型，利用極大似然法對Copula函數(shù)中未知參數(shù)進(jìn)行了估計(jì)。通過計(jì)算理論Copula與經(jīng)驗(yàn)Copula值間的歐式距離，AIC指標(biāo)和BIC指標(biāo)，對模型的擬合優(yōu)度進(jìn)行評價，數(shù)值越小，模型擬合度越好。所有站點(diǎn)的三種評價指標(biāo)均值都表明Frank-Copula擬合最好。所以本文選取Frank-Copula進(jìn)行聯(lián)合概率和重現(xiàn)期的分析。

3.1 站點(diǎn)尺度遼寧省干旱特征分析

經(jīng)過計(jì)算發(fā)現(xiàn)各站點(diǎn)的二維聯(lián)合分布和聯(lián)合重現(xiàn)期的變化趨勢基本一致。因此，本文以位于遼西地區(qū)的綏中和遼中地區(qū)的沈陽兩個站點(diǎn)為例，對比分析站點(diǎn)尺度識別出的干旱事件特征對應(yīng)的二維聯(lián)合分布以及重現(xiàn)期的變化。兩個站點(diǎn)有著相同的干旱歷時和干旱強(qiáng)度范圍，便于對比分析在相同的干旱歷時和干旱強(qiáng)度的情況下，不同站點(diǎn)的聯(lián)合概率分布和聯(lián)合重現(xiàn)期。

表1顯示了在單變量重現(xiàn)期分別為2、5、10、20和50年一遇的情況下，綏中站和沈陽站對應(yīng)的干旱歷時和干旱強(qiáng)度的值，以及在相同的干旱歷時和干旱強(qiáng)度的組合下，兩站點(diǎn)的特征變量的聯(lián)合重現(xiàn)期。

例如當(dāng)單變量重現(xiàn)期為2年時，綏中站對應(yīng)的干旱歷時為3.154月，沈陽站對應(yīng)的干旱歷時為3.033月，綏中站對應(yīng)的干旱強(qiáng)度為2.909，沈陽站對應(yīng)的干旱強(qiáng)度為2.807。在此干旱強(qiáng)度和干旱歷時的組合下，綏中站和沈陽站的聯(lián)合重現(xiàn)期分別為1.760年和1.744年。由表1可知，在相同的干旱歷時和干旱強(qiáng)度條件下，兩站點(diǎn)的聯(lián)合重現(xiàn)期都比單變量重現(xiàn)期小，這是由于計(jì)算單變量重現(xiàn)期只考慮單一因素，而聯(lián)合重現(xiàn)期考慮兩個因素共同的影響。

圖1為綏中站(圖1a)和沈陽站(圖1b)聯(lián)合概率分布等值線圖。隨著干旱歷時和干旱強(qiáng)度值的增大，兩者的聯(lián)合累積概率值在不斷增大，其中綏中站干旱歷時小于4個月，干旱強(qiáng)度小于3時，發(fā)生干旱的聯(lián)合累積概率增長較為迅速，在這種情況下，有大約70%的干旱事件發(fā)生；沈陽站干旱歷時小于4個月，干旱強(qiáng)度小于2時，發(fā)生干旱的聯(lián)合累積概率增長較為迅速，在這種情況下，有大約50%的干旱事件發(fā)生。

圖2為綏中站和沈陽站的聯(lián)合重現(xiàn)期及聯(lián)合重現(xiàn)期等值線圖。隨著干旱歷時和干旱強(qiáng)度的增大，兩個站點(diǎn)的干旱事件發(fā)生的重現(xiàn)期在逐漸增大，其中綏中站在干旱歷時大于4個月并且干旱強(qiáng)度大于4時，發(fā)生干旱事件的重現(xiàn)期增長迅速(圖2a和圖2c)；沈陽站在干旱歷時大于4個月并且干旱強(qiáng)度大于3時發(fā)生干旱事件的重現(xiàn)期增長迅速(圖2b和圖2d)。當(dāng)干旱歷時和干旱強(qiáng)度同為相同的條件下的最大值時，綏中站的重現(xiàn)期為25年，沈陽站的重現(xiàn)期為15年(圖2a和圖2b)，說明綏中地區(qū)在同等高強(qiáng)度長歷時的條件下發(fā)生重旱的頻率較低。綏中站的歷史干旱事件集中發(fā)生在干旱歷時小于3個月，干旱強(qiáng)度小于4的情況下，沈陽站的歷史干旱事件集中發(fā)生在干旱歷時小于5個月，干旱強(qiáng)度小于5的情況下(圖2c和圖2d)。

3.2 區(qū)域尺度遼寧省干旱特征分析

參考左冬冬等對于干旱歷時和干旱強(qiáng)度的劃分[27]，將遼寧省的干旱劃分為16種情況：月內(nèi)輕旱，月內(nèi)中旱，月內(nèi)重旱，月內(nèi)特旱，季內(nèi)輕旱，季內(nèi)中旱，季內(nèi)重旱，季內(nèi)特旱，跨季輕旱，跨季中旱，跨季重旱，跨季特旱，半年以上輕旱，半年以上中旱，半年以上重旱，半年以上特旱。

基于遼寧省各站點(diǎn)的Copula聯(lián)合累積分布，分別計(jì)算在16種情況下各個站點(diǎn)發(fā)生干旱事件的概率，利用ArcGIS中IDW插值實(shí)現(xiàn)空間化，得到遼寧省各區(qū)域在不同情境下的干旱事件發(fā)生的概率分布圖3。就整個遼寧省區(qū)域而言，季內(nèi)中旱、跨季中旱、月內(nèi)輕旱、半年以上重旱、月內(nèi)中旱和季內(nèi)輕旱的情況出現(xiàn)的概率較高。在這六種情況中，遼寧省發(fā)生干旱情況出現(xiàn)概率最高到最低的依次排序?yàn)椋杭緝?nèi)中旱、跨季中旱、月內(nèi)輕旱、半年以上重旱、月內(nèi)中旱、季內(nèi)輕旱。季內(nèi)中旱、月內(nèi)輕旱和月內(nèi)中旱的情況下，遼西地區(qū)的發(fā)生干旱的概率最高。

基于遼寧省各站點(diǎn)的聯(lián)合重現(xiàn)期，計(jì)算在各種情況下每個站點(diǎn)的最大聯(lián)合重現(xiàn)期和最小聯(lián)合重現(xiàn)期，利用ArcGIS中的IDW插值得到遼寧省各區(qū)域在不同情境下的最大聯(lián)合重現(xiàn)期分布圖4和最小聯(lián)合重現(xiàn)期分布圖5。干旱歷時半年以上，干旱強(qiáng)度大于2的情況下不存在最大聯(lián)合重現(xiàn)期。由圖4知，當(dāng)干旱強(qiáng)度屬于“輕旱”，隨著干旱歷時的增加，遼寧省整個區(qū)域的最大聯(lián)合重現(xiàn)期增加，而當(dāng)干旱強(qiáng)度分別屬于“中旱”和“重旱”的情況下，隨著干旱歷時的增加，遼寧省整個區(qū)域的最大聯(lián)合重現(xiàn)期沒有明顯變化。而當(dāng)干旱歷時分別屬于“月內(nèi)”、“季內(nèi)”和“跨季”時，隨著干旱強(qiáng)度的增加，遼寧省各區(qū)域的最大聯(lián)合重現(xiàn)期表現(xiàn)出明顯的增加趨勢。說明當(dāng)干旱歷時一致時，干旱強(qiáng)度對于各種情況下的最大聯(lián)合重現(xiàn)期更敏感。由圖5知，當(dāng)干旱強(qiáng)度屬于“輕旱”和“中旱”時，隨著干旱歷時的增加，遼寧省各區(qū)域的最小聯(lián)合重現(xiàn)期增加；當(dāng)干旱強(qiáng)度屬于“中旱”和“特旱”時，遼寧省各區(qū)域的最小聯(lián)合重現(xiàn)期沒有明顯變化；當(dāng)干旱歷時分別屬于“月內(nèi)”、“季內(nèi)”、“跨季”和“半年以上”時，隨著干旱強(qiáng)度的增大，遼寧省各區(qū)域的最小聯(lián)合重現(xiàn)期明顯增大。這說明當(dāng)干旱歷時一致時，干旱強(qiáng)度對于各種情況下的最小聯(lián)合重現(xiàn)期更敏感。

圖1 典型站點(diǎn)聯(lián)合概率分布等值線圖

圖2 典型站點(diǎn)聯(lián)合重現(xiàn)期和聯(lián)合重現(xiàn)期等值線

4 結(jié)論和討論

本文利用遼寧省24個氣象站點(diǎn)近30年逐月降水?dāng)?shù)據(jù)，SPEI-3基于游程理論從序列中分離出干旱事件，再利用Copula函數(shù)建立的二維聯(lián)合分布分析遼寧省干旱特征并得到結(jié)果如下：

圖3 不同情境下的干旱事件發(fā)生概率分布圖

圖4 不同情境下的最大聯(lián)合重現(xiàn)期分布圖

圖5 不同情境下的最小聯(lián)合重現(xiàn)期分布圖

(1)遼寧省干旱事件中的特征變量干旱歷時符合威布爾分布，干旱強(qiáng)度符合對數(shù)正態(tài)分布。通過多種Copula函數(shù)進(jìn)行了遼寧省干旱特征變量聯(lián)合分布特征的分析，F(xiàn)rank-Copula函數(shù)模擬效果最好。

(2)站點(diǎn)尺度的分析顯示所有站點(diǎn)的二維聯(lián)合概率分布和聯(lián)合重現(xiàn)期的趨勢總體保持一致。在相同的干旱歷時和干旱強(qiáng)度的條件下，各站點(diǎn)的聯(lián)合重現(xiàn)期都比單變量重現(xiàn)期小。干旱事件多集中在干旱歷時小于4個月，干旱強(qiáng)度小于3的情況下發(fā)生。

(3)區(qū)域尺度的分析顯示遼寧省發(fā)生干旱情況出現(xiàn)概率最高到最低的排序?yàn)椋杭緝?nèi)中旱>跨季中旱>月內(nèi)輕旱>半年以上重旱>月內(nèi)中旱>季內(nèi)中旱。在季內(nèi)中旱、月內(nèi)輕旱和月內(nèi)中旱的情況下，遼西地區(qū)發(fā)生干旱的概率最高。每個站點(diǎn)的最大聯(lián)合重現(xiàn)期和最小聯(lián)合重現(xiàn)期表現(xiàn)為：當(dāng)干旱歷時保持一致，隨著干旱強(qiáng)度的增加，其聯(lián)合重現(xiàn)期有明顯增大的趨勢，說明干旱強(qiáng)度在干旱歷時保持一致時，對各種情況的最大最小聯(lián)合重現(xiàn)期更敏感；當(dāng)干旱強(qiáng)度保持一致，干旱強(qiáng)度為輕旱時，隨著干旱歷時增加，最大最小聯(lián)合重現(xiàn)期增大，屬于中旱時，隨著干旱歷時增加，最小聯(lián)合重現(xiàn)期增加，為其它干旱強(qiáng)度等級時，最大最小聯(lián)合重現(xiàn)期沒有明顯變化。

在已有研究中，不同指數(shù)在干旱等級和變化特征中有著明顯差異，在區(qū)域尺度尤其明顯[28]。本文選取是3個月尺度的SPEI，主要表征季節(jié)尺度內(nèi)水分虧缺情況，適用于反映農(nóng)業(yè)干旱情況[29]。如果選取其他尺度的SPEI指數(shù)，干旱事件的識別將發(fā)生變化。同時，在基于游程理論的干旱事件識別與分析中，閾值選取沒有統(tǒng)一的規(guī)定，已有研究中，學(xué)者多采取單個閾值[27]，進(jìn)行干旱事件的判別，容易降低識別干旱事件的精確度，所以可以進(jìn)一步優(yōu)化閾值設(shè)置的方法[30]。本文僅針對干旱歷時和干旱強(qiáng)度的兩個特征變量進(jìn)行聯(lián)合，現(xiàn)有研究已有選用多個變量[31-32]，而隨著干旱特征變量的增加，多變量聯(lián)合函數(shù)的結(jié)構(gòu)也會更加復(fù)雜。干旱特征對于時間和空間有著較強(qiáng)的敏感性，不同的區(qū)域的時空特性使得Copula函數(shù)對于干旱的應(yīng)用問題具有時變性[2]。所以，干旱指標(biāo)選取、閾值優(yōu)化、Copula多維特征變量分析和時變性探究都將是未來的研究重點(diǎn)。