基于随机森林模型的国家重点保护药用植物物种优先保护区的识别*

张军， 石若利， 张雪， 谭雪平， 刘涛， 魏爱松

云南大学 建筑与规划学院， 昆明 650504

植物是人类赖以生存和发展的必要的物质基础， 人类在日常的衣食住行和医疗卫生等方面均离不开植物[1]． 药用植物生态资源种类是构成我国现代自然科学和社会土壤环境资源的重要组成部分． 研究表明， 我国具有悠久的多种药用植物资源种类和药用植物栽培技术． 《中国植物红皮书： 稀有濒危植物》第1册共收载植物354种， 药用植物168种， 其中稀有种38种， 渐危种84种， 46种药用植物面临灭绝的危险． 如何有效保护国家重点药用植物物种已成为全世界共同面临的环境问题之一．

我国在药材生产方面已经拥有了较高的价值和优势， 药材资源蕴藏量大， 但也有部分药材资源面临着市场供不应求的困难． 以我国6个区域进行排序， 能够供药用植物种类的分别是华东、 中南、 西南、 华北、 东北、 西北[2-3]． 目前我国境内有明确记载的药用植物共11 146种， 大多数传统中药材采用野生药用资源[4]， 属于濒危状态的植物已接近3 000种， 其中具有一定药用价值的植物约占60%～70%[5]． 目前， 《种子法》及各种相关中药法规在我国相继出台， 建立了对药用植物新品种的审定和认证机制， 并进一步出台了切实可行的法规制度以管理药用植物种子的生产经营等活动． 2013年， 中国提出了《生态保护红线管理办法》， 使生物多样性保护工作迈上了新的台阶． 生态保护红线是指对维护国家和区域生态安全及社会经济可持续发展具有重要战略意义的关键生态保护区域． 生态保护红线的划定能够维持关键物种、 生态系统与种质资源生存的最小面积， 有效保护了生物的多样性． 中国国家重点保护药用植物物种的优先保护区也是生态保护红线的重要组成部分． 如何准确识别和确定优先保护区是生物多样性保护研究的重要热点之一， 也是进行保护活动的前提条件． 准确识别和确定优先保护区能够指导资源的合理分配， 使保护效益最大化[6]． 当今世界， 对于生物多样性的优先保护区的确定研究， 在许多发达国家和地区已经进行了很多的实践和探讨， 并且取得了较为丰硕的成果．

准确有效地获取物种的适宜分布区对于识别和确定优先保护区至关重要， 物种分布预测技术和地理信息系统(GIS)为生物多样性优先保护区的识别提供了一种可靠的方法[7-8]． 该方法基于物种分布相关数据和环境因素， 利用物种分布预测模型来获取物种的适宜分布区域， 在全球范围内得到了广泛的应用． 钱灵颖等[9]以39种厦门市重点保护植物为研究对象， 通过物种分布模型MaxENT获得物种潜在的分布栅格图， 利用空间保护优先化定量工具 Zonation软件识别理论上既适宜重点保护植物生存又能够保证景观连通性的区域， 获得本地重点保护植物景观保护等级． 周来[10]利用林分数据和土壤剖面调查数据， 分析了杉木人工纯林和混交林的土壤理化性质特征， 基于修正的通用土壤流失方程对研究区的土壤保持情况进行了定量分析， 并结合采伐迹地植被恢复特征进行优先保护区域的识别． 王静[11]基于系统保护规划， 以海口市沼泽湿地优势植物潜在适生区为保护对象， 通过系统保护规划软件， 计算对海口市沼泽湿地保护具有不可替代价值的关键性区域， 筛选出最优保护目标， 进一步识别得到优先保护区． 陈龙等[12]以太行山生物多样性保护优先区北京区域为例， 提取待保护目标物种的分布范围， 识别目标物种集聚的热点区和空缺区． 陶国庆等[13]以滇西北三江并流区为例， 基于多准则决策分析方法(MCDA)， 对滇西北 24 种植被亚型进行保护价值评分同时分析了保护成本， 并对保护价值和保护成本进行叠加分析， 确定研究区植被的保护等级以及优先保护区． 周阳[14]以淇澳岛红树林保护区为对象， 通过研究该地区植物的物种优势度、 群落物种多样性、 种群空间等指标分析该地区的植物分布格局．

物种分布预测模型在对各生物物种保护区的识别应用上已十分广泛， 包括回归模型、 生态位模型和机器学习模型等， 但针对国家重点保护药用植物物种的优先保护区识别还很少见． 一些研究人员将现有模型与其他方法相比较， 认为随机森林模型在全球范围内做相应的物种分布研究工作所表现出的整体性能最好[15-16]．

本研究以126种国家重点保护药用植物物种为研究对象， 依托ArcGis 10.3， MATLAB等软件， 采用随机森林模型， 通过对药用植物物种进行建模， 预测各个物种的适应性生存率． 探讨以下问题： ① 以国家重点保护药用植物物种的生境性适宜指数作为本次研究的切入点， 进行空间自相关分析， 确定国家重点保护药用植物物种优先保护区的所在位置； ② 进行优先保护区空缺分析； ③ 针对未被保护的优先保护区进行其优先保护的顺序分层． 探讨此类问题能够对国家重点保护药用植物物种的发展提供稳定的空间， 也能为政府在制定管理决策时提供科学的依据， 让稀有的自然资源能够得到最大的保护[17]．

1 材料与方法

1.1 物种分布数据的获取及处理

药用植物物种类别主要结合《中国植物红皮书： 稀有濒危植物》[18]《Chinese traditional medicine resources》[19]以及全国中药资源普查成果进行选取． 《中国植物红皮书： 稀有濒危植物》中记录了354种濒危植物， 包含物种的历史分布地区、 受到威胁的程度、 在我国的生存和分布情况及其保护层次等相关信息． 药用植物物种分布点数据来源于中国科学院植物研究所(http： //www.iplant.cn/)， 因该网站的数据只有热力图， 故该文所使用的数据是热力图经矢量化处理后的数据， 并根据其中的重要性、 特殊性以及数据的相对完整性进行选取， 同时剔除了在分布地点上数据不全的各类濒危药用植物野生物种， 最终选择了作为国家重点保护的药用植物126种， 其中39种为国家一级保护植物， 37种为国家二级保护植物． 在本文选取的126种药用植物物种中， 被IUCN《受威胁物种红色名录》评定为易危(UV)物种的有 11 种， 濒危(EN)物种的有 76种， 极危(CR)物种的有1种， 受到不同程度威胁的物种占所选物种的69.84%．

1.2 环境变量及相关数据出处

环境变量是随机森林模型能够顺利完成模型的建立及结果预测的基本条件． 本研究选取了12个能够反应气候特征、 栖息地和人类影响的环境变量， 以此来预测物种的适宜分布区(表1， 图1和图2)．

表1 环境变量的出处及含义

1.3 模型与分析

随机森林模型是Breiman在2001年提出的一种分类预测方法， 是一种汇总式自学习数据处理方法、 现代分类与回归技术[20]． 随机森林模型是基于分类树算法进行模拟和迭代， 在变量和数据使用上进行随机化， 分成若干分类树， 再汇总分类树的结果进行判别[21-22]． 随机森林模型对多元共线性不敏感， 但对大量数据的处理效率较高， 可以预测多达几千个的解释变量， 在数据缺失时， 仍然有很高的精度． 本研究运用MATLAB作为模型运行平台完成相应的运算．

1.3.1 建构随机森林模型

以中国区域作为分析对象， 在研究中选取公里网格作为预测背景点， 点与点之间经度和纬度的间隔均为0.2， 此研究共有24 050个背景点被选择． 运用ArcGis 10.3软件提取各物种相应的分布点和背景点及其所有与环境变量相关的值， 随机生成3 000个不适合物种生长的背景点数据， 用训练集和测试集来总结物种的相关分布数据和背景点数据． 训练集由分布数据构成， 其数据分别占相应基数的75%， 剩余的25%则构成测试集． 训练集的作用是构建模型， 测试集的作用是检验模型的精确度[23-24]． 建模后， 对每个物种进行预测， 得到背景点物种的适生概率， 将每个背景点上所有物种的适生概率进行叠加并进行标准化， 得到国家重点保护药用植物物种的生境适宜性指数．

1.3.2 模型精度检验

设Ntree为随机森林模型中所有分类树的数目， 当Ntree较小时， 随机森林模型的分类错误率较高， 而性能较低． 随机森林模型的特点是它具有不过度拟合的性质， 所以它可以通过增加Ntree的数量， 来保证集成化分类器的种群多样性． 然而在构建随机森林模型时的操作复杂程度却与Ntree成正比，Ntree过高， 随机森林模型在构建过程及运算时所耗费的时间就会大大增加[25]． 因此，Ntree对随机森林模型的性能、 可解释性和复杂性之间的平衡具有重要意义． 本文通过测试集对Ntree进行选择， 当Ntree为1 500颗时， 测试集达到最大正确率84.42%， 说明模型预测的准确性较高．

审图号： GS(2019)1815号图1 生物气候变量情况

1.3.3 最终适宜度的计算

在生物气候变量(6个)和海拔数据的基础上进行所有数据标准化处理， 使用随机森林模型计算出背景点的基础适宜度， 插值生成1 km大小的栅格， 再根据人类足迹指数、 生态功能区、 植物净初级生产力、 土地利用类型计算出最终适宜度指数(Fls)：

Fls=(Fs-0.5×Footprint)×Efz×Npp×Landuse

(1)

式中，Fs为基础适宜度，Footprint为人类足迹指数，Efz为生态功能区，Landuse为生物群落区，Npp为植物净初级生产力．

结果以0.25，0.50，0.75为节点， 分为不适宜、 次适宜、 适宜和极佳生境． 去掉不适宜、 次适宜生境， 将适宜、 极佳生境视为优先保护区， 分别与自然保护区、 生态功能区叠加， 观察物种的受保护情况．

审图号： GS(2019)1815号图2 其余重要变量情况

1.3.4 优先保护区的识别与空缺分析

通过ArcGis 10.3平台结合反距离加权插值运算的方法计算生态环境的适宜性系数， 制作分辨率为1 km的生境适宜性指数图． 将生境适宜性指数图每一个像元都转换成为1 km×1 km的矢量网格， 1个网格代表1个评价单元．

本研究用GeoDa统计软件中的空间自相关分析方法探测各个生境适宜性指数的空间聚集区情况， 从而得到国家重点保护药用植物物种的优先保护区， 以确保药用植物物种能够长久地生存． 同时， 在ArcGIS 10.3软件中将识别出需优先保护的区域完成叠加操作， 叠加的对象分别是自然保护区和生态功能区． 自然保护区和生态功能区数据均来自中科院资源环境科学与数据中心(http： //www.resdc.cn/)． 为了降低空间尺度对空间自相关分析的影响， 根据近年来中国主要药用植物的具体地理学区域进行划分， 将我国重点保护药用植物物种生境适宜性指数图进行分类， 可分为岭南(平原)区、 华东区、 华中(平原)区、 华北(平原)区、 东北(平原)区、 华西区、 横断(太行)山区、 云贵区、 青藏(高原)区、 蒙新区， 共10个区作为中国主要药用植物的具体地理学区域．

全局和局部两个指标对于空间自相关性的分析能够起到很好的衡量作用． 首先利用全局指标对10个不同地理区域的生境适宜性指数的空间模型进行探测， 确定整个区域之间是否存在空间的自相关性和显著程度； 然后将10个不同的地理区域进行局部空间自相关分析， 进一步检验是否存在观测值的高值或低值局部空间聚集， 并识别高值区和低值区． Moran’sI指数通常被广泛应用在衡量各个空间元素的相互关系上， 该指数的大小范围在-1至1之间． 具有强烈空间正相关的要素用1来表示， 具有强烈空间负相关的要素用-1来表示， 不具有或者不存在空间相关性的则用0来表示[26]．

全局Moran’sI的计算公式：

(2)

局部Moran’sIi的计算公式：

(3)

2 结果与分析

2.1 识别优先保护区

通过对全局Moran’sI的计算公式进行综合计量， 获得10个地理分布区的生境适宜性指数(表2)． 结果显示， 各地理分布区内生境适宜性指数的全局空间自相关指数皆大于0.90(p<0.05)． 说明10个地理分布中评价单元的生境适宜性指数存在显著的空间自相关．

表2 10个药用植物地理分布区的全局Moran’s I指数

基于局部 Moran’sIi的计算公式进行综合计量， 获得10个地理分布区的生境适宜性指数的局部空间聚集情况． 再通过ArcGis 10.3软件对其进行描述， 完成一个具有较高值限的聚集区图像， 得出国家重点保护药用植物物种不适宜、 次适宜、 适宜和极佳生境的分布情况(图3)． 优先保护区是指定的适宜度指数值相对较大的保护区域， 该区域的特征主要是国家重点保护植物物种数量相对较多并且受到人类活动的干扰相对较少．

审图号： GS(2019)1815号图3 不同适宜性生境的分布情况

国家重点保护药用植物物种的优先保护区主要分布在我国的华西区、 华中区、 岭南区、 华东区、 云贵区、 横断山区， 其中包括长江以南的武夷山脉、 秦岭地区以及秦岭西北部至大巴山余脉山区西部、 长白山、 雅鲁藏布江下游河谷区、 滇西南、 海南岛中部． 此外， 西北部大多为不适宜区域， 次宜适应区域遍及云贵高原腹地． 通过计算， 优先保护区面积为78.85万km2， 约占我国国土面积的8.21%． 经统计得出优先保护区的斑块个数为3.40万个， 平均斑块面积为23.19 km2．

2.2 优先保护区空缺分析

将优先保护区与我国现有的自然保护区和生态功能区的面状数据进行叠加， 得到优先保护区目前的被保护情况(图4)． 优先保护区在自然保护区中的面积约为3.36万km2， 约占自然保护区总生境的3.16%， 占优先保护区总生境的4.26%； 优先保护区在生态功能区中的面积约为24.33万km2， 约占生态保护区总生境的10.78%， 占优先保护区总生境的30.86%． 综合考虑来看， 优先保护区被规划保护的土地总面积约为26.89万km2， 占我国优先保护区总生境的34.10%， 被规划保护的优先保护区总面积占比并不高， 说明我国对药用植物物种的保护还不够充分． 自然保护区以青藏区较多， 岭南区、 华东区较少， 而优先保护区主要分布在长江以南的岭南区、 华东区． 针对雅鲁藏布江下游河谷区、 滇西南、 东南丘陵区， 可适当增强保护区建设力度， 通过制定相关法规， 对国家重点保护药用植物物种适宜生长的地区加以保护．

审图号： GS(2019)1815号图4 自然保护区、 生态功能区和优先保护区叠加

2.3 优先保护顺序分析

对未被保护的优先保护区进行优先保护分层， 首先， 最应被优先保护的区域是云贵区、 岭南区、 华中区、 华东区、 横断山区， 其中包括武夷山脉、 南岭周边区域、 西双版纳、 海南岛中部等． 主要是因为该区域的药用植物相对较丰富， 生境适宜性普遍较高， 植被类型较复杂， 面积相对较小．

其次， 应该被保护的区域是东北区、 华北区、 华西区、 青藏区、 蒙新区， 其中包括长白山地区等， 主要原因是该区域内的优先保护区占地面积很小、 分布比较散、 覆盖率较低．

3 讨论与结论

国家重点保护药用植物物种的优先保护区在南部山脉地区的斑块所占面积较大、 分布也较集中． 这可能是由于该地区具备气候适宜， 有着大部分药用植物生长发育的光照与温度等条件， 有利于药用植物生长． 另外， 山区受人类干预的强度也相对较低， 因此这些地区的生物多样性极为丰富， 拥有大量特殊的植物物种． 根据《Flora of China》及相关统计资料表明， 省级以上药用植物的特有生态品种在我国西部呈现较为集中． 第一是与西南的一些省份有着紧密联系的地区， 如云南东南部、 南部和西北部， 四川西部(西北部、 西南部) 至中部， 西藏东部和东南部等； 第二是我国的东南部和地处华南的地区， 包括海南等； 第三是从华中到秦岭地区， 包括湖北西部至陕西南部等． 本文的研究结果表明， 目前国家重点保护药用植物物种的优先保护区和3个省级以上药用植物重点地区在中国境内存在部分重叠． 但由于这些地区大多位于一些经济发达地区， 受到现代人类的干扰大； 这些地区也是鸟类迁徙的主要连接地区， 对于迁移鸟类完整的自然生活历程具有重要的意义． 另外还存在东北的大、 小兴安岭， 东北至华南沿海及我国长江流域中下游等地区分布相对较少且分散、 破碎化的程度严重等状况． 导致这种特殊情况的原因主要是这些地区的自然海拔度和温度相对变化小， 不具备形成多样的地形及气候条件． 因此， 针对这些优先级保护区加强环境生态保护工作， 对于我国自然资源的生物多样性保护意义重大．

通过对国家重点药用植物物种优先保护区与国家生态功能区及自然保护区进行对比分析， 发现他们许多地方还存在着重叠， 特别是在生物多样性相对较高的横断山区、 云贵区、 岭南区、 华西区和华中区等， 说明这些地区在其生态机制功能的发挥以及对其生物多样性的维护等方面具有很大的保护价值． 由于划定目的和方法的不同， 优先保护区与重点生态功能区及全国生态功能区虽有重叠， 但重叠的范围并不大． 以我国南岭地区为例， 南岭北部山脉不仅被认为是一个属于天然生态水源综合涵养的功能区， 还属于生态功能区， 占生态功能区与自然保护区区域总面积的7.59万km2与1.23万km2， 相应的优先保护区占生态功能区和自然保护区的面积各为0.36万km2与0.83万km2．

此外， 我国中东、 中西部也是保护需求量大的地方， 而药用植物物种优先保护区在生态功能区、 自然保护区的覆盖范围较小． 长江以南地区， 特别是武夷山脉地区， 东南丘陵区的优先保护区范围较大， 应适当增加对这些区域的保护力度．

同时， 中国国家重点保护药用植物物种的优先保护区也不是一成不变的， 本文共选取126种药用植物物种， 其中在优先保护区范围内的药用植物物种有124种， 保护率高达98.41%． 对于未被优先保护区保护的物种以及由于缺少数据未被纳入到分析中的物种， 可以结合这些物种的分布情况， 对优先保护区作适当的补充和调整． 除此之外， 还需要针对影响这些物种分布的环境因素进行调查和更新， 如气候变化、 人类活动干扰、 社会经济发展水平的变化等， 才能更加真实地对优先保护区进行稳定且长久的保护和管理．

本研究选取随机森林模型预测126种国家重点保护药用植物物种的生境适宜率， 将适宜、 极佳生境作为国家重点保护药用植物物种的优先保护区， 分别与自然保护区、 生态功能区进行叠加计算得到国家重点保护药用植物物种的生境适宜度指数， 采用空间自相关的研究方法确定国家重点保护药用植物物种的优先保护区， 进一步明确了在保护工作中亟需重视的区域， 为相关政府部门更有效地分配和利用有限的保护资源提供科学的依据．