基于MaxEnt和ArcGIS预测广东杉木良种在湖北省同一适宜引种生态区*

胡 超 于 静

（1.湖北省林业局林木种苗管理总站，湖北 武汉 430079； 2.岭南生态文旅股份有限公司，湖北 武汉 430062）

林木良种是有适宜生态区域要求的，如果自然条件不适宜，再好的良种也达不到丰产、稳产。因林木良种不适应引种区自然条件而造成巨大损失的教训是深刻的：20 世纪70 年代，各地在油茶Camellia oleifera生产发展过程中调购种子比较随意，较多地方因为超地理区域引种造林，引种前没有进行科学预判，盲目性的引种，导致幼林生长不良、成林产量很低，在人力、物力等方面都造成了较大的损失[1]。

传统的林木良种引种适宜生态区凭主观经验判断较多，如，宜林范围内每个按水平分布的气候带和垂直气候带都分布着特有类型的森林植被。经纬度由北向南，由西向东调运范围大于相反方向的范围，海拔不超过300～500 m，但是，1958年，湖北引种广东、福建马尾松Pinus massoniana种子成功，用事实改变了过去专家认为“马尾松南种北移的幅度不能超过2°～3°”的定论[2]。1979年李传志[2]论证马尾松一次北移6°～7°育苗是可以成功的。所以，温度、降水、土壤等主要环境因子相似，即为林木良种同一适宜引种生态区。

杉木Cunninghamia lanceolata是湖北省主要造林树种之一。传统的杉木良种引种适宜生态区也是凭主观经验判断较多。杉木良种数量较多，且生长周期长，像农作物良种一样，对所有杉木良种都进行引种试验的可行性不大，所以，采用较为可靠的预测方法对引种的可行性提供数据参考尤为必要。目前，可用方法包括了MaxEnt 和ArcGIS 分析法等[3]。鉴于此，本研究拟广东杉木良种为引种对象，应用MaxEnt 生态学模型[4-5]，以100×100=10 000 m2，即1 hm2为单元，用34个环境因子预测广东杉木良种在湖北省同一适宜引种生态区，以期为引种应用提供参考。

1 数据与方法

1.1 基本信息来源

广东省审定杉木良种信息来源于湖北省林业局林木种苗管理总站（表1）。

表1 广东省杉木产区审定杉木良种Tab. 1 Cunninghamia lanceolata superior varieties approved by Cunninghamia lanceolata production area in Guangdong province

34 个环境因子数据获取于中国气象科学数据共享服务网、中国科学院资源环境科学数据中心、国家青藏高原科学数据中心、中国西部环境与生态科学数据中心（表2）。

表2 广东省杉木良种在湖北省同一适宜引种生态区环境因子Tab. 2 Environmental factors of the Hubei province identical suitable introduction ecological distribution of superior varieties of Cunninghamia lanceolata in Guangdong province

中国行政区划数据、中国海拔（DEM）数据

获取于中国科学院资源环境科学数据中心和湖北省林业调查规划院。

1.2 信息数据处理

1.2.1 分布数据处理 为避免样点数据在某个地理空间上过度聚集，在广东省杉木适生范围内，用ArcGIS 的Create fishnet 工具生成空间为40 行×40 列的格网数据，以1 个格网作为1 个采样单元对杉木良种的分布数据进行采样（图1）[6]。根据选育单位确定的杉木良种适宜的自然地理环境条件范围，如，适宜海拔范围为1 500 m 以下，在Excel 表中，剔除高程小于0 m、高程大于1 500 m、土壤厚度小于30 cm 和异常值的采样点，全部采样分布记录共320 条。按照MaxEnt 软件的“Samples”的要求整理数据，将分布点以“物种+经度+纬度（西经、南纬的值为负，经纬度为十进制小数格式。）”另存为CSV 格式文件。

图1 广东省杉木产区采样点分布Fig. 1 Distribution of sampling points in Cunninghamia lanceolata production area of Guangdong province

1.2.2 环境因子处理 地形因子（经度、纬度、高度）与环境因子有较好的回归关系，利用中国2 160 个基本、基准地面气象观测站的观测数据（1981—2010 年），推算模拟无测站区域的环境资源分布情况。建立Bio1～Bio10、Bio13～Bio21、Bio24～Bio27 等23 个环境因子的空间分布模型，其表达式为：

式中，Y为环境因子要素；λ为经度；φ为纬度；h为海拔(m)；函数f(λ,φ,h)为气候学方程；ε为残差项，可视为小地形因子(坡度、坡向等)及下垫面对环境的影响。将f(λ,φ,h)展成三维二次趋势面方程[7]。

式中b0～b9为待定系数，利用SAS9.4 建立逐步回归优化回归模型，模拟23 个环境因子的宏观趋势项，分别建立23 个环境因子的小网格推算模型（表3）。

表3 环境因子的小网格推算模型Tab. 3 Small grids reckoning models of regionalization indexes of environmental factors

在中国海拔（DEM）数据支持下，在ArcGIS里，用23 个环境因子的小网格推算模型，将环境因 子Bio1～Bio10、Bio13～Bio21、Bio24～Bio27 分别插值为100 m × 100 m 网格的基础数据[8-9]。用IDW 法分别插值其残差项为100 m × 100 m 网格的修正数据。用Spatial Analyst 工具→数学→逻辑→加，将每个环境因子的基础数据和修正数据叠加相加为环境因子栅格数据。23 个环境因子栅格数据用投影栅格工具统一为地理坐标系D_WGS_1984。以湖北省和广东省矢量边界为掩膜，裁剪出这23 个环境因子栅格数据图层。最后，用栅格转ASCII 工具将这23 个环境因子栅格数据转换保存为MaxEnt 所需要的ASCII 格式文件。

在ArcGIS 里，将下载的Bio11、Bio12、Bio22、Bio23、Bio28～Bio34 等11 个环境因子数据通过重采样工具使其像元大小与Bio1～Bio10、Bio13～Bio21、Bio24～Bio27 等23 个环境因子一致[10]。11 个环境因子数据统一为地理坐标系D_WGS_1984。以湖北省和广东省矢量边界为掩膜，裁剪出这11 个环境因子栅格数据图层。最后，用栅格转ASCII 工具将这11个环境因子栅格数据转换保存为MaxEnt 所需要的ASCII 格式文件。

1.3 模型构建

1.3.1 MaxEnt 软件建模 （1）物种数据：将之前导出的杉木良种分布数据(csv 格式)的文件，通过Browse 加载到MaxEnt 软件“Samples”模块。

（2）环境数据：把ASCII 格式文件的34 个环境数据加载到MaxEnt 软件“Environmental layers”模块。

（3）参数设置：使用auto features 选项，根据自动特征规则进行计算，所有的要素类型都将用到。结果以comulative 类型和ASCII 格式输出，并定义其输出位置。设置界面的选择 settings里‘Random test percentage’ 设 置 为25，随 机 选取75%的样本点数据作为训练数据[11]，‘replicates’选择3 次重复作为平行试验，最大迭代次数设为500 次，收敛阀值设为0.00 001，取值范围0～100[12]。 选 择‘Do jackknife to measure variable importance’衡量所有变量的重要性，MaxEnt 软件分别对每一个环境影响因子进行刀切图绘出。

1.3.2 ROC 曲线绘制 绘制响应曲线(response curves)评价模型精度。ROC 曲线以真阳性率为纵坐标(敏感性，实际存在且被预测为存在的比率)，以假阳性率(1-特异性，实际不存在但被预测为存在的比率)为横坐标，AUC 值指 ROC 曲线与横坐标围成的面积值，值域为0～1。AUC 值越大表示与随机分布相距越远，环境因子变量与预测的物种地理分布之间的相关性越大，即模型预测效果越好，反之说明模型预测效果越差。AUC 值在 0.5～0.6，0.6～0.7，0.7～0.8，0.8～0.9，0.9～1 分别表示模拟效果失败、较差、一般、好、非常好[4,8]。34 个环境因子预测模型的训练样本和测试样本的AUC 值 达 到0.802 和0.739( 图2)，AUC 均 值 在0.7～0.8 之间，说明模型预测效果一般。

图2 初始模型的ROC 曲线分析及AUC 值Fig. 2 ROC curve and AUC value for the initial model

1.4 模型优化

在使用MaxEnt 模型进行较大空间范围的杉木良种同一适宜引种生态区预测时，如果环境因子变量过多、变量空间共线性过强，将导致模型的复杂性增加，随机误差增大。所以，过多低贡献率的环境因子变量会导致模型运行结果的准确性降低。因此，需要对环境因子进行筛选或降维[13]。

1.4.1 筛选贡献率高的环境因子变量 在34 个环境因子中，对于杉木良种同一适宜引种生态区分布贡献较大的环境因子变量有：Bio3、Bio6、Bio8、Bio10、Bio11、Bio17～Bio19、Bio31，累计贡献率为93.8%。Bio1、Bio 2、Bio 4、Bio 5、Bio 7、Bio 9、Bio12～Bio16、Bio20～Bio30、Bio32～Bio34 等25 个环境变量的贡献率都小于1%（表4），对杉木的种植分布影响有限，对这25 个环境因子变量进行剔除[14]。

表4 各环境因子变量的贡献率Tab. 4 Contribution rate of each environmental factor variable

1.4.2 筛选正规化训练增益高的环境因子变量刀切法(jackknife test)测定各环境因子变量权重。刀切法就是每次都忽略一个环境因子变量，然后基于剩下的环境因子变量来对杉木良种同一适宜引种生态区进行预测，然后MaxEnt 软件自带程序画出柱形图作为依据评估环境因子变量的重要性。红色条带代表所有变量的贡献；深蓝色的条带越长，说明该变量越重要；浅蓝色的条带长度代表除该变量以外，其他所有变量组合的贡献。Bio3、Bio6、Bio8、Bio10、Bio11、Bio17～Bio19、Bio31 等6 个贡献较大的环境因子变量中，Bio31 对应的深蓝色条带很短（图3），说明它本身的增益值几乎接近于0，表明它对预测杉木良种同一适宜引种生态区并不是重要环境因子变量，所以，剔除环境因子变量Bio31。

图3 刀切法的环境因子变量重要性分析Fig. 3 Analysis of the importance of environmental factors variables in the Jackknife method

1.4.3 筛选多重共线的环境因子变量 用GIS 软件的值提取至点工具提取有效分布点的环境因子变量数值，用SPSS 软件对贡献较大的Bio3、Bio6、Bio8、Bio10、Bio11、Bio17～Bio19 等9 个主导环境因子进行Spearman 相关分析（表5），检验环境因子变量之间的多重共线性，Bio8 分别与Bio3、Bio6、Bio10 和Bio11 的相关系数|r|≥0.8，分别对比初始模型中二者的贡献率，Bio3、Bio6、Bio10 和Bio11 贡献率较小，所以，剔除贡献率较小的变量Bio3、Bio6、Bio10 和Bio11，提高模型模拟的精度[14]。

表5 关键环境因子变量的相关系数Tab. 5 Correlation coefficient of key environmental factor variables

2 结果与分析

2.1 模型的有效性

用剩余的Bio8、Bio17～Bio19 等4 个主导环境因子变量重新建模，重建模型的训练样本和测试样本的AUC 值达到0.7563 和0.750(图4)，AUC均值在0.7～0.8 之间，表明重建模型适用性及模拟精度均达到合格水平，与主导环境因子变量之间的相关性较大，预测同一适宜引种生态区的结果合格，可以据此进行引种推广。

图4 重建模型的ROC 曲线分析及AUC 值Fig. 4 ROC curve and AUC value for the reconstruction model

2.2 杉木良种同一适宜引种生态区预测及适生等级划分

MaxEnt 进行3 次重复试验，选取重复试验中，AUC 值最高的图层导入ArcGIS 软件进行适宜等级划分和可视化表达（图5）。MaxEnt 模型输出的数据为ASC Ⅱ格式，用ArcGIS 的ASCII to Raster 功能，输出数据类型选FLOAT，使该结果可在 ArcGIS 中显示[14]。利用“Reclassify”功能，划分分布值等级及相应分布范围，并使用不同颜色表示，划分标准为： 存在概率＜0.05 为不适生区；0.05 ≤存在概率＜0.33 为低适生区； 0.33 ≤存在概率＜0.66 为中适生区； 存在概率≥0.66 为高适生区[4,8]。整体来看，广东省杉木产区的杉木良种在湖北省的低适生区面积为1 294 849 hm2，主要分布在：鄂东南的黄梅县、武穴市、蕲春县、阳新县、大冶市、通山县、崇阳县、通城县、咸安区和赤壁市；鄂西南的长阳县、秭归县、兴山县、巴东县、建始县、鹤峰县、恩施市、宣恩县和来凤县。低适宜区域在引种杉木良种时，需要选择适宜的小生境。

图5 广东省杉木良种在湖北省同一适宜引种生态区分布Fig. 5 Distribution of the Hubei province identical suitable introduction ecological distribution of superior varieties of Cunninghamia lanceolata in Guangdong province

2.3 杉木良种地理分布与生物气候变量的关系

用刀切法(Jackknife Test)检测4 个主导环境因子变量对于分布增益的贡献，结果(表6)表明：年日最低气温≤0.0℃日数(Bio8)对杉木分布的增益最大，分布值随年日最低气温≤0.0℃日数的升高而减小(图6)。月最长连续降水日数(Bio18)也对杉木分布的影响较大，当月最长连续降水日数为14～21 日，分布值随月最长连续降水日数的升高而增大 (图7)。月最长连续降水量(Bio19)也对杉木分布的影响较大，当月最长连续降水量为118.971～280 mm，分布值随月最长连续降水量的升高而增大；当月最长连续降水量为280～883.572 mm，分布值随月最长连续降水量的升高而减小(图8)。

图6 年日最低气温≤0.0℃日数(Bio8)反馈曲线Fig. 6 Days with annual daily minimum temperature below 0℃ (Bio8) feedback curve

图7 月最长连续降水日数(Bio18)反馈曲线Fig. 7 Monthly longest continued precipitation days(Bio18) feedback curve

图8 月最长连续降水量(Bio19)反馈曲线Fig. 8 Monthly maximum continued precipitation (Bio19)feedback curve

表6 主导环境因子变量的贡献率Tab. 6 Contribution rate of dominant environmental factor variable

3 讨论与结论

如果没有做到“适地区适良种”，则可能导致林木良种育苗和林木良种造林的失败。主要原因是不同地理、气候和土壤等环境因子的质或量不同，因而对该林木良种所要求的生态条件的满足程度不同。所以，只有在适生地域内，在适宜的立地条件下，选择适合的林木良种，才能发挥林木良种造林的优良特性，实现速生丰产。否则，林木良种表现不好，甚至不如一般的当地品种。

基于MaxEnt 生态位模型的同一适宜生态区研究中，环境因子数据常来自于世界气候—全球气候数据库网站，仅有19 个环境因子，空间分辨率仅为5arc-min[5-6,10,14]。本研究选取34 个重要环境因子，用中国2 160 个基准地面气象观测站的观测数据，推算模拟无测站区域的环境资源分布情况，相较于全球气候数据库环境因子数量更多，数据点更密集。

传统的杉木良种引种同一适宜生态区都是以乡镇、县、市、省等行政单位为单元。然而，影响杉木成活生长的光、热、水、气等环境因子，受太阳辐射、大气环流的影响，在地面上呈地带性的分布。由于山体起伏，垂直森林地带在实际上并不都是连续的，而是由断断续续地呈孤岛状分布的地块组成。所以，本研究是以100×100=10 000 m2，即1 hm2为单元。

传统的林木引种是以单个树种划出同一适宜生态区[5-6,10,14-16]。然而，随着自然条件演变和科学技术的发展，转抗性基因育种、种间和远缘杂交育种等遗传改良工作在广泛开展，每年都有新的林木良种通过审定。在相同的立地条件下，同一树种，不同良种之间的生长好坏是有显著差异的。为了获得精准的引种效果，本研究以单个良种划出同一适宜生态区精准预测广东省杉木良种在湖北省同一适宜生态区。

通过运用MaxEnt 生态位模型对广东省杉木良种在湖北省同一适宜引种生态区进行分析，证明了MaxEnt 模型在林木良种引种应用方面的可行性以及可信度，同时结合刀切法探讨对杉木良种生长影响最显著的环境因子，这对广东省杉木良种适生性分析提供了更进一步的技术支撑。如同诸葛亮需要精准的预测天气，才能草船借到箭一样。还如同哥伦布需要精准的指南针，才能航海发现新大陆一样。需要基于MaxEnt 和ArcGIS 精准预测广东省杉木良种在湖北省同一适宜引种生态区，才能避免引种广东省杉木良种的盲目性，从而获得良好的引种效果，值得推广应用。