基于组合物种分布模型（Ensemble Model）的厚朴适宜生境分布模拟

朱 妮

（1.安康学院旅游与资源环境学院，陕西安康 725000；2.秦巴国土资源利用与环境保护协同创新中心，陕西安康 725000；3.安康市汉江水资源保护与利用工程技术研究中心，陕西安康 725000）

厚朴（Magnolia officinalis Rehd.et Wils）为木兰科落叶乔木，在我国主要分布于南方各省[1-2]。该物种是中医药体系常用药，主要药用部位为其树皮，药用价值极高。厚朴在我国的药用历史已有两千余年，现代药理学研究证实，该药材具有抗菌、抗病毒、抗过敏等功效，同时对于胃溃疡，应激性胃功能障碍等症状有一定治疗和缓解作用[3-4]。同时实验表明厚朴的主要药用成分厚朴酚、和厚朴酚等对癌细胞具有抑制作用[5]。厚朴经济价值高，但是由于药用厚朴的生长周期较长，一般需要15年其树皮才能达到药用标准，由于不可持续的采收方式以及森林的砍伐，目前该物种野生资源急剧减少，生境破碎度较高[6]。另一方面，厚朴在地区生态恢复以及生态建设方面也有优势，该树种根系较为发达，生长速度快，且具有一定观赏价值，可以在一定区域进行推广。鉴于此，预测厚朴适宜生境分布，可以为厚朴的人工种植，野生资源保护等提供理论依据。物种分布模型（species distribution model，SDM）即生态位模型，是利用统计学算法在多维环境空间中模拟物种生态位范围，并将物种生态位范围抽象成统计学规则，模拟物种在不同时空环境下的适宜生境潜在分布的数学算法集合[7-8]。近年来，随着统计技术与生态位理论的发展，物种分布模型在生态学中的应用迅速增加，目前此类模型已经广泛地用于生物地理学、保护生物学、全球变化生物学与生境或物种管理等方面[9-13]。目前已有大量的统计算法被用于物种分布模型构建，在建模过程中针对不同的建模目的，物种生态位特征以及建模数据基础，研究者需要选择不同模型算法[14]。目前这一选择主要是靠一种或几种（如：kappa 系数、TSS 值、ROC 等）统计学评价指标来确定，但是最优的模型算法不一定能在所有建模条件下统计学精度能达到最优[15]，克服模型应用的不确定性是亟待解决的问题[16]。研究表明组合物种分布模型（ensemble model）能够在一定程度上解决这一问题。该模型可以将模型的主要趋势（即平均值、中位数或其他百分位数）和所有模型的总体变化（以及不确定性）映射起来；也可以集成不同模型的其他方面，如变量的重要性或模型响应曲线[17-18]。

BIOMOD 是目前最为流行的物种分布单模型集合预报平台，2003年W.Thuiller[19]提出了一个新型计算机平台框架BIOMOD（BIOMOD：BIOdiversity MODelling），其目的是利用不同类型的统计建模方法最大化当前物种分布的预测精度和增强物种未来潜在分布预测的可靠性。但是早期BIOMOD只集成了4 种单模型算法即广义线性模型（GLM）、广义相加模型（GAM）、分类与回归树模型（CART）以及人工神经网络模型（ANN）。2009年，W.Thuiller[20]进一步完善了BIOMOD 的计算框架，同时也发布了R 环境下免费且开源的BIOMOD2 程序包，目前BIOMOD 中集成在原有的4 种算法上又增加了5种算法即：表面分布区分室模型（SRE）、多元适应回归样条函数（MARS）、柔性判别分析（FDA）、推进式回归树（BRT）、随机森林（RF），同时该平台还集成了不同模型初始条件的设置，不同的模型参数化条件，例如：回归算法中设置不同阶数的多项式和光滑样条，在分类算法中设置分类树中的节点，以及在神经网络算法中设置隐藏层的数目以及公式等[20]。

本文基于大量文献确定物种分布记录，同时收集并筛选相关环境因子，采用组合物种分布模型，模拟我国南方厚朴的潜在适宜生境分布。组合分布模型要求使用较多的单模型算法构造以寻求物种分布模拟的最佳解决方案，减少模型不确定性，因此在本文中我们使用BIOMOD2 建模平台，构建其支持的全部9 个单模型算法（表面分布区分室模型、分类树分析、柔性判别分析、广义线性模型、广义加性模型、多元适应回归样条函数、人工神经网络、推进式回归树以及随机森林），依照模型TSS 值进行模型结果筛选与定权，构造组合物种分布模型（ensemble model）探究厚朴的分布规律，并比较不同物种分布模型的精度以及空间表现。

1 材料和方法

1.1 厚朴分布位点与研究区的确定

精确的物种分布数据是物种分布模型构建的基础，本研究中为了保证物种分布数据的数据真实可靠，所有的厚朴样点数据都来源于公开发表的文献[1-2，21-23]，同时为了保证厚朴地理位置的精确度及模型构建的科学性，本研究对获得的厚朴采样点数据按照以下条件进行了筛选：①只选取2000年之后的精确位置信息的数据；②删除经纬度重复的数据，并且样点间的距离大于10 km。经筛选处理后参与建模的厚朴采样点数据共143 个（图1）。

图1 研究区和种群分布点的地理位置Figure 1 Study area and geographic locations of populations

1.2 环境图层数据搜集以及数据预处理

完备且准确的环境因子是物种分布模型构建的基本要求，建模时关键因子的缺失会造成无法解释的方差，同时也可能导致建模过程中对物种生态位的量化产生偏差。在本文中我们采用物种分布模型构建最为常用的Bio（bioclimatic variables）数据系列，该数据系列一共包括19 个具有生态学意义的气候变量[24-25]。数据来源于自世界气候数据库（WOR LDCLIM version 2.0，http：//www.worldclim.org）。该数据由世界各地气象站1970—2000年记录的气温降水等信息结合地形信息插值生成。同时本文收集了研究区的高程数据，数据源于中国科学院地理科学与资源研究所资源环境数据云平台（http：//www.resdc.cn/data.aspx?DATAID=200），利用 ARCGIS 10.3 空间分析功能我们得到了坡度，坡向2 个环境因子。另外，为了进一步体现厚朴潜在分布的实际区域，本试验下载了中国1∶100 万植被类型数据图，数据同样源于资源环境数据云平台。本研究采用的环境因子空间分辨率均为30″（1 km2）。

前人研究表明，如BIO 数据系列存在较强的数据共线性，环境因子的相关性会造成建模时的统计信息的冗余，同时也会对建模的结果产生不利的影响[26]。因此必须对环境因子进行筛选，本文采用相关分析法[26-28]以及厚朴生境特征分析筛选参与建模的环境因子。首先对于收集到的22 个环境因子（19 个气候变量和3 个地形变量）求取相关系数，对于相关系数＞±0.7 的多组环境变量进行筛选，每组环境变量只选取一个与厚朴生境特征关系最为密切的变量，参与模型构建，同时检查随后筛选得到的环境变量的相关系数，确保其小于0.7[27-28]。最后本文选取：年平均气温（bio1）、平均日较差（bio2）、温度年较差（bio7）、年降水量（bio12）、降水量的季节性变化（bio15）、最冷季降水量（bio19）、坡度（slope）、坡向（aspect）共8 个环境变量参与本文模型构建（表1）。

表1 厚朴潜在地理分布评价所需的地理环境变量Table 1 Geographical elements index used for predicting potential geographic distribution of M.officinalis

1.3 BIOMOD单模型构建

在BIOMOD2 运行时候需要物种的存在点位数据，也需要物种的不分布点位数据（伪存在点数据）。在本文中BIOMOD2 的设置如下：首先根据厚朴存在点位数据与背景点位（研究区范围内没有存在点位的栅格）数据随机生成3 组伪存在点数据，每组伪存在点数目均为500。然后将参与建模的样点数据（包括存在点与伪存在点）随机分为2 个部分，70%的数据用为训练集，剩余数据为测试集，同时我们设置在模型运行以及模型评价时，采样点的总权重与伪采样点权重相等。然后将数据输入并运行9 个单模型算法，每个模型算法重复10 次，本研究将一共产生270 个单模型运算结果（3 组样点数据*9 个模型算法*10 次重复）。最后，采用TSS 和AUC 值对模型进行评价，TSS 的值是指模型正确预测存在点的百分比与正确预测伪存在点百分比之和减去1，TSS 值域范围均为[-1-1]，取值越接近1，模型精度较高，一般情况下，取值＞0.8，表示模型结果较为理想，取值0.5 以下表示模型结果较差[27]。受试者工作特征曲线（receiver operating characteristic curve，ROC 曲线）是目前物种分布模拟结果验证使用最为广泛的指标。AUC（area 大于0.8 时，表示模型结果较为精确[28]。

1.4 组合物种分布模型构建

为了减少模型本身以及数据生成（主要是伪采样点）过程中的不确定性，本研究采用组合物种分布模型来预测厚朴适宜生境的空间分布。综合考虑本文中所构建的所有270 个单模型运算结果，本文中采用TSS 值进行筛选构成最终数和模型的单模型结果，只保留TSS＞0.8 的单模型运算结果采用加权平均法构建组合物种分布模型，单模型运算结果权重由其TSS 值决定，公式如下（1）：

其中，wj为第j 个单模型结果的权重，rj表示第j 个单模型的 TSS 值，h 表示 TSS＞0.8 的单模型结果数目。

最后将单个模型归一化后的结果与对应权重相乘，然后求和，以此构建组合模型并计算研究区厚朴潜在适宜生境指数yi，计算公式为（2）：

式中：yi是第i 个栅格的厚朴适宜生境指数；wj为第j 个模型的权重；xij为第j 个模型中的第i 个栅格的数值。yi的取值范围[0，1]，yi越接近 1，表明该栅格单元的环境条件条件越适宜厚朴生长。

为了更直观地显示物种分布预测结果，需要将本文连续型的适宜性指数分布结果通过定义阈值转换为布尔型结果，即超过阈值为1，低于阈值为0，本文中采用最大TSS 值法进行阈值设定，TSS 值的计算依赖于阈值的设定，不同的阈值对应不同的TSS 值，本文中采用使得组合模型TSS 值为最大的分类阈值对建模结果进行分类。

同时本文提取中国1∶100 万植被类型数据图中的阔叶林、针阔叶混交林、落叶果树园以及亚热带经济林等植被类型作为厚朴分布限制范围与组合模型的二值化运算结果求交集，以确定厚朴实际潜在分布区的位置以及面积。

2 结果与分析

2.1 不同模型预测结果精度

本文采用TSS 和ROC 这2 个模型评价指标对参与建模的9 个单模型算法进行精度检验（图2），结果表明，对于不同的评价指标，最优模型不同，就ROC 而言，最优模型为RF，其平均 AUC 值为0.95，而9 个模型平均TSS 最大的模型为MARS 其平均TSS 值为0.87，这在一定程度上说明了单模型统计精度的不确定性，不同的统计精度指标给出了不同的结果。综合2 个评价指标来看，本研究中性能最优的模型为RF 与MARS，预测性能较优的模型为：GAM、GBM、FDA 以及 GLM。相比以上模型，ANN 的精度较差，其平均TSS 值为0.82，模型精度最低的是SRE 与CTA，其TSS 值均为超过0.8，这表明这两个模型并未参与组合模型的构建。同时图2 中每个模型在不同评价指标数据轴上延伸的范围代表其多次建模精度评价的稳定性，我们可以看出9 种模型算法中精度稳定性最好的是ANN，其次为RF、MARS 与CTA，模型精度评价稳定性最差的是SRE与GLM，这也在一定程度上说明了模型结果的不确定性。我们计算了组合物种分布模型的TSS 与AUC值，结果表明，组合模型能够在一定程度上提高模型精度，其 TSS 值为 0.905、AUC 值为 0.975。

2.2 不同模型算法的空间表现

本研究中为了展现不同算法的空间表现，将所有单模型算法的结果分别相加，然后归一化处理，得到不同单模型对厚朴适宜性指数的空间分布（图3）。由图中我们可以看出，模型结果具有很高的一致性，所有的模型结果都显示厚朴高适宜指数分布区位于四川东部、甘肃南部、陕西南部、湖北西部同时在贵州、湖南、江西、福建、浙江等省份也有大量的分布区存在。但是在细节上，模型结果也有较大差异，由高适宜指数区域分布面积大小来看，SRE的模型结果最为保守，由于其模型结果为二值化结果，所以其分布区面积较为集中，而GLM 的结果高适宜指数分布区面积最大，几乎包括了研究区的大部分区域，其分布北界达到陕西中部、甘肃中南部，分布西界延伸到了四川西部以及云南省，分布南界到达广东南部，甚至在台湾南部也出现了一定面积的高适宜生境指数分布区域。

从算法分类上看回归模型（GLM、GAM、MARS）的中GLM 与MARS 结果较其他模型的结果都较大，回归模型的结果不同适宜性指数分布区的分布细节较为丰富。分类模型（CTA、FDA、RF）中 FDA 的空间模拟结果最差，高适宜生境区域被明显的分为两个隔离的区域，分类模型的结果中高适宜指数分布区集中分布，游离于主要集中分布区的小板块状的高适宜分布区较少。而复杂的机器学习算法如：ANN、GBM、RF 等空间模拟结果比较收敛，提供了生境分布的细节，具有较好的空间表现，同时依据模型性能评价（图2）大多数此类模型具有较高的模型精度。因此我们可以初步得出结论：在相同的数据基础上，复杂的机器学习模型有着较好的统计学模拟精度，同时其空间表现也优于简单的传统统计模型。综上所述我们可以得出结论：基于同组数据，同样的模型设置，不同模型的结果在模型精度以及空间表现上具有一定的差异，因此模型空间表现的不确定性无法避免。

图2 不同模型预测结果精度Figure 2 Accuracy of prediction results of different models

2.3 基于组合分布模型的厚朴适宜生境分布预测

由组合模型阈值二分类结果可知（图4-A），厚朴适宜生境面积较大，主要分布在陕西南部，甘肃东南南部以及河南部分地区的秦岭山地，还包括与秦岭地区毗邻的四川东部、重庆以及湖北西部、其分布区向东南延伸还包括贵州大部、广西北部、湖南西部、南部、江西大部、福建中部、北部，以及浙江南部。

当考虑植被类型限制的时，厚朴的适宜生境面积急剧减少，经过计算，我们得出适宜生境的面积为0.53×106km2。由图4-B 可知，厚朴适宜生境集中分布区主要位于秦岭地区以及四川西部，重庆市等湿润的温带山区，范围包括甘肃天水南部、陇南市南部；陕西汉中，安康以及商洛市；四川广元、巴中、绵阳、德阳、南充、达川、成都、广安、雅安、仁寿、宜宾等地区；重庆大部分地区。该区域应该作为厚朴原生产地保护以及人工种植的重点区域。

3 讨论与结论

图3 不同模型算法的结果Figure 3 Results of different model algorithms

成功物种分布预测一般需要如下条件：恰当数量的种群采样点数据；准确、可靠、能够描述适宜目标物种生境环境特征的环境因子数据，合理且能够描述物种生境特征的模型策略[31]。在物种分布研究中，最常用的物种分布数据来源于各种网络标本馆，这些数据容易获得且较为全面，但是此类数据往往缺乏精确的地理定位信息，同时分标本馆物种采样点数据采集自于二三十年前，其时效性难以保证，这些都会导致物种分布模型结果的不确定性。而本文采用了2000年之后具有精确地理定位信息的物种采样点数据，同时我们对原始数据进行了筛选，最大限度地避免了分布数据空间自相关的影响，保证了本文中参与建模的厚朴点位数据的准确性与代表性。植物生长和繁殖需要适宜的温度条件和充足的水分，因此环境因子的适宜性也是保证物种分布研究结果合理性的关键之一[9，18，31]。前人研究表明，19 个bio 生物气候变量对物种分布研究是有用的，并且是应用最为广泛的数据之一[31]。在本文中我们采用相关系数对19 个生物气候变量、3 个地形变量共22 个环境因子进行筛选，以避免因子相关性的影响。同时以上过程中我们发现年平均气温与海拔的相关性较高大于0.7，因此必须要对这两个因子进行筛选，这两个因子中海拔属于间接变量，此类变量与物种生存不具有直接生理相关性，只是通过影响其他因子间接的影响物种分布。对于海拔而言，他不直接作用于物种分布，而是通过影响区域水热的再分配来影响物种分布，在不同的时空尺度下，海拔与气温降水因子的关系不是恒定的，因此其与物种分布的关系也是变化的。因此我们选择能够更加直观表达区域热量情形的平均气温作为构建模型的候选因子之一[32]。目前已有大量的文献阐述了物种分布模型构建中不同模型算法的优劣，但是没有任何一种或者一类模型算法在所有物种分布研究中都能获得最高的统计精度，适合于所有生态位特征的最优模型算法并不存在[9]。组合物种分布模型提供了一个更加合理的模型建模方案。其理论基础是：单一模型算法不能在所有情况下都能够得到最优结果，所有的模型算法都有缺陷，但可以提供有用信息。大多数现有研究已经表明，通过集合模型，物种分布预测的准确性可以显著提高[33]。

图4 组合模型结果Figure 4 Results of Ensemble Mode

本研究综合考虑地形、气候和植被因素等环境条件对重要药材树种厚朴适宜生境的影响，研究结果表明影响厚朴适宜性分布的主要因子是温度年较差（Bio 7），年平均气温（Bio 1），年降水量（Bio 12）以及平均日较差（Bio 2），在本研究中根据环境因子响应曲线，我们可以得出厚朴适宜生境的Bio7 范围为 25～30 ℃，Bio1 的范围为 10～21 ℃，Bio12 的范围为 1 000～2 000 mm，Bio2 的范围为 Bio2＜9 ℃，这些因子范围描述了温暖湿润，且气温具有季节性变化特征的温带以及亚热带气候条件。前人研究表明，厚朴喜光，性喜凉爽、潮湿的气候，宜生于雾气重，相对湿度较大的地方。分布区年平均温度10～20 ℃，1月平均温度3～9 ℃，年降雨量800～1 800 mm 最优为1 400 mm[15-16]，本研究的有结果与前人研究相符，也说明了模型的合理性。

本研究基于BIOMOD 平台构建了9 种单物种分布模型，但是基于模型的统计学精度以及空间表现我们可以看出，模型的不确定性确实存在，统计学精度较高的模型空间表现不一定最好，如在本研究中的MARS 模型的统计学精度较高，但是空间表现都一般。可以确切的得出结论的是：传统的包络模型SRE 模型效果最差，模型结果无法给出分布细节，精度较低。以GLM 和CTA 为代表的传统回归模型和分类模型在厚朴适宜生境分布时效果一般，而以GBM、RF 为代表的复杂机器学习模型能带来更好的模拟效果。本研究再一次证实了不存在适应于所有建模条件的最优模型，集成多模型结果的组合物种分布模型能够在一定程度了减少模型的不确定性，增加模型的精度。