程晋昕 段长春 闫生杰
1)(云南省气候中心, 昆明 650034)2)(云南省气象科学研究所, 昆明 650034)3)(云南省大理白族自治州祥云县气象局, 大理 672100)
薄壳山核桃又名美国山核桃,原产于美国和墨西哥北部,主要种植区位于美国佐治亚、新墨西哥、露易斯安那、得克萨斯等州[1]。云南省1974年首次将薄壳山核桃作为果用树种引种[2],并成为全国最早开展薄壳山核桃产业化种植的省份,也是种植面积和产量最大的省份[3]。近年薄壳山核桃作为一种兼具经济、社会和生态效益的特色作物,其种植产业受到很大关注[4],开展气候适宜区划对云南省薄壳山核桃种植产业的规划、扩大与调整具有指导意义。
针对薄壳山核桃种植的气候适宜性,张日清等[5-7]结合薄壳山核桃原产地气候生态位宽度与国内引种地的前期引种效果,对全国范围内的种植适宜性进行了区划研究;Sparks[8]从物种环境适应性的角度,分析薄壳山核桃对气候条件的需求;董润泉[9]对云南省薄壳山核桃的引种和发展进行总结,对不同气候带的种植效果进行简要分析。上述研究多依赖于专家经验和站点观测,主观性较强、精细化程度低,应用于地形与气候条件复杂多样的云南省具有一定局限性。近年来,GIS与机器学习技术在农业气象领域得到广泛应用[10-13]。其中,最大熵模型(MaxEnt)因其客观、定量的特点及良好的性能成为农业气候区划的重要工具,并成功应用于玉米、苹果、猕猴桃、当归等作物的气候适宜性和风险区划[14-17]。然而,上述研究需要研究区域具备数量较多、代表性强的种植点作为数据支撑,云南省薄壳山核桃种植产业规模有限,种植点空间分布数据相对稀缺。
本研究以最大熵模型和GIS技术为基础,提出一种农业气候适宜性区划方法。利用美国本土(contiguous United States)的薄壳山核桃种植点和气候数据,实现主导气候因子的筛选和MaxEnt模型训练,在改进的基础上构建气候区划模型,基于中国云南省气候数据开展薄壳山核桃种植气候适宜性区划,为种植产业的规划和发展提供参考。
本研究所用数据主要包括气候数据、薄壳山核桃种植点分布数据、地形数据和地理信息数据等。其中,气候数据包括年平均气温、30年极端最低气温、无霜期、活动积温(指日平均气温稳定通过10.0℃ 的积温,下同)、月平均气温、年降水量、年日照时数,涉及中国云南省和美国本土两个区域。其中,中国云南省气候数据采用1981—2010年全省125个国家级气象站累年统计结果,来源于云南省气象信息中心的全国综合气象信息共享平台(CIMISS)[18]。
美国本土气候数据主要用于薄壳山核桃原生境气候条件分析和区划模型的构建,包括两个数据来源:年平均气温、年降水量、月平均气温数据来源于世界气候网(WORLDCLIM,http:∥www.worldclim.org)提供的空间网格数据集[19],该数据集在农业气候适应性区划、物种分布预测和气候变化等研究领域得到广泛认可[20-22]。由于WORLDCLIM未提供无霜期、日照时数和活动积温数据,故采用美国气象站点逐日数据进行累年值统计,数据来源于美国国家海洋和大气管理局国家气候数据中心(NOAA National Climatic Data Center)提供的全球历史气候网络日数据集(GHCN-DAILY)[23]。为保证数据的可靠性,本研究仅保留数据记录完整的世界气象组织站点,包含176个气温观测站点和138个日照观测站点,部分气温站点与日照站点重叠。
薄壳山核桃种植点分布数据用于种植区气候要素提取,来源于全球生物多样性信息数据库网络(GBIF,http:∥www.gbif.org)提供的薄壳山核桃物种采样信息和谷歌地图提供的薄壳山核桃种植点标记信息。为避免局部区域分布点过于密集影响最大熵模型效果,采用0.5°×0.5°的网格对美国本土进行空间划分,每个网格中仅保留1个样本点,最终得到有效样本274个。地形数据用于气候要素的插值,采用空间分辨率为30″(约1 km2)的航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission,SRTM)数字高程模型产品(DEM)[24]提供海拔高度信息,并利用ArcGIS软件提取坡度和坡向。地理信息数据中,中国云南省采用1:25万中国空间矢量数据集,美国本土则采用1:100万全球空间矢量数据集,两个数据集均来源于中国云南省气象信息中心。
1.2.1 潜在气候因子的选取
热量条件是影响薄壳山核桃种植的重要因素。张日清等[5]指出,薄壳山核桃原生境中心产区年平均气温为13.2~20.3℃,夏季平均气温为24.0~30.0℃,7月平均气温为25.6~29.4℃,无霜期为180.0~280.0 d,生长季中大于等于10℃积温为2500~3700℃·d,热量条件的不足易导致结实不良。董润泉[9]总结云南省各地产的坚果质量,得到类似结论:暖温带所产坚果不具备食用经济价值,北亚热带、中亚热带、南亚热带坚果品质依次递增。低温也是影响薄壳山核桃的春季发芽的重要条件,需冷量因品种而异[5],但过度低温意味着生长季过短,进而对结实造成不利影响[8]。极端最低气温对薄壳山核桃的分布也有一定限制,不同品种差异较大,一般南方品种不能低于-18.0℃[5],部分北方品种能够承受-40.0~-34.0℃的极端低温[8]。薄壳山核桃的整个生长发育周期需要充足的水分条件,张日清等[5]提出,年降水量为1000.0~1600.0 mm的地区较为理想。在美国西部的干旱和半干旱地区,依靠灌溉和地下水的补充,年降水量下限可低至600.0 mm[8]。在降水时间分布方面,Sparks[8]指出,薄壳山核桃在芽膨大、果实膨大和果仁发育3个阶段对水分胁迫较为敏感,分别与3—5月、6—7月和8—9月3个时期的降水条件密切相关。在光照条件方面,薄壳山核桃属于喜光作物[25],对全年日照有较高要求。此外,4—5月正值花期到幼果成熟期,日照充足、温暖少雨有利于坐果率提高[26]。结合已有研究成果,本文选取年平均气温、活动积温、无霜期、30年极端最低气温、1月平均气温、7月平均气温、年降水量、3—5月降水量、6—7月降水量、8—9月降水量、年日照时数和4—5月日照时数12个要素作为MaxEnt模型构建的潜在气候因子,开展气候适宜性分析。
1.2.2 气候因子插值
考虑到云南省气候和地形条件复杂,本研究以站点数据和DEM数据为基础,采用1 km水平分辨率进行气候因子的插值。对于平均气温、活动积温、无霜期、30年极端最低气温、1月平均气温、7月平均气温,采用趋势面法进行插值,其基本思路为气候要素与当地的经度、纬度、海拔等地理因子关系密切[27],可表示为
y=F(λ,φ,β1,β2,β3,…)+ε。
(1)
式(1)中,y表示气候要素值;函数F体现地形的趋势特征,在样本较少的情况下通常采用线性函数,可通过站点数据进行拟合;λ和φ分别表示经度和纬度;β1,β2和β3表示地形因子,如海拔高度、坡度和坡向;ε为残差项,体现局部尺度的综合影响,可采用反距离权重法插值。
趋势面法对于降水量和日照时数推算适用性较差[28],故采用薄盘光滑样条函数法对年降水量、3—5月降水量、6—7月降水量、8—9月降水量、年日照时数和4—5月日照时数进行空间插值,地形因子的影响则通过协变量引入。为减小MaxEnt模型采样偏差,基于薄壳山核桃种植点分布数据空间范围,在完成插值后,按照24.0°~42.0°N,70.0°~125.0°W的范围对美国本土气候要素网格数据进行空间裁切。
1.2.3 基于MaxEnt模型的气候适宜性区划
最大熵原理[29]是统计学和机器学习领域的一种经典理论,在所有满足已知约束条件的概率模型中,熵最大模型与真实概率分布最接近。以该理论为基础,Phillips等[30]和Elith等[31]提出用于物种分布预测的最大熵模型(MaxEnt),利用已知物种分布数据和环境因子的空间分布,对物种存在概率的空间分布进行估计。物种的存在概率越高,代表环境因子的适宜性越强。在农业气候适宜性区划应用中,一般采用研究区域内若干个种植区样本点完成MaxEnt模型的训练,将覆盖该区域的气候数据输入模型,对模型输出结果划分适宜性等级[14,16-17]。
本文采用美国自然历史博物馆网站提供的开源软件MaxEnt(3.3.3k版本)[32]构建最大熵模型,主要包括以下3个步骤:①利用美国本土种植点潜在气候因子对MaxEnt模型进行多次训练,通过交叉验证评价各气候因子对模型的影响,筛选得到模型主导因子;②利用种植点主导因子重新训练模型,并对其进行精度评价;③基于覆盖中国云南省的主导因子插值数据,将MaxEnt模型进行外推。
1.2.4 气候适宜性区划模型的改进
作为引进树种,薄壳山核桃在云南省种植规模相对较小,种植点分布数据较少,缺乏经纬度信息,难以满足MaxEnt模型构建需求。而美国本土作为薄壳山核桃的原产地和主要种植区[3],种植点分布样本充足,并具有较强代表性。因此,本文基于美国本土的种植点样本训练模型,并将模型外推至中国云南省。Elith等[33]在针对外来物种的分布预测研究中指出,如果外推区域存在部分环境因子明显偏离训练样本取值范围,则应谨慎对待MaxEnt外推结果,同时引入多元环境相似度指标对模型的可靠性进行量化。云南省气候类型多样,模拟可靠性问题不容忽视。借鉴上述思路,本文提出调整型最大熵适宜性指数,利用外推区域与训练样本气候变量的偏离程度对MaxEnt模拟结果进行非线性抑制,其表达式如下:
Ia=Ip·R。
(2)
式(2)中,Ia为调整型最大熵适宜性指数,取值范围为0~1;Ip为存在概率因子,即MaxEnt模型输出值,取值为0~1;R为模型可靠性因子,可通过以下方法计算:
(3)
(4)
R=minRi。
(5)
假定模型包含n个气候因子,利用式(3)计算外推格点每一个气候因子与训练样本中该因子取值范围的偏离程度,其中,Vi为因子i在某一格点的取值,Vmax,Vmin分别为因子i在训练样本中的最大值和最小值。在此基础上,采用式(4)将Di映射为S型曲线。式中,Ri为气候因子的可靠性指标,随Di单调递增。当因子值Vi未超出训练样本取值范围,则Di=0,Ri≈1;当Vi略偏离训练样本取值范围时,Ri缓慢降低,随着偏离程度的增加,Ri下降幅度先逐渐增大,后逐渐收敛于0,如图1所示。根据式(5)在所有气候因子中取最小的Ri作为模型可靠性因子R。通过因子R的引入,在所有气候因子均未明显偏离训练样本的外推格点,Ia与Ip取值接近;而对于气候因子明显偏离训练样本的格点,Ia将被充分抑制。
图1 气候因子值(Vi)与可靠性因子(Ri)的映射关系示意图Fig.1 The mapping relation between climatic variable(Vi) and reliability factor(Ri)
气候要素间的自相关性会为MaxEnt模型带来冗余信息,影响模型的训练效果和运算效率。将薄壳山核桃种植点样本随机划分为两个子集,包括训练集样本204个和测试集样本70个。将训练集样本和覆盖美国本土的12个潜在气候因子插值结果输入MaxEnt软件构建初始模型,采用5折(5-fold)交叉验证[34]的方法进行检验,通过贡献率、置换重要性和刀切法(Jackknife)[35]分析结果,对影响模型的主导因子进行筛选。其中,贡献率体现因子对模型的影响程度,置换重要性体现因子的不可替代性(表1)。刀切法则采用仅保留单个因子和剔除单个因子两种策略对模型的增益变化进行比较,仅保留单个因子的情况下,增益较大的因子包含更丰富的有效信息;剔除单个因子情况下,增益较小因子包含更多其他因子没有承载的信息(图2)。结果表明:在气温相关的因子中,7月平均气温对模型的主导意义最明显,各项评价指标均优于其他11个因子,体现出夏季高温对薄壳山核桃种植区分布具有显著影响。年平均气温、活动积温和无霜期都是热量条件的综合体现,在各项评价中水平接近,均表现出相对较强的主导性。根据Pearson相关系数计算结果,年平均气温、活动积温和无霜期3个因子之间相关显著,故仅保留贡献率相对较高的年平均气温。30年极端最低气温和1月平均气温贡献率仅为0.4%,置换重要性处于中等水平,对模型影响较小。然而,上述因子均包含冬季低温条件和越冬风险相关信息,具有重要的生理学意义,不宜直接舍弃。因二者同样具有显著相关性,根据置换重要性和刀切法分析结果,仅保留30年极端最低气温。
在降水方面,3—5月降水量和年降水量分别在贡献率和置换重要性方面优势明显,对薄壳山核桃种植具有显著影响。6—7月降水量和8—9月降水量在各项评价中表现一般,在外推区域与年降水量具有显著相关性,对训练和外推过程均影响较小,故将其剔除。在日照方面,年日照时数贡献率达20.1%,仅次于7月平均气温,其置换重要性和刀切法分析结果也处于较高水平。4—5月日照时数在剔除单个因子的刀切法分析结果中水平较高,但置换重要性仅为2.6%。考虑到年日照时数综合反映光照资源,而4—5月日照时数侧重于授粉条件,生理学意义难以相互包含,二者应同时保留。
表1 潜在气候因子贡献率和置换重要性Table 1 Contribution and permutation importance of each potential climatic variables
图2 潜在气候因子的刀切法分析图Fig.2 Jackknife analysis of each potential climatic variables
综上,用于构建模型的主导因子最终包括7个:7月平均气温、年平均气温、30年极端最低气温、年降水量、3—5月降水量、年日照时数和4—5月日照时数。
基于上述7个主导气候因子,采用204个训练集样本重新训练最大熵模型,并利用70个测试集样本对模型进行验证。模型主要选项包括随机生成背景样本为10000个,最大迭代次数为500次,收敛阈值为10-4,启用clamping机制,输出类型采用对数变换结果。采用受试者工作特征曲线(ROC)下的面积(A)对建模结果进行评价(图3),评价标准如下:A<0.60(失败),0.60≤A<0.70(较差),0.70≤A<0.80(一般),0.80≤A<0.90(好),A≥0.90(非常好)[14]。结果表明:训练集A值为0.871,达到 “好”的标准,体现模型训练效果良好,可信度较高。测试集A值为0.820,略低于训练集,同样达到“好”的标准,表明模型具有良好的泛化能力,对训练样本以外的数据仍具有较高的模拟精度。
图3 最大熵模型ROC曲线Fig.3 ROC curve of MaxEnt model
利用种植点经纬度坐标提取气候要素值,并对样本中7个主导因子的最小值与最大值进行统计。结合云南省气候要素插值数据,根据式(3)~式(5)逐格点计算可靠性因子R(图4)。结果表明:云南省大部区域R>0.9,7个主导因子均未明显偏离训练样本取值范围,MaxEnt模拟结果具有较好的可靠性。R≤0.5的区域主要出现在滇西北、滇东北和滇东南局部区域,主要原因在于滇西北西部、滇东北北部和滇东南东部地区年日照时数、4—5月日照时数远未达到训练样本下限,而滇西北北部高海拔地区年平均气温、7月平均气温也明显低于训练样本取值范围。
图4 云南省可靠性因子(R)空间分布Fig.4 The distribution of reliability factor(R) in Yunnan Province
将覆盖云南省全境的7个主导气候因子输入训练后的MaxEnt模型,生成存在概率因子Ip,即改进前的适宜性指数。在此基础上,利用式(2)计算调整型最大熵适宜性指数Ia,即改进后的适宜性指数(图5)。由图5可以看到,Ip与Ia总体趋势相似,受立体气候特征影响空间分布较为破碎。云南省北部的高原气候带和温带区域气温类要素较种植点样本的平均水平低,Ip与Ia取值普遍低于0.1。随着海拔高度的逐渐下降,温度条件的改善对上述适宜性指数的增大起促进作用,日照与降水对适宜性指数影响也越来越明显。在图5a中,滇东北北部和滇东南东部的河谷地带,因热量、水分条件良好,MaxEnt模型输出较高的Ip值。根据2.3节分析,上述区域与薄壳山核桃原生境相比,光照条件明显偏差,不利于薄壳山核桃生长发育,直接应用MaxEnt模型的跨区域外推结果Ip进行区划,可能导致部分区域适宜性的误判。相比之下,通过引入可靠性因子R,改进后的气候适宜性指数Ia(图5b)在上述区域受到充分抑制,从作物生理学角度,改进结果更合理。
根据相关文献[9]报道,对云南省薄壳山核桃的主要种植县(市、区)进行整理,如图6所示。相比图5可见,Ia值较高的区域主要分布于滇中及以西、以南地区,与主要种植县(市、区)在空间分布具有较好的对应关系。在滇东北北部和滇东南东部地区,Ip较Ia明显偏高,且出现较大范围Ip>0.4的区域,但该地区鲜有种植案例。由此可见,改进后的适宜性指数Ia更符合实际种植情况,进一步体现改进方法的有效性。
图5 改进前(a)与改进后(b)气候适宜性指数的对比Fig.5 Comparison of climate suitability index before improvement(a) and after improvement(b)
图6 云南省薄壳山核桃主要种植县(市、区)分布Fig.6 Distribution of major pecan planting counties in Yunnan Province
针对云南省全境调整型最大熵适宜性指数计算结果,采用自然断点法将薄壳山核桃种植气候适宜性划分为4个等级:Ia≥0.35为最适宜,0.20≤Ia<0.35为适宜,0.08≤Ia<0.20为次适宜,Ia<0.08 为不适宜(图7)。由图7可以看到,薄壳山核桃种植的气候适宜性在空间分布上与热量资源关系最密切,最适宜区和适宜区多位于云南省南部和中部海拔较低的区域,次适宜区和不适宜区主要出现在海拔较高或日照时数偏低的区域。
最适宜区占云南省总面积的5.3%,沿金沙江、元江、阿墨江、大盈江、瑞丽江和南定河等流域及其支流零散分布。该区域年平均气温为15.8~21.9℃,7月平均气温为20.8~26.4℃,30年极端最低气温为-6.4~2.0℃,热量条件优越,霜冻害风险小,休眠期低温条件较好。3~5月降水量为65.7~350.6 mm,年降水量为636.2~1633.9 mm,均具有较大的变化范围,年降水量高于Spark[8]提出的原生境降水量下限。年日照时数为1633.5~2570.5 h,4—5月日照时数为333.3~490.3 h,光照条件良好。分布于云南省各地的最适宜区均具有热量资源丰富的特点,但其他条件存在较大差异,如北部金沙江及其支流龙川江沿线的最适宜区降水偏少,但日照充足;与之相反,南部的最适宜区降水充沛,但日照条件略差。从地貌的角度,最适宜区多分布于河谷和盆地,但元江河谷中部气候干热,冬季低温不足和降水偏少同时制约河谷地带的适宜性,河谷两侧海拔1800 m以下的山地更适宜薄壳山核桃种植。
图7 云南省薄壳山核桃种植气候适宜性区划Fig.7 Climate suitability regionalization of pecan plantation in Yunnan Province
适宜区占云南省总面积的23.1%,主要分布于滇中及以南低海拔地区,部分区域沿河谷向北延伸。该区域年平均气温为13.7~23.9℃,7月平均气温为19.1~26.5℃,30年极端最低气温为-8.0~1.5℃,年降水量为769.0~1804.2 mm,3—5月降水量为78.3~342.2 mm,年日照时数为1503.4~2449.7 h, 4—5月日照时数为298.6~476.8 h。适宜区气候类型丰富多样,总体而言能为薄壳山核桃种植提供较为良好的条件,但与最适宜区相比,气候要素综合特征与原生境差异较大。部分区域受到单方面气候条件不利影响,如分布于滇中的适宜区热量条件较差,滇南一带的适宜区休眠期低温条件较差,滇东南东部则主要受到日照偏少的影响。
次适宜区占云南省总面积的30.7%,分布范围较广,根据气候条件可分为以下两类:一类分布于北亚热带和中亚热带,在纬度较高的滇西北、滇东北和滇中地区位于河谷和平原地带,在滇中以南地区则以山地和丘陵为主,上述区域年平均气温为13.0~18.0℃,虽然薄壳山核桃植株仍能存活,但热量条件不足易导致结实年限长、坚果质量不佳等问题[9],引种果用树种应慎重考虑。另一类次适宜区分布于滇东南、滇西南东部,以及滇西北南部的怒江河谷,具备良好的热量条件,但日照不足明显制约该区域的气候适宜性;部分区域处于北热带,冬季过于温暖也对发芽具有不利影响。
不适宜区占全省总面积的40.9%,集中分布于滇东北大部、滇中东北部和滇西北西部,也出现在其他区域的山脊地带。其中,以青藏高原南延和滇东北为代表的高海拔地区,热量条件差,夏季高温不足,越冬风险较大,不利于薄壳山核桃生长发育。云南省东部和西北部边缘海拔较低的区域,光照条件与原生境差距较大,也不适合发展种植产业。
目前,薄壳山核桃种植产业在大理和玉溪已得到较好发展,形成多个产量较高的种植点,如漾濞瓦厂乡、永平杉阳和新平华森果业公司等[9],高产种植点在分布上与最适宜区和适宜区基本吻合。普洱、德宏、临沧、楚雄、丽江和红河等州市具有较大范围的最适宜区和适宜区,部分地区已推广种植[3,9],可适当扩大种植规模。此外,云南省地形复杂,降水有效性不高[36],引种过程中还应充分考虑灌溉条件。
本文基于MaxEnt模型,利用气候与种植点分布数据计算薄壳山核桃适宜性指数,并对中国云南省开展适宜性分区,结果表明:
1) 通过交叉验证,影响薄壳山核桃种植的主导气候因子包括7月平均气温、年平均气温、30年极端最低气温、年降水量、3—5月降水量、年日照时数和4—5月日照时数。经验证,基于上述因子构建的MaxEnt模型具有较高的精度和泛化能力。
2) 通过引入可靠性因子,根据模拟区域和训练样本气候因子值域的偏离程度,对MaxEnt模型输出结果进行非线性调整,实现区划模型的改进。对于部分气候因子明显偏离训练样本取值范围的外推区域,该方法能够有效抑制适宜性的高估,并提高跨区域外推结果的可靠性。
3) 云南省薄壳山核桃种植气候最适宜区和适宜区分别占云南省总面积的5.3%和23.1%,主要分布于热量资源丰富、日照相对充足并具备较好冬季低温条件的亚热带地区和热带地区边缘。受云南省复杂地形与气候条件影响,区划结果呈现破碎化特征,在薄壳山核桃引种过程中,应该充分考虑立体气候的影响,尤其要重视干热河谷地带气候资源的合理利用。
本文直接采用薄壳山核桃物种种植点分布信息开展气候适宜性研究,但该物种对气候条件的适应能力较强,随着品种选育和种植技术的不断发展,薄壳山核桃对气候条件的适应范围逐渐扩大,且不同品种对气候条件的需求存在一定差异[8,37]。在未来研究中,有必要对品种进行区分,将气候适宜性区划细化到部分主栽品种。此外,还应关注气候变化对农业气候资源的影响[38],对未来气候变化背景下薄壳山核桃气候适宜性布局开展预测研究。