基于Copula函數的陜西蘋果晚霜凍特征分析

姜 元，茹曉雅，羅 琦，李美榮，王 釗，王景紅，馮 浩，張 東，蘇寶峰，于 強，何建強

·農業生物環境與能源工程·

基于Copula函數的陜西蘋果晚霜凍特征分析

姜 元1,2，茹曉雅1,2，羅 琦1,2，李美榮3，王 釗3，王景紅3，馮 浩2,4，張 東5，蘇寶峰6，于 強4，何建強1,2※

（1. 西北農林科技大學旱區農業水土工程教育部重點實驗室，楊凌 712100；2. 西北農林科技大學中國旱區節水農業研究院，楊凌 712100；3. 陜西省農業遙感與經濟作物氣象服務中心，西安 710015；4. 中國科學院水利部水土保持研究所黃土高原土壤侵蝕與旱地農業國家重點實驗室，楊凌 712100；5. 西北農林科技大學園藝學院，楊凌 712100；6. 西北農林科技大學農業農村部農業物聯網重點實驗室，楊凌 712100）

為探討利用Copula函數對蘋果晚霜凍進行特征分析的適用性，該研究首先基于陜西省蘋果主產區7個氣象站點1971－2018年的逐日最低氣溫（min）數據集，提取出晚霜凍事件的歷時和強度兩個特征變量。然后，基于6種不同的Copula函數構建晚霜凍特征變量的聯合分布，并進行擬合優度評價。最后，利用優選的Copula函數分析晚霜凍發生的概率及重現期。結果表明：陜西蘋果產區各站點1971－2018年受晚霜凍的影響在空間分布上由東南向西北方向加重，各站點晚霜凍的歷時和強度之間均具有顯著的正相關關系。當晚霜凍強度和歷時增大時，其聯合累積概率也相應增大，且增大趨勢變緩。各站點聯合重現期代表的“或”事件比同現重現期所代表的“且”事件更容易發生。當單變量重現期取值較小時，可將聯合重現期和同現重現期視為單變量重現期的兩種極端情況，對其實際范圍進行估計。總體而言，陜西蘋果產區各站點發生長歷時且高強度晚霜凍事件的概率較小，但是位于延安果區的站點相較于其他果區站點更容易發生高強度或長歷時的晚霜凍事件，以及高強度長歷時同時發生的晚霜凍事件，需要重點加以關注。該研究可為陜西蘋果產區應對晚霜凍災害提供理論依據。

氣象；災害；晚霜凍；Copula函數；頻率分析；重現期；蘋果

0 引 言

蘋果在中國的栽培面積和產量均位居世界首位[1-2]。蘋果的廣泛栽培使得蘋果生產更易受到各種氣象災害的影響，晚霜凍是中國北方蘋果產區面臨的首要極端氣象災害[3-5]，會對蘋果產量和質量有嚴重影響[4,6]，給當地蘋果產業造成重大經濟損失[4]。目前，全球氣候變暖的形勢日趨嚴峻[7]，冬季增溫幅度也更為明顯[8]，暖冬導致果樹春季萌芽及開花提前[9-10]，使其遭受晚霜凍的風險日益增加。陜西地區的大陸性季風氣候特征十分明顯，全年約有40%的強降溫天氣發生在3－4月份。然而，該時期正處于蘋果開花期，頻繁的冷空氣過程導致晚霜凍災害風險較高[11-13]。

目前，關于蘋果晚霜凍的研究多集中在建立蘋果晚霜凍指數等級并探究晚霜凍事件時空變化特征和成災因子。例如，李健等[14]基于氣象站最低氣溫構建了蘋果果樹花期凍害預警指標，并分析了主果區花期凍害的時空變化特征。屈振江等[13]構建了凍害風險災損率對蘋果花期凍害風險進行了評估。柏秦鳳等[15]基于花期凍害風險指數對蘋果花期凍害進行降尺度風險區劃。王景紅等[12]利用人工氣候箱模擬低溫并結合歷史災情調查對蘋果花期霜凍指標進行了修訂優化。張琪等[3]基于物候模型探究了蘋果霜凍發生特征。邱美娟等[1]結合晚霜凍氣象指標，對蘋果花期晚霜凍氣候風險進行了評估。王明昌等[11]采用線性傾向法和風險指數法分析了陜西省禮泉和旬邑兩縣蘋果花期凍害風險情況。

上述研究大多基于晚霜凍指標進行蘋果晚霜凍災害風險的分析，但是評估霜凍等氣象災害事件風險的較佳方法是通過詳細了解以持續時間和低溫強度為特征的霜凍事件的發生頻率[16]。李美榮等[17]建立了蘋果花期嚴重凍害最低氣溫的概率分布模型及重現期預測，但是只使用了單變量頻率分析方法且僅考慮了影響霜凍的一個變量，無法捕捉到極端氣象災害問題的復雜性。

另一種評估氣象災害事件風險更為有效的的方法是建立多元聯合分布[18]。然而，許多傳統多元聯合分布要求各變量的邊緣分布屬于同一類型[19]。Copula函數作為一種構造靈活的模型不限制各變量邊緣分布的選擇，能夠通過邊緣分布和相關性結構兩部分來構建多維聯合分布，形式靈活多樣，被廣泛應用于降雨頻率[20]、干旱頻率[21]和洪水頻率[22]等水文和災害事件分析。Chatrabgoun等[16]首次采用Copula函數對葡萄園霜凍進行了概率評估，通過對霜凍持續時間和嚴重程度這2個特征變量建立聯合分布，探究了位于伊朗馬來爾地區葡萄園霜凍現象的風險和影響，表明Copula函數可以作為構建霜凍多維特征變量聯合分布函數的有效工具。然而，將Copula函數應用于蘋果晚霜凍特征分析的研究目前還鮮有報道。

本研究以陜西省蘋果主產區7個地面氣象觀測站1971－2018年的氣象數據集為基礎數據，采用Copula函數對陜西蘋果產區的晚霜凍災害進行概率評估，研究目標包括：1）逐個氣象站提取晚霜凍事件發生的連續天數和日最低溫度最小值的絕對值分別作為表征晚霜凍歷時和強度的2個特征變量；2）對2個特征變量的單變量邊緣分布進行擬合，并分析變量間相關性；3）基于6種Copula函數，建立2個特征變量的多元聯合分布模型并進行擬合優度檢驗；4）分析計算聯合累積概率和重現期，從而對陜西蘋果產區各站點的晚霜凍災害風險進行評估。本研究的結果有望在站點尺度上為陜西蘋果晚霜凍災害防御和蘋果產業生態管理提供科學依據。

1 材料和方法

1.1 研究區概況與數據來源

陜西蘋果主產區位于黃土高原地區，是符合蘋果生長氣象指標的優質產區之一。根據氣候特點及物候差異，將陜西蘋果主產區劃分為延安果區、渭北東部果區、渭北西部果區和關中西部果區[23]（圖1）。本研究使用的1971－2018年陜西省氣象站點的逐日最低氣溫（min）數據來源于中國氣象科學數據共享服務網（http://cdc.cma.gov.cn/）。使用該氣象數據前需對各個站點進行檢查，刪除數據序列長度過短或缺失的站點，并對數據異常值進行處理。最終選擇研究區域內的延安市寶塔區、洛川、長武、旬邑、銅川、白水和禮泉等7個氣象站作為研究站點（圖1）。由于陜西蘋果花期物候觀測數據不夠完整，所以本文使用的蘋果花期物候數據主要來自王潤紅等[24]基于物候模型對陜西省蘋果花期進行模擬所得的數據。

1.2 研究方法

1.2.1 蘋果晚霜凍事件識別

劉映寧等[26]根據蘋果晚霜凍氣象指標，結合災情調查資料，發現蘋果花期內不同程度的低溫會對蘋果生長發育產生不同的影響：日最低溫度低于?4 ℃時，出現嚴重晚霜凍，中心花受凍率高達80%，減產30%以上；日最低溫度在?2至?4 ℃時，出現中度晚霜凍，中心花受凍率達60%至80%，對產量、品質、商品率產生嚴重影響；日最低溫度在0至?2 ℃時，出現輕度晚霜凍，中心花受凍率達30%至60%，對產量、品質、商品率產生明顯影響。此外，相關文獻也給出了蘋果花期凍害的臨界溫度。李美榮等[26]將日最低氣溫min為0、?2、?4 ℃作為蘋果在花期內不同時期受凍的臨界溫度；劉映寧等[26]提出把日最低氣溫分別低于0、?2、?4 ℃作為陜西蘋果花期凍害農業保險的3個等級；王景紅等[27]則采用0、?2、?4 ℃作為陜西蘋果不同等級花期凍害的臨界指標。因此，本研究將日最低氣溫min為0 作為發生蘋果晚霜凍事件的臨界溫度，每個晚霜凍事件中日最低溫度最小值的絕對值被作為該期間的晚霜凍強度。蘋果晚霜凍不僅與低溫強度有關，也與低溫的持續時間相關，強度越大或持續時間越長受凍率越高[12]，因此將晚霜凍發生的連續天數定義為晚霜凍歷時。

圖1 陜西蘋果產區亞區劃分及氣象站點分布

1.2.2 構建晚霜凍歷時和晚霜凍強度的邊緣分布函數

在基于Copula函數構建晚霜凍歷時和強度聯合分布之前，首先要確定晚霜凍歷時和強度這2個特征變量的邊緣分布函數，同時要考慮它們之間的相依性。本文選用正態分布、指數分布、伽馬分布、威爾布分布、柯西分布、邏輯斯諦分布、對數正態分布7種概率分布函數構建這2個特征變量邊緣函數。使用R語言的“fitdistrplus”包，利用極大似然法估計相關參數，對晚霜凍歷時和強度這2個特征變量的邊緣分布函數進行擬合，并通過Kolmogorov-Smirnov（K-S）檢驗選取最優的單變量邊緣分布函數。采取3種常用的相關系數來進行2個晚霜凍特征變量間的相關性度量，包括皮爾遜積矩相關系數（Pearson product-moment correlation coefficient）、斯皮爾曼秩相關系數（Spearman’s rank correlation coefficient）、肯德爾秩相關系數（Kendall rank correlation coefficient）。

以陜西省禮泉站為例，因為使用逐日最低氣溫來識別晚霜凍事件，所以晚霜凍歷時數據均為整數天，導致晚霜凍歷時存在較多的相同值，因此晚霜凍歷時的經驗分布呈條帶狀（圖2）。在對晚霜凍歷時邊緣分布函數進行擬合過程中，這種條帶狀的離散分布會導致樣本個數明顯少于實際值，進而影響擬合效果。參考韓會明等[28-29]引入干旱歷時離散變量連續化處理方法，本研究對整數的晚霜凍歷時加上[?0.5, 0.5]的均勻隨機分布變量，從而使離散化的數據序列連續化作為處理后的經驗分布。此外，Michele等[30]證明，對樣本添加均勻分布的擾動不會改變原始的統計量信息。

圖2 禮泉站晚霜凍歷時邊緣分布擬合情況

1.2.3 構建晚霜凍歷時和晚霜凍強度的Copula聯合分布函數

Copula函數可用于構造邊緣分布不同的多個變量的聯合分布函數。根據Sklar定理[31]，設=(<)和=(<)分別為特征變量晚霜凍歷時和強度的邊緣分布函數，其聯合分布函數為(,)（式（1））。

(,)=(<,<)=[(),()]（1）

式中為特征變量取值，為Copula函數，為邊緣分布函數。

確定單個變量的邊緣分布函數后，要優選一種Copula函數連接單變量分布函數。本文選用了6種常用的二維Copula聯合分布函數（表1），包括Archimedean Copula函數中的Gumbel Copula，Frank Copula，Clayton Copula以及Joe Copula，和橢圓Copula函數中最常見的Normal Copula和T Copula[16,20-21]。其中，參數估計均選用極大似然法進行估計，擬合優度的檢驗選用赤池信息準則（Akaike Information Criterion, AIC；式（2））和貝葉斯信息準則（Bayesian Information Criterion, BIC；式（3）），以赤池信息量AIC值和貝葉斯信息量BIC值最小作為擬合最優的檢驗標準。

AIC=?2ln()+2（2）

BIC=?2ln()+ln()·（3）

式中為參數的數量；為似然函數；為數據數量。

表1 本研究中所采用的6種二維Copula函數

注：為特征變量取值，為邊緣分布函數，為關聯參數，為生成函數自變量。

Notes:are characteristic variable values;,are the marginal distribution functions;denotes to the association parameters;is the generator function argument.

1.2.4 基于晚霜凍歷時和強度的晚霜凍事件重現期確定

重現期作為水文分析中評估風險的一個通用標準，被廣泛用于確定重大災害事件發生的概率，極具現實意義和應用價值。晚霜凍是一種對蘋果生長極具威脅性的氣象災害，采用重現期對晚霜凍事件進行評估，有助于了解晚霜凍事件的發生規模與發生頻率之間的關系，對晚霜凍事件的發生進行預測，有助于陜西蘋果產區更好地防御晚霜凍災害。一般情況下，重現期的分析僅限于單變量，所謂單變量重現期（式（4））是指變量超過某一特定值出現一次的間隔時間，然而由于極端氣象災害問題的復雜性，導致單變量重現期分析結果具有一定的不準確性[17]。晚霜凍事件是晚霜凍歷時和強度兩個因素相互作用的結果，所以在晚霜凍事件特征分析中應該重點考慮重現期的聯合屬性和條件屬性。對于多特征變量的聯合分布而言，聯合重現期（Joint return period；式（5））代表多個特征變量中某一個特征變量大于給定閾值出現一次的“或”事件的間隔時間，而同現重現期（Co-occurrence return period；式（6））代表多個特征變量中每一個特征變量均大于給定閾值出現一次的“且”事件的間隔時間。

傳統基于單變量重現期的計算公式為

式中F()為變量的邊緣分布函數；單變量同理。

二維聯合重現期和同現重現期的計算公式分別為[32]：

2 結果與分析

2.1 陜西省蘋果主產區晚霜凍事件特征分析

本研究利用陜西省蘋果主產區各個氣象站點數據，提取出1971－2018年間該區域蘋果花期發生晚霜凍事件的晚霜凍歷時和晚霜凍強度（表2），以及每次晚霜凍發生的起始和終止日期。該區域不同果區的站點發生晚霜凍的頻率存在明顯差異，其中位于延安和渭北西部果區的4個站點共發生晚霜凍頻率高，共計266次，且晚霜凍的歷時和強度值也相應較大；渭北東部和關中西部的站點發生晚霜凍共89次，頻次相對較少，且歷時和強度值也較小。在陜西蘋果產區各氣象站中，延安站不僅嚴重晚霜凍和總晚霜凍發生次數最多，而且晚霜凍強度平均值也最大；旬邑站的晚霜凍歷時和強度的最大值均為最大。

基于晚霜凍歷時的起止時間進一步分析晚霜凍歷時過程（圖3），可見總體上陜西省蘋果主產區晚霜凍事件歷時多分布于3月31日至4月16日。關中西部果區的站點（圖3a）明顯發生晚霜凍的次數少、歷時短、強度低；渭北西部果區（圖3d和圖3e）和延安果區（圖3f和圖3g）的站點發生的晚霜凍次數多、歷時長、強度高。陜西蘋果主產區各站點在1980年晚霜凍最多發生達17次；在2013年發生晚霜凍共12次，為近10年晚霜凍發生次數最多的年份。

表2 陜西省蘋果主產區1971－2018年晚霜凍事件統計

圖3 陜西蘋果產區7個氣象站的晚霜凍歷時分布

2.2 晚霜凍歷時和強度邊緣分布函數的擬合

本研究采用正態分布、指數分布、伽馬分布、威布爾分布、柯西分布、邏輯斯諦（Logistic）分布、對數正態分布等7種概率分布，對陜西蘋果產區晚霜凍歷時和強度的邊緣分布函數進行擬合，得到晚霜凍歷時和強度邊緣分布的K-S檢驗結果以及變量之間的相關系數（表3）。選取擬合度最高的概率分布類型（顯著性水平=0.05），最終確定各站點晚霜凍歷時邊緣分布函數為對數正態分布。對于晚霜凍強度這一特征變量的邊緣分布擬合，各站點情況不一，其中禮泉站、白水站和旬邑站為Logistic分布，銅川站、洛川站和延安站為指數分布，長武站為伽馬分布。相關性檢驗結果表明各站點Pearson、Spearman和Kendall這3個相關系數均大于0.4且都通過了顯著性水平=0.05的顯著性檢驗，說明各站點的晚霜凍歷時和強度具有一定的相關性，可基于Copula函數建立聯合分布函數進行進一步分析。

表3 單變量邊緣分布的K-S檢驗和變量間的相關系數

注：*表示通過了顯著性水平= 0.05的顯著性檢驗，**表示通過了顯著性水平= 0.01的顯著性檢驗。

Notes: * is statistically significance level of α= 0.05; ** is statistically extreme significant level of= 0.01.

2.3 基于Copula函數晚霜凍事件聯合分布函數的建立

本研究基于擬合的晚霜凍歷時和強度的邊緣分布函數，建立了晚霜凍歷時和強度之間的6種Copula分布函數，利用極大似然法對其中未知參數進行了估計。此外，還計算了相應的AIC、BIC值（式（2）和式（3）），對Copula分布函數擬合優度進行評價（圖4）。根據AIC和BIC準則判斷可知，Normal Copula函數對禮泉站和旬邑站晚霜凍事件的擬合效果最好；Clayton Copula函數對白水站和銅川站晚霜凍事件的擬合效果最好；對洛川站晚霜凍事件擬合最優的是Frank Copula函數；對延安站和長武站晚霜凍事件擬合優度最高的是Joe Copula函數。

圖4 陜西蘋果產區晚霜凍歷時和強度的6種Copula分布函數擬合優度檢驗

2.4 晚霜凍歷時和強度的聯合概率

進一步進行本研究中各站點晚霜凍特征變量二維聯合概率分布分析（圖5），可知當晚霜凍強度和歷時值增大時，其聯合累積概率也相應增大。總體而言，各站點所在區域發生短歷時低強度、短歷時高強度、長歷時低強度的晚霜凍事件的概率較大，而同時滿足長歷時和高強度的晚霜凍事件的發生概率較小。具體地，在禮泉站（圖5a），晚霜凍強度低于3 ℃、晚霜凍歷時在0.5至2 d時，其聯合累積概率隨特征變量的增大而迅速增大；在白水站（圖5b），晚霜凍強度高于3 ℃、晚霜凍歷時高于3 d時，其聯合累積概率增大趨勢明顯變緩；在銅川站（圖5c），晚霜凍強度不超過4℃、晚霜凍歷時在0.5至3 d時，其聯合累積概率隨特征變量的增大而迅速增大；在旬邑站（圖5d），晚霜凍強度不超過2 ℃、晚霜凍歷時不超過4 d時，其聯合累積概率隨特征變量的增大而迅速增大；在長武站（圖5e），晚霜凍強度超過1 ℃、晚霜凍歷時超過2 d時，其聯合累積概率增大趨勢明顯變緩；在洛川站（圖5f），晚霜凍強度不超過2 ℃、晚霜凍歷時不超過3 d時，其聯合累積概率隨特征變量的增大而迅速增大；在延安站（圖5g），晚霜凍強度超過1 ℃、晚霜凍歷時超過2 d時，其聯合累積概率增大趨勢明顯變緩。

圖5 陜西蘋果產區7個氣象站晚霜凍歷時和強度聯合概率分布

2.5 晚霜凍歷時和強度的重現期

以等值線圖的形式對晚霜凍事件的聯合重現期進行描繪（圖6）。根據晚霜凍兩個特征變量值垂直相交的交點所在等值線可查出聯合重現期，可發現隨著晚霜凍特征變量不斷增大，兩變量的聯合重現期也隨之增大。根據等線圖的疏密情況可明顯看到，不同站點間晚霜凍事件的聯合重現期存在一定的差異，各站點中禮泉站最不容易發生高強度或長歷時的晚霜凍事件（圖6a），而白水站和銅川站較不易發生高強度或長歷時的晚霜凍事件（圖6b和6c）。對于其他站點，若霜凍強度達到8 ℃或者晚霜凍歷時達到8 d時，在旬邑站其聯合重現期約為150 a（圖6d），在長武站其聯合重現期約為280 a（圖6e），在洛川站其聯合重現期約為150 a（圖6f），在延安站其聯合重現期約為50 a（圖6g）。這說明相較于旬邑站、長武站和洛川站，延安站更容易受到高強度或長歷時晚霜凍事件的威脅，也進一步表明聯合重現期等值線圖可以表征不同區域間的晚霜凍事件發生概率，能夠服務于減災預防工作。

通過單變量重現期（式（4））以及構建的單變量邊緣分布的逆函數，分別計算當單變量重現期為3、5、10、20、50、100 a時，各站點晚霜凍歷時和強度值，再計算聯合重現期（式（5））和同現重現期（式（6））對應值（表4）。結果發現聯合重現期和同現重現期均隨著晚霜凍特征變量的增大而增大，這與等值線圖（圖6）所示的聯合重現期變化規律相一致。在單變量取值同等增幅條件下，同現重現期的增幅要明顯高于聯合重現期，這說明該區域各站點聯合重現期代表的“或”事件比同現重現期所代表的“且”事件更容易發生。然而，不同站點的重現期變化存在一定的差異，例如，當單變量重現期增幅相同時，延安站對應的同現重現期增幅明顯大于旬邑站，這說明延安站更易受到長歷時且高強度晚霜凍事件的威脅。各站的單變量重現期、聯合重現期和同現重現期具有相同的大小順序，單變量重現期總介于聯合重現期和同現重現期之間，因此可以通過計算聯合重現期和同現重現期來估計單變量重現期的范圍，即可得到不同晚霜凍特征變量所代表的晚霜凍事件發生的頻率，但是僅當單變量重現期較小時估計的范圍才比較準確。

圖6 陜西蘋果產區7個氣象站晚霜凍歷時和強度聯合重現期等值線

表4 陜西蘋果產區7個氣象站晚霜凍歷時和強度的單變量重現期和聯合分布重現期

續表

3 討 論

在對極端氣象災害事件進行特征分析時，首先要明確氣象災害指標以及事件相關的特征變量，要能夠對極端氣象災害事件進行綜合性評價。以往對晚霜凍災害的研究重點是探究低溫強度的影響，然而張曉煜等[33]提出低溫和低溫持續時間是對霜凍危害程度起決定作用的關鍵因子；王景紅等[12]表示蘋果晚霜凍不僅與低溫強度有關，也與低溫持續時間密切相關，溫度越低或持續時間越長，蘋果受凍率越高；袁佰順等[6]發現0 以下溫度持續時間和低溫強度都需作為判斷霜凍受災和受災程度的因素。因此，本研究選取晚霜凍歷時和晚霜凍強度作為特征變量，基于Copula函數建立特征變量的多元聯合分布，克服了單變量頻率分析方法不能捕捉到極端氣象災害復雜性的不足，為陜西蘋果產區晚霜凍特征研究提供了新的思路和方法。

本研究通過對陜西蘋果產區晚霜凍事件的特征分析可知，在空間上陜西蘋果產區受晚霜凍的影響總體上呈現東南向西北加重的趨勢，這與前人對陜西蘋果晚霜凍的探究結果基本一致[13-14,17,25-26,34]。該空間分布趨勢是由于延安果區和渭北果區所處的經緯度、高海拔及地勢起伏大所導致的。李美榮等[17]對旬邑站進行單變量重現期計算，得到10 a重現期對應的晚霜凍強度為5.2 ℃，而本研究為3.81 ℃，這可能是由于研究所選的時間尺度不同、識別晚霜凍所用的氣象指標不同，以及擬合分布類型不同所導致的。此外，本研究采用Copula函數通過對晚霜凍歷時和強度這2個特征變量建立多元聯合分布評估陜西蘋果晚霜凍風險，相比于李美榮等[17]僅考慮低溫強度而使用單變量頻率的分析方法，更能反映晚霜凍事件受晚霜凍歷時和強度兩個特征變量共同影響的復雜性，從而為精準評估陜西蘋果晚霜凍風險提供更為全面和科學的信息。

準確的蘋果物候數據是評估蘋果花期凍害風險的重要依據。然而，目前國內關于蘋果物候觀測數據的積累十分有限，存在時間序列短、缺測嚴重等問題[11]，并且現有的驅動物候模型模擬的氣象數據也是有限的，從中國氣象科學數據共享服務網所獲得的氣象數據基本只能到縣域尺度，缺乏空間密集的氣象觀測數據，難以對蘋果晚霜凍事件進行區域精準識別，所以未來研究要加強對陜西蘋果產區蘋果花期的多站點和時序觀測，采用質量較高的氣象數據產品來驅動物候模型進行模擬。蘋果屬于多年生作物，不同地區、不同樹齡、不同品種的蘋果對氣候響應有所不同，然而目前晚霜凍指標沒有進行針對性的細化，使得目前對晚霜凍事件的識別工作具有一定的局限性。本文僅對晚霜凍歷時和強度這2個特征變量進行聯合，后續研究可以增加對晚霜凍特征量的選取，但是多元聯合分布函數也會隨著變量維數增多變得愈發復雜。所以，蘋果花期模擬、晚霜凍指標細化、晚霜凍事件精準識別、基于Copula函數的更高維特征變量分析等都是今后進行蘋果晚霜凍災害研究需要關注的重點方向。

4 結 論

本研究根據日最低氣溫閾值提取1971–2018年間陜西蘋果產區各氣象站晚霜凍事件的歷時和強度，基于Copula函數建立了晚霜凍歷時和強度的聯合分布，并分析了研究區內晚霜凍事件的重現期概況，得到以下主要結論：

1）位于延安和渭北西部果區的站點共發生晚霜凍266次，是主要的晚霜凍災害頻發區，并且歷時長、強度高；渭北東部和關中西部的站點發生晚霜凍共89次，頻次相對較少，而且歷時短、強度低。陜西蘋果產區各站點受晚霜凍的影響在空間分布上呈由東南向西北加重的趨勢。

2）陜西蘋果產區各站點晚霜凍歷時的最優邊緣分布均為對數正態分布，晚霜凍強度的最優邊緣分布隨站點而各異。各站點晚霜凍的歷時和強度之間的相關系數均大于0.4，具有顯著的正相關關系。

3）本研究采用各站點擬合優度最高的Copula 函數建立陜西蘋果產區各站點晚霜凍歷時和強度的聯合分布函數，分析晚霜凍特征變量的聯合累積概率可知，同時滿足長歷時和高強度的晚霜凍事件的發生概率較小。

4）不同站點間晚霜凍事件的重現期存在一定的差異，若霜凍強度達到8℃或者晚霜凍歷時達到8 d時，延安站聯合重現期約為50 a，相較于其他站點百年以上的聯合重現期，延安站更容易受到高強度或長歷時，以及高強度且長歷時晚霜凍事件的威脅。陜西蘋果產區各站點聯合重現期代表的“或”事件比同現重現期所代表的“且”事件更容易發生。當單變量重現期的值較小時，可將聯合重現期和同現重現期可看作單變量重現期的兩種極端情況，對單變量重現期的實際范圍進行估計。

總之，本文基于Copula函數建立的聯合分布模型能夠同時描述晚霜凍事件的歷時和強度之間的關系，可以利用重現期評估陜西蘋果產區站點尺度未來晚霜凍發生的風險，從而為研究該區域各站點晚霜凍事件的發生規律，提高該區域各站點晚霜凍的預報能力提供新的科學方法和依據。

[1] 邱美娟，劉布春，劉園，等. 中國北方主產地蘋果始花期模擬及晚霜凍風險評估[J]. 農業工程學報，2020，36(21)：154-163.

Qiu Meijuan, Liu Buchun, Liu Yuan, et al. Simulation of first flowering date for apple and risk assessment of late frost in main producing areas of northern China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(21): 154-163. (in Chinese with English abstract)

[2] 萬煒，師紀博，劉忠，等. 棲霞市蘋果園氮磷養分平衡及環境風險評價[J]. 農業工程學報，2020，36(4)：211-219.

Wan Wei, Shi Jibo, Liu Zhong, et al. Nitrogen and phosphorus nutrient balance and environmental risk assessment of apple orchard in Qixia city[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(4): 211-219. (in Chinese with English abstract)

[3] 張琪，李榮平，張曉月，等. 綏中蘋果霜凍災害特征及其對氣候變暖的適應分析[J]. 北方園藝，2016(13)：30-33.

Zhang Qi, Li Rongping, Zhang Xiaoyue, et al. Analysis of apple frost disaster characteristics and adapt to the climate warming in Suizhong[J]. Northern Horticulture, 2016(13): 30-33. (in Chinese with English abstract)

[4] 馬延慶，劉長民，朱海利，等. 陜西咸陽渭北旱塬地區優質蘋果基地生態氣候特征分析[J]. 干旱地區農業研究，2008(1)：146-153.

Ma Yanqing, Liu Changmin, Zhu Haili, et al. The analysis of ecosystem characteristic of high quality apple producing base on Weibei arid plain region in Xianyang of Shaanxi Province[J]. Agricultural Research in the Arid Areas, 2008(1): 146-153. (in Chinese with English abstract)

[5] 李化龍，趙西社. 陜西黃土高原果業氣候生態條件研究及應用[M]. 北京：氣象出版社，2010.

[6] 袁佰順，許彥平，姚曉紅，等. 基于春季晚霜凍害的天水蘋果氣候風險區劃[J]. 干旱區資源與環境，2015，29(2)：185-189.

Yuan Baishun, Xu Yanping, Yao Xiaohong, et al. The compartment of later frost climate risk for apple in Tianshui[J]. Journal of Arid Land Resources and Environment, 2015, 29(2): 185-189. (in Chinese with English abstract)

[7] Organization W M. WMO Provisional Statement on the State of the Global Climate in 2019[EB/OL]. http://www.cma.gov.cn/en2014/climate/featutes/202003/P020200312816904145935.pdf, 2022-6-23.

[8] Cook B I, Wolkovich E M, Parmesan C, Divergent responses to spring and winter warming drive community level flowering trends[J]. 2012, 109(23): 9000-9005.

[9] Ge Q, Wang H, Rutishauser T, et al. Phenological response to climate change in China: A meta-analysis[J]. Global Change Biology, 2015, 21(1): 265-274.

[10] 郭佳，張寶林，高聚林，等. 氣候變化對中國農業氣候資源及農業生產影響的研究進展[J]. 北方農業學報，2019，47(1)：105-113.

Guo Jia, Zhang Baolin, Gao Julin, et al. Advances on the impacts of climate change on agro-climatic resources and agricultural production in China[J]. Journal of Northern Agriculture, 2019, 47(1): 105-113. (in Chinese with English abstract)

[11] 王明昌，劉布春，劉園，等. 陜西蘋果主產縣花期凍害風險評估[J]. 中國農業氣象，2020，41(6)：381-392.

Wang Mingchang, Liu Buchun, Liu Yuan, et al. Assessment on the freezing injury risk during apple flowering in Liquan and Xunyi[J]. Chinese Journal of Agrometeorolog, 2020, 41(6): 381-392. (in Chinese with English abstract)

[12] 王景紅，劉璐，高峰，等. 陜西富士系蘋果花期霜凍災害氣象指標的修訂[J]. 中國農業氣象，2015，36(1)：50-56.

Wang Jinghong, Liu Lu, Gao Feng, et al. Revision on meteorological indices of florescence frost disaster for Fuji apple in Shaanxi Province[J]. Chinese Journal of Agrometeorology, 2015, 36(1): 50-56. (in Chinese with English abstract)

[13] 屈振江，劉瑞芳，郭兆夏，等. 陜西省蘋果花期凍害風險評估及預測技術研究[J]. 自然災害學報，2013，22(1)：219-225.

Qu Zhenjiang, Liu Ruifang, Guo Zhaoxia, et al. Study of risk assessment and prediction of apple blooming freezing injury in Shaanxi Province[J]. Journal of Natural Disasters, 2013, 22(1): 219-225. (in Chinese with English abstract)

[14] 李健，劉映寧，李美榮，等. 陜西果樹花期低溫凍害特征及防御對策[J]. 氣象科技，2008(3)：318-322.

Li Jian, Liu Yingning, Li Meirong, et al. Low temperature and freezing injury to fruit trees at bloom stage in Shaanxi and countermeasures[J]. Meteorological Science and Technology, 2008(3): 318-322. (in Chinese with English abstract)

[15] 柏秦鳳，王景紅，梁軼，等. 基于縣域單元的降尺度蘋果花期凍害風險區劃[J]. 中國農學通報，2013，29(16)：153-158.

Bai Qinfeng, Wang Jinghong, Liang Yi, et al. The apple flowering frost damage risk zoning down-scaling based on county level[J]. Chinese Agricultural Science Bulletin, 2013, 29(16): 153-158. (in Chinese with English abstract)

[16] Chatrabgoun O, Karimi R, Daneshkhah A, et al. Copula-based probabilistic assessment of intensity and duration of cold episodes: A case study of Malayer vineyard region[J]. Agricultural and Forest Meteorology, 2020, 295: 108150.

[17] 李美榮，李星敏，柏秦鳳，等. 蘋果極端氣象災害氣溫極值的分布及重現期預測[J]. 干旱地區農業研究，2012，30(3)：257-261.

Li Meirong, Li Xingmin, Bai Qinfeng, et al. Probability distribution and recurrence interval prediction of the extreme value of air temperature in extreme meteorological disasters of apples[J]. Agricultural Research in the Arid Areas, 2012, 30(3): 257-261. (in Chinese with English abstract)

[18] Pouteau R, Rambal S, Ratte J P, et al. Downscaling MODIS-derived maps using GIS and boosted regression trees: The case of frost occurrence over the arid Andean highlands of Bolivia[J]. Remote Sensing of Environment, 2011, 115(1): 117-129.

[19] Kurowicka D, Cooke R M. Uncertainty Analysis with High Dimensional Dependence Modelling[M]. Chichester: John Wiley & Sons, 2006.

[20] 郭愛軍，暢建霞，王義民，等. 近50年涇河流域降雨-徑流關系變化及驅動因素定量分析[J].農業工程學報，2015，31(14)：165-171.

Guo Aijun, Chang Jianxia, Wang Yimin, et al. Variation characteristics of rainfall-runoff relationship and driving factors analysis in Jinghe river basin in nearly 50 years [J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(14): 165-171. (in Chinese with English abstract)

[21] 王曉峰，張園，馮曉明，等. 基于游程理論和Copula函數的干旱特征分析及應用[J]. 農業工程學報，2017，33(10)：206-214.

Wang Xiaofeng, Zhang Yuan, Feng Xiaoming, et al. Analysis and application of drought characteristics based on run theory and Copula function[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(10): 206-214. (in Chinese with English abstract)

[22] 藺文慧，宋松柏. 基于混合Copula函數的洪水聯合概率計算方法[J]. 水力發電學報，2021，40(12)：40-51.

Lin Wenhui, Song Songbai. Study on calculation method of flood joint probability based on mixed Copula function[J]. Journal of Hydroelectric Engineering, 2021, 40(12): 40-51. (in Chinese with English abstract)

[23] 鄔定榮，霍治國，王培娟，等. 陜西蘋果花期機理性預報模型的適用性評價[J]. 應用氣象學報，2019，30(5)：555-564.

Wu Dingrong, Huo Zhiguo, Wang Peijuan, et al. The applicability of mechanism phenology model to simulating apple flowering date in Shaanxi provice[J]. Journal of Applied Meteorological Science, 2019, 30(5): 555-564. (in Chinese with English abstract)

[24] 王潤紅，茹曉雅，蔣騰聰，等. 基于物候模型研究未來氣候情景下陜西蘋果花期的可能變化[J]. 中國農業氣象，2021，42(9)：729-745.

Wang Runhong, Ru Xiaoya, Jiang Tengcong, et al. Based on the phenological model to study the possible changes of apple flowering dates under future climate scenarios in Shaanxi province[J]. Chinese Journal of Agrometeorology, 2021, 42(9): 729-745. (in Chinese with English abstract)

[25] 劉映寧，賀文麗，李艷莉，等. 陜西果區蘋果花期凍害農業保險風險指數的設計[J]. 中國農業氣象，2010，31(1)：125-129，136.

Liu Yingning, He Wenli, Li Yanli, et al. Chinese Journal of Agrometeorolog[J]. A study on the risk index design of agricultural insurance on apple florescence freezing injury in Shaanxi fruit zone, 2010, 31(1): 125-129,136. (in Chinese with English abstract)

[26] 李美榮，朱琳，杜繼穩. 陜西蘋果花期霜凍災害分析[J]. 果樹學報，2008(5)：666-670.

Li Meirong, Zhu Lin, Du Jiwen. Analysis of frost damage during apple bloom period in Shaanxi Province[J]. Journal of Natural Disasters, 2008(5): 666-670. (in Chinese with English abstract)

[27] 王景紅. 陜西主要果樹氣候適宜性與氣象災害風險區劃圖集[M]. 西安：陜西科學技術出版社，2012：85-95.

[28] 韓會明，劉喆玥，劉成林，等. 基于Copula函數的贛江流域氣象干旱特征分析[J]. 水電能源科學，2020，38(8)：9-13.

Han Huiming, Liu Zheyue, Liu Chenglin, et al. Analysis of meteorological drought characteristics in Ganjiang River basin based on Copula function[J]. Water Resources and Power, 2020, 38(8): 9-13. (in Chinese with English abstract)

[29] 張卓群，馮冬發，侯宇恒. 基于Copula函數的黃河流域干旱特征研究[J]. 干旱區資源與環境，2022，36(1)：66-72.

Zhang Zhuoqun, Feng Dongfa, Hou Yuheng. Study on drought characteristics in Yellow River basin based on Copula function[J]. Journal of Arid Land Resources and Environment, 2022, 36(1): 66-72. (in Chinese with English abstract)

[30] Michele C D, Salvadori G, Vezzoli R, et al., Multivariate assessment of droughts: Frequency analysis and dynamic return period[J]. Water Resources Research, 2013, 49(10): 6985-6994.

[31] Sklar A, Fonctions de repartition an dimensions et leurs marges[J]. Publ. inst. statist. univ. Paris, 1959, 8: 229-231.

[32] Salvadori G. Bivariate return periods via 2-Copulas[J]. Statistical Methodology, 2004, 1(1): 129-144.

[33] 張曉煜，馬玉平，蘇占勝，等. 寧夏主要作物霜凍試驗研究[J]. 干旱區資源與環境，2001(2)：50-54.

Zhang Xiaoyu, Ma Yuping, Su Zhansheng, et al. Experiments on frost injuries of main crops in Ningxia Province[J]. Journal of Arid Land Resources and Environment, 2001(2): 50-54. (in Chinese with English abstract)

[34] 王景紅，柏秦鳳，梁軼，等. 2013年陜西蘋果花期凍害氣象條件分析及受凍指標研究[J]. 果樹學報，2015，32(1)：100-107，174.

Wang Jinghong, Bai Qinfeng, Liang Yi, et al. Study on freezing index and weather conditions causing Shaanxi apple florescence freezing injury in 2013[J]. Journal of Natural Disasters, 2015, 32(1): 100-107, 174. (in Chinese with English abstract)

Analysis of the characteristics of apple later frost risks in Shaanxi Province based on Copula functions

Jiang Yuan1,2, Ru Xiaoya1,2, Luo Qi1,2, Li Meirong3, Wang Zhao3, Wang Jinghong3, Feng Hao2,4, Zhang Dong5, Su Baofeng6, Yu Qiang4, He Jianqiang1,2※

(1.712100,; 2.712100,; 3.710015,; 4.,712100,; 5.,712100,; 6.,,712100,)

Late frost is one of the most destructive meteorological disasters in the Loess Plateau of China. A great threat can be posed to the sustainable production of apples, leading to great economic losses in the apple industry. Thus, it is a high demand to explore the occurrence of late frost events for the prevention of apple late frost disasters in the decision-making of the local apple industry. In this study, the late frost return periods were investigated on the duration and severity of late frost events using the Coupla functions. The reliability of the model was then verified to analyze the characteristics of apple late frost. The study area was taken as the apple-producing areas of Shaanxi Province in western China. The meteorological datasets were collected from the seven weather stations from 1971-2018. The daily minimum temperature (min) of 0℃ was taken as the critical temperature for the occurrence of apple late frost events, in order to extract the two characteristic variables of duration and severity of late frost events. These characteristic variables of late frost events were first fitted by seven common distribution functions, respectively. Kolmogorov Smirnov (K-S) test was then carried out to verify the model. The joint distributions of late frost characteristic variables were constructed to evaluate the goodness-of-fit using six Copula functions. The occurrence probability and return period of late frost events were analyzed with the optimized Copula functions. The results showed that the severity of late frost risks generally increased from the southeast to northwest in the study area from 1971 to 2018. The optimal marginal distribution of late frost duration was in the log-normal distribution, while there was a great difference in the optimal marginal distributions of late frost severity. A significant positive correlation was found between the duration variables and the severity of late frost at each station. The joint cumulative probability increased significantly, as the severity and duration of late frost increased, but the increasing trend was much slower than before. A much more significant increase was observed in the co-occurrence return period under the same increase of univariate value, compared with the joint. The univariate return period was always between the joint and co-occurrence return period. Once the univariate return period was small enough, the optimal range of the actual univariate return period can be estimated, according to the joint and the co-occurrence return period. In general, a low probability was found in the late frost events with long duration and high severity at the weather stations in apple-producing areas in Shaanxi Province of China. However, the stations in the Yan’an area were more susceptible to the late frost events with high severity or long duration, as well as the late frost events with both high severity and long duration. Thus, more attention can be paid to late frost risks in the Yan’an area. This finding can provide a theoretical base to deal with the late frost disaster in apple production.

meteorological; disasters; late frost; Copula function; frequency analysis; return period; apple

10.11975/j.issn.1002-6819.2022.23.018

S166

A

1002-6819(2022)-23-0170-11

姜元，茹曉雅，羅琦，等. 基于Copula函數的陜西蘋果晚霜凍特征分析[J]. 農業工程學報，2022，38(23)：170-180.doi：10.11975/j.issn.1002-6819.2022.23.018 http://www.tcsae.org

Jiang Yuan, Ru Xiaoya, Luo Qi, et al. Analysis of the characteristics of apple later frost risks in Shaanxi Province based on Copula functions[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(23): 170-180. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2022.23.018 http://www.tcsae.org

2022-08-21

2022-11-26

國家重點研發計劃（2021YFD1900700）；陜西省重點研發計劃重點產業創新鏈（群）-農業領域項目（2019ZDLNY07-03）；陜西省氣象局秦嶺和黃土高原生態環境氣象重點實驗室2019年開放研究基金課題（2019Z-5）；西北農林科技大學人才專項資金（千人計劃項目）；高等學校學科創新引智計劃（111計劃）(B12007)

姜元，博士生，研究方向為農業生態系統模擬。Email：yuanjiang59@nwafu.edu.cn

何建強，教授，博士生導師，研究方向為農業生態系統模擬。Email：jianqiang_he@nwafu.edu.cn