张 超，谢佰承
（湖南省气象科学研究所 气象防灾减灾湖南省重点实验室，湖南 长沙 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 油茶测产资料
油茶测产资料为湖南省林业科学院提供的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.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.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 ℃时间、有日照时间、雨日时间等对油茶花期影响不明显的因子进行了剔除，这种客观与主观相结合的方法，比利用逐步线性回归自动筛选因子的方法更有针对性，更准确。
由于油茶样地测产样本数量有限，在构建油茶花期低温阴雨灾损模型时，所用到的样本为油茶测产全体样本，并未区分哪些减产样本主要是由花期低温阴雨灾害造成的，减产样本与花期低温阴雨灾害未一一对应起来，因此，在建模时可能会出现选取的关键影响因子不够准确，致使关键影响因子对某些样本减产贡献不明显，甚至对增产有贡献的现象，对模型的预测准确性造成一定影响，可通过增加样本并对样本按灾害分类的方式来进一步提升模型预测的准确率。
虽然花期的气象条件对油茶产量的影响最为显著，模型的预测结果可用于油茶灾损的预估，但花期并不能代表整个生育期，其预估结果与油茶的实际灾损情况会有一定偏差，为提升预测结果的准确性，后期可采用构建油茶花期灾损模型的方法，筛选全生育期的关键影响因子，构建油茶全生育期的灾损模型。
猜你喜欢阴雨样地油茶仁怀市二茬红缨子高粱的生物量及载畜量调查贵州畜牧兽医(2022年3期)2022-06-28油茶种植市场展望及营造林技术探究农村科学实验(2021年7期)2021-12-24扶绥县春耕春播期低温阴雨的变化规律和发展趋势广西农学报(2021年4期)2021-12-20额尔古纳市兴安落叶松中龄林植被碳储量研究现代园艺(2021年23期)2021-12-01持续阴雨期间鳜鱼养殖池常见问题及管理措施当代水产(2021年6期)2021-08-13垦荒东坡赤壁诗词(2021年1期)2021-03-24基于角尺度模型的林业样地空间结构分析新农业(2020年18期)2021-01-0715 年生鹅掌楸林分生长差异性研究花卉(2020年24期)2021-01-04家乡的油茶故事作文·低年级(2020年7期)2020-07-28中国油茶花开的春天中国林业产业(2019年11期)2019-01-09
