湖南油茶花期低溫陰雨災損評估模型構建

張 超，謝佰承

（湖南省氣象科學研究所 氣象防災減災湖南省重點實驗室，湖南 長沙 410118）

油茶Camellia oleifera是我國特有的重要木本油料樹種[1]。目前，我國油茶主要分布在長江流域及其以南的14個省份，其中，湖南省的油茶種植面積、茶油產量及其產值均居全國首位[2]。同所有露天生產的農作物一樣，氣象條件是影響油茶產量、品質、成本的重要因素[3-8]。油茶的營養生長和生殖生長互相影響、交錯，從春梢生長到果實成熟每個物候期的氣象條件均會影響油茶的生長發育，進而影響到油茶產量。

油茶是我國獨有的樹種，因此關于油茶生長與氣象條件的關聯性研究主要集中在國內。20世紀80年代后，許多學者相繼開展了油茶種植的氣候適應性研究[3,9-15]，找出了利于油茶生長的氣象條件以及影響油茶生長的主要氣象災害[5,16-18]，通過分析氣象條件與油茶產量的關系[8,19-21]，建立了油茶氣象產量預測模型[4,22-24]，并探討了氣象條件與油茶經濟性狀、含油率等的關系[6-7,25-31]，但關于油茶氣象災害損失評估技術的相關研究鮮見報道，因此在前人關于氣象災害對油茶產量的影響研究成果的基礎上，開展油茶氣象災害災損評估模型的研制，對農業氣象防災減災具有重要作用，可提升氣象為農服務水平。

遍歷油茶各物候期，油茶開花授粉時的氣象條件是影響油茶產量的主要因素之一，而對油茶開花授粉影響最大的是花期低溫陰雨[9]。油茶花期一般在10月上旬—12月下旬，是影響油茶產量的關鍵物候期，花期若氣溫低、降雨多、日照少，會造成油茶大量落花、坐果率低[16]。當最低氣溫為0 ℃時，油茶花瓣開始有輕度凍害[19]；盛花期出現小于-2.5 ℃的低溫，會凍壞花蕾，凍死傳粉地蜂，造成次年減產2～4成[20,32-33]，若最低氣溫更低且持續時間延長，減產幅度一般在3成以上，甚至可達9成[33]。油茶花期遇陰雨天或出現低于0 ℃的低溫天氣，坐果率極低，甚至為零[34]。連陰雨可使整個花期推遲而導致減產，且對油茶產量的影響隨著盛花期內降雨時間和降雨總量的增加而加大[35]。在油茶盛花期：雨日大于13 d會影響豐產[3]；雨日超過50%，對授粉結實不利，產量降低；雨日超過70%時，較難授粉結實，產量極低[36]。

目前，花期低溫陰雨災害發生頻率有增大趨勢[37]，使得油茶減產風險進一步增大，而花期對油茶產量的影響最為顯著[22]。鑒于此，本研究中以花期低溫陰雨災害為研究對象，以判別是否會造成災損為目標，基于篩選出的花期低溫陰雨災害的關鍵氣候影響因子，采用極端梯度提升算法來構建湖南油茶花期低溫陰雨災損評估概率模型，旨在為湖南油茶氣象災害監測和評估業務提供參考。

1 資料與方法

1.1 資料來源

1.1.1 油茶測產資料

油茶測產資料為湖南省林業科學院提供的2006—2016年湖南油茶種植區的31塊油茶樣地（圖1）的油茶鮮果單位面積產量數據，該數據由湖南省林業科學院技術人員選擇各油茶樣地代表性地塊采摘鮮果后測定。不同樣地的油茶種植年份、品種及土壤等均有差異，這些因素均會對油茶產量造成影響，因此在使用油茶測產資料前對其進行處理。

圖1 湖南省31個油茶樣地分布Fig. 1 Distribution of 31 C. oleifera sample plots of Hunan

油茶良種林（指經國家林業局林木品種審定委員會審定的油茶良種）和非良種林（指其他未審定的品種）產量差異較大，因此從31塊油茶樣地中剔除11塊非良種林油茶樣地的產量數據；油茶從種植到進入豐產期需7～8 a的時間，在進入豐產期前，隨著油茶樹不斷生長，產量遞增，因此剔除未進入豐產期的油茶產量數據。不同良種林之間油茶產量的差異較大，為減少不同區域、不同時期生產力水平對油茶產量的影響，采用相對氣象產量（yr）來表征油茶產量，相對氣象產量為正值表示增產，負值表示減產，計算公式為yr= [ (ya-yt) /yt]×100%。式中：ya為實際產量，yt為趨勢產量。利用5點滑動平均法計算趨勢產量yt，趨勢產量主要反映農業技術水平的提高對產量的影響，具有漸進性和相對穩定性。

經上述處理，共得到83組油茶樣地相對氣象產量資料，其中增產資料41組，減產資料42組。以相對氣象產量為指標，分離了油茶種植年份、品種及土壤對油茶產量資料的干擾，可認為油茶產量僅受氣象因素的影響。

1.1.2 氣象資料

氣象資料來源于湖南省信息中心，包括湖南97個氣象站和自動觀測站2006—2016年油茶花期逐日平均氣溫、最高氣溫、最低氣溫、降水量和日照時長資料。為了使獲取的油茶樣地氣象要素數據更加客觀真實，采用氣象要素小網格推算技術，將湖南各站點數據插值為500 m×500 m分辨率的格點數據，選取離油茶樣地最近的格點值作為各油茶樣地的站點值。

1.2 研究方法

1.2.1 極端梯度提升

極端梯度提升（extreme gradient boosting，XGBoost）屬監督模型，是應用機器學習領域一個強有力的工具，具有防止過擬合、算法效率高、損失計算精確等特點。作為Boosting集成學習算法家族中的一員，XGBoost屬于樹集成模型，將K個CART回歸樹的結果進行求和，作為最終的預測值。

XGBoost相對于線性模型在進行預測時往往有更好的精度，但是也失去了線性模型的可解釋性，通常被認為是黑箱模型。2017年，Lundberg等[38]提出了SHAP（shapley additive explanation）值這一廣泛適用的方法，可用來解釋XGBoost等各種模型。SHAP屬于模型事后解釋的方法，其核心是計算其中每個特征（即氣候影響因子）的SHAP值。對于每個預測樣本，模型均產生1個預測值，SHAP值即該樣本中每個特征所分配到的數值。假設第i個樣本為xi，第i個樣本的第j個特征為xij，模型對第i個樣本的預測值為yi，整個模型的基線（通常是所有樣本的目標變量的均值）為yb，那么SHAP值服從等式：

式中：f(xij)為xij的SHAP值。直觀上看，f(xi1)就是第i個樣本中第1個特征對最終預測值yi的貢獻值，當f(xi1)＞0，說明該特征提升了預測值，有正向作用；反之，說明該特征使得預測值降低，有反作用。SHAP值類似回歸系數，有正負之分，有大小之分。

根據計算得到的各影響因子SHAP值的大小，判斷各影響因子在建模中的貢獻大小，從而實現對關鍵影響因子的篩選，并利用篩選出的關鍵影響因子來構建油茶花期低溫陰雨災損評估概率模型。

1.2.2 XGBoost模型檢驗

采用受試者工作特征曲線（receiver operating characteristic，ROC）對構建的概率模型進行檢驗。以曲線下面積（area under the curve，AUC）的值來判斷模型預測的準確度，若AUC值大于0.5，且越接近于1，說明模型的擬合效果越好。

XGBoost模型中，判斷閾值（臨界值）是評判模型預測概率的關鍵。通常情況下，以系統默認值0.5作為評判標準，但這樣有可能造成較大偏差，因此采用ROC中的敏感度和特異性值計算約登指數，從而得到最佳臨界值，進而對概率進行分類。如果模型的預測概率大于該最佳臨界值，則認為有災損發生，小于該最佳臨界值，則認為無災損發生。

2 結果與分析

2.1 油茶花期低溫陰雨災害關鍵影響因子

構建油茶花期災損模型時，首先要對花期內不利于油茶生長的氣候因子實行全覆蓋，然后從中找出關鍵氣候影響因子?；ㄆ诘闹饕獨庀鬄暮榈蜏仃幱?，通過查閱相關文獻資料、結合前期研究成果及借鑒專家經驗，找出了16個與花期低溫陰雨相關聯、能反映災害程度、表征災害影響的因子（表1）。根據計算得到的各影響因子SHAP值的大小，判斷各影響因子在建模中的貢獻大小，并結合相關分析和方差分析實現對關鍵影響因子的篩選。

表1 油茶花期低溫陰雨的主要氣候影響因子以及方差分析、相關分析結果?Table 1 The main climatic influencing factors of low temperature and overcast rain damage at florescence of C. oleifera,results of variance analysis and correlation analysis

2.1.1 主要氣候影響因子的相關分析和方差分析

將油茶樣地的相對氣象產量與氣候影響因子做相關分析，結果見表1。從表1可以看出：僅日平均氣溫不小于10 ℃活動積溫、日平均氣溫小于10 ℃積溫、平均溫度3個因子與油茶的相對氣象產量具有較好的相關性，通過α=0.05的顯著性檢驗；平均最高溫度、平均最低溫度2個因子僅通過α=0.1的顯著性檢驗；其他因子與油茶的相對氣象產量的相關性不明顯。

將油茶樣地增產和減產2類樣本所對應的主要氣候影響因子分別做方差分析，找出2類樣本間具有顯著差異的氣候因子，結果見表1。從表1可以看出，日平均氣溫不小于10 ℃活動積溫、日平均氣溫小于10 ℃積溫、平均溫度、平均最高溫度、平均最低溫度、有日照時間（d）這6個因子的P值均小于0.05，累積日照時間（h）、無日照時間（d）和寡照時間（d）這3個因子的P值均小于0.1，說明在油茶增產與減產的年份間這幾個因子具有較顯著的差異，其他因子的差異性不明顯。

2.1.2 關鍵影響因子的篩選

采用XGBoost算法計算16個因子的SHAP值，16個因子在各樣本中SHAP值平均絕對值的排序如圖2所示，這可以看作因子重要性的排序。從圖2可以直觀看出，平均最低氣溫對建模至關重要，該因子所有樣本的SHAP值的平均絕對值最大，以SHAP值的平均絕對值大于0.5作為篩選因子的標準，選取平均最低氣溫、日平均氣溫小于10 ℃積溫、寡照時間（d）和極端最低氣溫為主要的影響因子。

圖2 各油茶樣地樣本中16個因子SHAP值的平均絕對值排序Fig. 2 Ranking of average absolute values of 16 factor SHAP values in each C. oleifera sample plot

以SHAP值分析結果為基礎，結合方差和相關分析的結果，選取相關性較好、方差檢驗顯著且對建模較重要的因子作為花期低溫陰雨的關鍵影響因子。在按SHAP值排序選取的4個因子中，平均最低氣溫和日平均氣溫小于10 ℃積溫與樣本相關性好，且增產與減產樣本的差異顯著；寡照時間（d）雖無明顯的相關性，但其增產與減產樣本的差異通過了顯著性檢驗；無論是相關性分析還是方差檢驗，極端最低氣溫均未通過顯著性檢驗，但根據前期的研究成果，0 ℃以下低溫的出現會對油茶產量造成較大的影響。經綜合分析，最終選取平均最低氣溫、日平均氣溫小于10 ℃積溫、寡照時間（d）和極端最低氣溫4個因子來構建花期低溫陰雨的災損模型。

2.2 油茶花期低溫陰雨災損模型

2.2.1 災損模型的訓練與檢驗

基于篩選出的4個花期低溫陰雨關鍵氣候影響因子，采用XGBoost算法來構建油茶花期低溫陰雨災損概率模型。選取總樣本的80%來訓練模型，用剩下的20%的樣本來進行模型驗證。XGBoost算法包括正則化系數、學習率、n估計量和樹的最大深度等眾多超參數。超參數定義了模型的復雜度或學習能力等特定基本屬性，是在開始學習過程之前需要確定的參數，調節超參數的意義在于最小化期望風險，使模型的優化度與復雜度達到平衡，盡可能同時避免欠擬合和過擬合[39]。

網格搜索是應用最廣泛的建立在交叉驗證基礎上的超參數搜索算法，這種窮舉式調參算法通過循環遍歷嘗試每種參數組合的可能性，找出表現最好的組合，找到全局最大或最小值。在此過程中采用了5折交叉驗證，即將訓練集5等分，取其中1份為驗證集，其余4份為新訓練集，經過5次在不同驗證集上的測試，取最優結果所對應的超參數組合。經驗證，得到的模型預測效果最好的參數組合見表2，準確率最高，為82.4%。

表2 XGBoost最優超參數配置Table 2 XGboost optimal hyperparameter configuration

同時采用ROC曲線和AUC值對模型進行評估，結果如圖3所示。圖3中實線為ROC曲線，該曲線越靠攏（0，1）點、越偏離45°對角線（虛線）表示模型預測效果越好。AUC是ROC曲線下的面積，AUC越接近1表明模型預測效果越好。所構建的油茶災損概率模型的ROC曲線偏離45°對角線較遠，較靠近（0，1）點，AUC值達到了0.82，說明該模型的預測效果較好。

圖3 油茶花期低溫陰雨災損概率模型的ROC曲線和AUC值Fig. 3 ROC curve and AUC value of the loss estimation model of low temperature and overcast rain damage at florescence to C. oleifera

根據ROC中的敏感度和特異性值計算約登指數，從而得到最佳臨界值（圖3中的圓點）為0.8。當模型預測的災損概率值超過0.8時認為有災損發生，當概率值低于0.8時則認為無災損發生。

2.2.2 關鍵影響因子對建模的影響

計算4個關鍵因子的SHAP值，分析各因子在構建災損模型時貢獻的大小，結果如圖4所示。從整體影響來看，極端最低氣溫對模型預測結果的貢獻最大，其次是平均最低氣溫、日平均氣溫小于10 ℃積溫、寡照時間。

各油茶樣地樣本中4個關鍵影響因子SHAP值的分布如圖5所示。由圖5可以看出：當極端最低氣溫小于0 ℃時，SHAP值多為負值，對預測災損發生有貢獻；平均最低氣溫與油茶產量有一定的正相關趨勢，當平均最低氣溫小于9 ℃時SHAP值多為負值，對預測災損發生有貢獻；日平均氣溫小于10 ℃積溫與油茶產量有一定的負相關趨勢，當日平均氣溫小于10 ℃積溫超過150 ℃時，SHAP值多為負值，對預測災損發生有貢獻；寡照時間與油茶產量也有一定的負相關趨勢，當寡照時間超過25 d時，SHAP值多為負值，對預測災損發生有貢獻。

圖4 各油茶樣地樣本中4個關鍵影響因子SHAP值的平均絕對值排序Fig. 4 Ranking of average absolute values of 4 key factor SHAP values in each C. oleifera sample plot

圖5 各油茶樣地樣本中4個關鍵影響因子SHAP值的分布Fig. 5 Distribution of 4 key factor SHAP values in each C. oleifera sample plot

2.2.3 災損模型的應用

將油茶花期低溫陰雨災損模型應用于2020年油茶花期氣象條件評估，結果如圖6所示。由圖6可以看出：根據最佳臨界值0.8，湖南全省大部分地方無災損發生，僅湘西北、湘西南的部分地區和湘東北局地有發生災損的風險，與花期氣象條件分析結果（湘西、湘北和湘南有發生低溫陰雨災害的風險）基本吻合，且評估結果較客觀，該模型可用于油茶花期低溫陰雨災損的評估。

3 結論與討論

3.1 結 論

在湖南油茶花期低溫陰雨災損評估模型的構建過程中，相關分析、方差分析和SHAP值分析結果表明，在16個油茶花期低溫陰雨災害氣候影響因子中，極端最低氣溫、平均最低氣溫、日平均氣溫小于10 ℃積溫和寡照時間這4個因子與油茶產量相關性較好，在增產與減產樣本間的差異顯著，在模型構建中的貢獻較大，是花期低溫陰雨災害的關鍵影響因子。

通過進一步分析這4個關鍵影響因子對油茶災損模型預測值貢獻的大小，得到了可判別災損發生的臨界值。當日平均氣溫小于10 ℃積溫不小于150 ℃，或寡照時間不小于25 d，或平均最低氣溫不大于9 ℃時，往往會伴有低溫陰雨災害出現，這些指標對災損的發生具有一定的指示意義，在實際工作中可用于花期低溫陰雨災害的監測，也可對花期是否有災損發生進行初步預判。

使用建模時預留的20%的樣本，對模型的預測準確率進行驗證，結果表明基于XGBoost算法建立的油茶花期低溫陰雨災損模型具有較好的預測效果。該模型的預測準確率達到了82.4%，ROC曲線偏離45°對角線較遠，比較靠近（0，1）點，AUC值高達0.82，這些評估指標均說明該模型具有較好的預測能力。

圖6 2020年湖南省油茶花期災損發生概率分布Fig. 6 Probability distribution of Hunan disaster loss at florescence of C. oleifera in 2020

將該模型應用于2020年油茶花期低溫陰雨的災損評估，模型的評估結果與花期氣象條件綜合分析的結果基本吻合，說明該模型具有較好的適用性，可用于油茶花期低溫陰雨災損評估業務工作中。

3.2 討 論

油茶花期是影響油茶產量的關鍵物候期，基于花期氣象條件建立的油茶產量模型，可用于開展油茶產量預估[22]，同樣，基于氣象條件建立的花期災損評估模型可為油茶災損進行預估，由于花期具有較長的時間提前量，可為后期采取防災減災措施提供依據。在篩選油茶花期關鍵影響因子時，首先通過查閱相關文獻資料、結合前期研究成果以及借鑒專家經驗，對花期內不利于油茶生長的氣候因子實行全覆蓋，然后對所有影響因子進行了系統分析，剔除了對油茶影響相對較小的因子，最終選取極端最低氣溫、平均最低氣溫、日平均氣溫小于10 ℃積溫和寡照時間來構建模型，這些均是反映低溫、陰雨程度的因子，與前人研究中[19,32-34]所提到的0 ℃以下的低溫、連陰雨天氣等造成油茶減產的重要因素具有一致性。蔣元華等[22]在利用油茶樣地測產數據構建花期油茶產量逐步回歸預測模型時，排在前10位的關鍵氣象因子主要有連續有日照時間、最低溫度小于0 ℃時間、平均最低氣溫、有日照時間、雨日時間等。本研究中在分析因子的重要性時，將最低溫度小于0 ℃時間、有日照時間、雨日時間等對油茶花期影響不明顯的因子進行了剔除，這種客觀與主觀相結合的方法，比利用逐步線性回歸自動篩選因子的方法更有針對性，更準確。

由于油茶樣地測產樣本數量有限，在構建油茶花期低溫陰雨災損模型時，所用到的樣本為油茶測產全體樣本，并未區分哪些減產樣本主要是由花期低溫陰雨災害造成的，減產樣本與花期低溫陰雨災害未一一對應起來，因此，在建模時可能會出現選取的關鍵影響因子不夠準確，致使關鍵影響因子對某些樣本減產貢獻不明顯，甚至對增產有貢獻的現象，對模型的預測準確性造成一定影響，可通過增加樣本并對樣本按災害分類的方式來進一步提升模型預測的準確率。

雖然花期的氣象條件對油茶產量的影響最為顯著，模型的預測結果可用于油茶災損的預估，但花期并不能代表整個生育期，其預估結果與油茶的實際災損情況會有一定偏差，為提升預測結果的準確性，后期可采用構建油茶花期災損模型的方法，篩選全生育期的關鍵影響因子，構建油茶全生育期的災損模型。