基于物种分布模型的画稿溪国家级自然保护区桫椤保护现状及影响因素*

许斌，钟悦，张焱

(1.内江师范学院 地理与资源科学学院,四川 内江 641100;2.青海师范大学 地理科学学院,青海省自然地理与环境过程重点实验室,青海 西宁 810000;3.叙永县画稿溪国家级自然保护区管理处,四川 泸州 646400)

自然保护区在保护生物多样性、自然景观及自然遗迹，维护国家和区域生态安全，保障我国经济社会可持续发展等方面发挥着重要的作用，是生态文明建设的重要载体，是保护生物多样性最有效的方式之一[1-2]。除典型的森林生态系统外，国家重点保护的野生动、植物物种及其栖息地也可能是森林类保护区的重点保护对象。四川省泸州市叙永县画稿溪于2003年被列为国家级自然保护区，以亚热带原始常绿阔叶林生态系统和桫椤群落及其伴生的珍稀野生动植物为主要保护对象。

桫椤(Alsophilaspinulosa)是隶属于桫椤科(Cyatheaceae)桫椤属(AlsophilaR.Br.)的树形蕨类、第四纪冰川劫难中幸存的孑遗植物，对研究古植被演化和地质变迁等具有重要参考价值[3]。桫椤喜暖湿环境，多见于我国南方省区[4-5]。已有关于桫椤的研究多见于大尺度分布[4-5]或基于基因水平分析其遗传信息[6]等，关于小尺度自然保护区桫椤生境的研究较少[3,7]。画稿溪国家级自然保护区内桫椤种群分布相对集中，本文以保护区内的桫椤分布资料，运用MaxEnt模型预测桫椤的潜在适生区、结合地理探测器分析影响桫椤生栖地因素，分析保护空缺和人类活动对桫椤潜在适生区的影响，为类似的动、植物保护区生物多样性监测和保护提供参考。

1 材料与方法

1.1 研究区概况

画稿溪国家级自然保护区(105°03′～105°40′E、27°42′～28°31′N)位于四川省泸州市叙永县境内，东、北面与贵州省赤水市、四川省古蔺县接壤，面积23 454 hm2，海拔范围320～1 780 m。其地形处于云贵高原黔山地北缘与四川盆地中部低山丘陵的过渡地带，地势南高北低，有画稿溪、墩梓河等河流由南向北流经保护区；其地处亚热带湿润季风气候区，四季分明、热量充足、降水丰沛、雨热同季、日照偏少、无霜期长，可分为干热河谷(海拔800 m以下)、暖湿平丘(海拔270～500 m)、暖湿低山(海拔500～1 000 m)和冷湿中山(海拔1 000～2 000 m)4个气候类型；年平均气温17.3 ℃，最热月(7月)均温27 ℃，最冷月(1月)均温6.9 ℃，≥10 ℃积温5 802 ℃；年平均降雨量1 285.6 mm，年平均相对湿度80%；年平均日照1 170 h以下，日照不足27%；无霜期350 d左右。土壤以红冷沙土、黄冷沙土、豆面土和灰包土4个土种为主。其植物区系处于泛北极植物区和古热带植物区的汇集和分界地带，植物多样性较高，种子植物中属热带成分295属，温带成分219属。

1.2 数据来源与预处理

2021年7月22—25日开展保护区内桫椤分布数据调查。由于调查期间泸州市普降大雨，墩梓河水大涨而致公路阻断，因此墩梓河沿岸自水尾镇官斗村村委会驻地上游部分、以及灵官梯河谷的桫椤分布采用鞠文彬等[7]的调查数据，其余桫椤分布数据通过实地调查获取。气候数据(月均温、月降水量、月光照时间、月平均相对湿度)来源于泸州市统计年鉴(http://tjj.luzhou.gov.cn/)，包含泸州市5个气象站点2011—2020年统计值。高程数据(DEM，数字高程模型)空间分辨率为30 m，来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)；泸州市行政区划数据地图资料来源于国家基础地理信息中心(http://www.ngcc.cn)；画稿溪保护区边界、土地利用、植被等资料来源于画稿溪国家级自然保护区管理处。

气候数据(月均温、月降水量、月光照时间、月平均相对湿度)通过ArcGIS 10.2 软件克里金插值为30 m空间分辨率的栅格数据；画稿溪保护区与赤水保护区毗邻，处于同一个植被区[7]，因此参照赤水保护区内桫椤物候特征(3月上旬叶芽开始抽发、5月发叶量出现第1个高峰期、6—8月发叶量降低、9—10月发叶量出现第2个高峰期，11月至第2年3月桫椤处于休眠或半休眠状态[8])对4—6月、7—8月、9—10月、11—3月、整年(1—12月逐月)的气候数据进行处理——温度、光照时间、相对湿度求取均值，降水量进行累积求和得到64个气候数据。方圆90 m内如有多个桫椤分布点，则随机选择其中1个保留以降低样本点空间上的自相关，最终得到桫椤分布点210个(图1)。通过ArcGIS 10.2 软件，由DEM生成坡度数据，参照文献[9]获取样本点采样偏差数据。选用ArcGIS 10.2 软件自然断点法将环境变量重分类为4类、转换成点后提取分类后的值，用于地理探测器分析；将环境变量转换成.asc格式用于MaxEnt模型预测。

为避免气候变量过多而致模型过度拟合，随机选择桫椤分布数据作训练集与测试集并交叉验证。以桫椤分布数据、64个气候数据、高程、坡度为输入环境变量，运行最大信息熵模型(MaxEnt，version 3.4.1)3次，检验各气候变量对桫椤潜在适生区的预测贡献。运用ArcGIS软件中band collection statistics工具对气候变量进行相关性检验，剔除相关系数(R)高于0.8的2个因子中对 MaxEnt 模型预测结果贡献较小的变量[10]，最后筛选出适合画稿溪桫椤生长的5个环境变量：海拔(DEM)、9月降水量(prep_9)、10月降水量(prep_10)、11月—翌年3月降水量(prep_11-3)、11月光照(sun_11)。

1.3 桫椤的潜在适生区预测

MaxEnt模型是一种基于最大熵理论提出的生态位模型[11]，由已知物种分布点的经度、纬度信息进行建模，计算物种分布栅格环境变量约束条件，在研究区所有栅格内寻找该约束条件下最大熵的可能分布，在物种分布点数量不定和环境变量之间相关性不明确的情况下仍具有较好的预测效果[12]，预测物种潜在分布结果较好[13]，能够对物种未来潜在地理分布进行定量分析[10,14-16]。采用随机选择桫椤分布数据作训练集与测试集，辅以样本点采样偏差数据，勾选软件自带交叉验证功能进行验证，以MaxEnt模型10次计算结果平均值作为桫椤的潜在适生区。采用受试者工作特征曲线(receiver operating characteristic curve，ROC)下的面积值(area under curve，AUC)检验 MaxEnt 模型预测精度：AUC值在0.5～0.6即模型不合格；0.6～0.7较差；0.7～0.8一般；0.8～0.9良好；0.9～1.0优秀[10,17]。AUC值由MaxEnt模型输出。训练集、验证集AUC值为0.981，表明模型预测水平可信(图2)，应用MaxENT模型预测桫椤潜在适生区分布可行。MaxEnt 模型预测桫椤潜在适生概率(P)取值为0.0～1.0，使用最大训练灵敏度加特异性(MSS)逻辑阈值区分模型中合适和不合适的生境，运用ArcGIS软件依据MSS判定的合适和不合适阈值将预测结果划分为5个等级，即：适生概率P<0.1为非适生区(I)；0.10.7为高适生区(V)[9]。

图2 受试者工作特征曲线Fig.2 ROC curve generated by MaxEnt model

1.4 影响桫椤潜在适生区因素分析

地理探测器(geographical detector)是研究地理现象空间分异性，并揭示其背后驱动因素的一种统计学方法。其核心思想是：假设研究区分为若干子区域，如果子区域的方差之和小于区域总方差，则存在空间分异性；如果两变量的空间分布趋于一致，则两者存在统计关联性[18]。该方法在自然科学[19]和人文社会科学[20]领域均得到广泛应用。因子探测器用于分析因变量Y的空间分异性，并探测某个因变量X对因变量Y的空间分布的解释力，用q值度量[21]；交互探测器可识别不同因子之间的交互作用，即评估多因子共同作用时是否会增加或减弱对因变量的解释力，或这些因子对因变量的影响是相互独立的，交互方式如表1；风险探测器用于判断两个子区域间的属性均值是否有显著的差别，用t统计量来检验；生态探测用于比较两因子对因变量的空间分布的影响是否有显著的差异，以F统计量来衡量[18]。

表1 地理探测器交互作用方式Tab.1 Interaction mode of geographical detector

在MaxEnt模型计算结果中，选取对画稿溪桫椤潜在适生区贡献率和重要性均不为0的环境变量，包括：海拔(DEM)、1月降水量(p_1)、9月降水量(p_9)、10月降水量(p_10)、11月至翌年3月累积降水量(p_11-3)、2月光照(sun_2)、5月光照(sun_5)、10月光照(sun_10)、11月光照(sun_11)、11月至翌年3月累积光照时长(sun_11-3)、7月平均温度(temp_7)、12月平均温度(temp_12)共12个环境变量，在ArcGIS 10.2中重采样为90 m、采用自然断点法进行重分类为4类作为地理探测器的输入变量，分析桫椤潜在适生区影响因素。

1.5 画稿溪桫椤潜在适生区人为干扰分析

画稿溪保护区内有县级以下公路、居民点和耕地。因耕地多在居民点周围，因此本研究参照已有研究结果[22-23]，借助ArcGIS 10.2软件缓冲区分析工具，分别以公路、耕地建立100 m缓冲区，进行合并，整合合并后的零星图斑，分别以公路缓冲区、耕地缓冲区、耕地公路缓冲区并集划分人为干扰区域，以衡量桫椤潜在适生区受干扰情况。

2 结果与分析

2.1 模型预测精度及关键环境变量

各变量对预测桫椤潜在适生区空间分布贡献排序依次为：DEM>9月降水量>11月光照>10月降水量>11月至翌年3月累积降水量。桫椤适生概率高于0.5时对环境需求：海拔(DEM)低于710 m，9月降水量(p_9)低于108 mm，11月光照(sun_11)在50.8～52.6 h之间、10月降水量(p_10)88～99 mm，11月至翌年3月累积降水量(p_11-3)165～174 mm(图3，表2)。环境变量重要性由MaxEnt模型输出。

图3 环境变量重要性刀切法检验Fig.3 Jackknife test of environment variable importance

表2 环境变量范围及贡献率Tab.2 Range and contribution of environmental variables

2.2 画稿溪桫椤潜在适生区分布

画稿溪桫椤II—V级潜在适生区沿墩梓河、画稿溪分布，其中：V级、IV级适生区主要沿墩梓河、画稿溪两岸呈狭长分布，面积分别为53.28 hm2和262.53 hm2；边缘适生区(III)与低适生(II)区面积分别为427.32 hm2和1 217.79 hm2(图4、表3)。

图4 画稿溪桫椤潜在适生区分布Fig.4 The potential distribution of A.spinulosa in Huagaoxi NR

表3 画稿溪桫椤潜在适生区面积Tab.3 Potential distribution area of A.spinulosa in Huagaoxi Natural Reserve

从垂直方向看，画稿溪桫椤适生程度随海拔升高而降低，II级潜在适生区总体分布在海拔约650～850 m范围内、III级潜在适生区分布在约560～680 m区域、IV级潜在适生区总体在约490～580 m范围、V级潜在适生区总体分布在约440～520 m区间，不同等级潜在适生区空间分布垂直差异明显(图5)。

图5 画稿溪桫椤潜在适生区海拔分布箱式图Fig.5 Boxplot of the altitude of A.spinulosa potential distribution in Huagaoxi Natural Reserve

2.3 画稿溪桫椤潜在适生区与保护区功能区空间的关系

画稿溪核心区、缓冲区、实验区II—V级桫椤潜在适生区分布总面积分别为89.49、178.71、1 692.71 hm2，占比分别为4.56%、9.11%、86.32%。其中：适生区和高适生区(IV、V级潜在适生区，其适生概率P>0.5)在核心区、缓冲区、实验区的面积分别为4.81、17.24、293.75 hm2，占适生区和高适生区总面积比例分别为1.52%、5.46%、93.02%；边缘适生区在核心区、缓冲区、实验区的面积分别为19.41、23.28、384.62 hm2，占比分别为4.54%、5.45%、90.01%；低适生区在核心区、缓冲区、实验区的面积分别为65.26、138.19、1 014.34 hm2，占比分别为5.36%、11.35%、83.29%(表3)；画稿溪桫椤潜在适生区约86.32%面积分布在实验区，缓冲区和核心区分布偏少，IV、V级潜在适生区在核心区分布仅为1.52%。

2.4 画稿溪桫椤潜在适生区与人类活动空间的关系

在保护区公路、耕地0.1 km缓冲区内，分布有大多数桫椤潜在适生区(图6a)。公路缓冲区内桫椤高适生区、适生区、边缘适生区、低适生区面积分别为43.49、168.65、121.14、233.04 hm2，占各等级适生区的面积比例则是81.63%、64.24%、28.35%、19.14%；其中适生区和高适生区面积为212.15 hm2，占适生区和高适生区总面积315.81 hm2的67.18%。耕地缓冲区内适生区面积分别为0.97、2.95、101.00、515.54 hm2，占各等级适生区的面积比例则是1.81%、1.12%、23.63%、42.33%；其中适生区和高适生区的面积为适生区和高适生区总面积的1.24%。公路和耕地缓冲区并集内桫椤高适生区、适生区、边缘适生区、低适生区面积分别为36.99、169.10、189.71、582.33 hm2，占各等级适生区的面积比例则是69.43%、64.41%、44.39%、47.82%；公路和耕地缓冲区并集内适生区和高适生区的面积为206.09 hm2，占适生区和高适生区总面积的65.26%。在保护区，公路主要影响桫椤高适生区、适生区的分布，耕地则限制边缘适生区分布。

画稿溪不同植被类型空间内，不同等级桫椤潜在适生区分布比例差异较大(图6b、表4)。竹林分布区域桫椤潜在适生区面积最大，其次为常绿阔叶林、水田、杉木林、常绿与落叶阔叶混交林，再次为山地灌丛和耕地。竹林主要分布于画稿溪和敦梓河两侧区域，竹林内高适生区、适生区、边缘适生区、低适生区面积分别为35.37、154.26、212.04、582.57 hm2，占各等级潜在适生区面积比例分别为66.39%、58.76%、49.62%、47.84%。在人类活动较为强烈的竹林、水田、坡耕地分布区域内，桫椤潜在高适生区、适生区、边缘适生区、低适生区面积总计为35.37、156.78、243.54、779.22 hm2，占各等级潜在适生区面积比例分别为66.39%、59.72%、56.99%、63.99%。总体来看，画稿溪桫椤潜在适生区与人类生产生活活动及经济来源所依赖的耕地、竹林等空间冲突较为明显，受人类活动干扰较为突出。

图6 画稿溪桫椤潜在适生区与人类活动空间的关系Fig.6 Spatial relationship between potential suitable area of A.spinulosa and human activities in Huagaoxi Natural Reserve

2.5 画稿溪桫椤潜在适生区影响因素分析

因子探测结果显示，不同环境变量对画稿溪桫椤潜在适生区分布影响的解释力具有差异，其中解释力最大的环境变量是海拔为22.4%，其余环境变量解释力在2%～3.4%之间，最小为10月份平均温度(2%)，这一结果说明除海拔外，其余环境变量对画稿溪桫椤潜在适生区分布影响较小(表5)。进一步对海拔变量进行风险探测分析，结果显示不同海拔区间内桫椤适生概率存在显著差异，且通过显著性水平为0.05的t检验(表6)。这亦表明海拔变量不仅控制画稿溪桫椤潜在适生区的分布，且在不同海拔梯度的影响具有差异性。

表5 画稿溪桫椤潜在适生区影响因素P值和q值Tab.5 P and q value of A.spinulosa potential suitable habitat influencing factors in Huagaoxi Natural Reserve

表6 画稿溪桫椤潜在适生区海拔因子差异性Tab.6 Elevation difference of A.spinulosa potential suitable habitat in Huagaoxi Natural Reserve

地理交互探测结果(表7)表明任意两两因子叠加后解释水平均高于单一因子，其中DEM与9月份降水变量叠加后对桫椤分布的解释水平最高(0.462)，其次为DEM与2月份光照变量叠加的解释水平(0.363)、DEM与5月份光照变量叠加的解释水平(0.357)、DEM与11月至翌年3月份累积光照时长叠加的解释水平(0.347)。这表明海拔和9月份降水量对画稿溪桫椤潜在适生区分布的影响最为明显，海拔与光照变量叠加亦对画稿溪桫椤潜在适生区空间分布具有较强的解释能力。

表7 画稿溪桫椤潜在适生区影响因素交互作用下的q值Tab.7 q value of influencing factor interactive of A.spinulosa potential suitable habitat in Huagaoxi Natural Reserve

3 讨论与结论

3.1 讨论

(1)自然因素对画稿溪桫椤潜在适生区分布的影响 MaxEnt模型预测结果显示画稿溪桫椤潜在适生区主要沿河谷分布，最为明显的影响因素为海拔，海拔与9月份降水量、11月份光照等变量累积贡献率达91.5%；地理探测器分析结果亦表明，海拔对桫椤潜在适生区的解释力最大，与9月份降水变量、2月、5月、11月光照变量分别叠加后的解释力亦较大；均表明海拔因素是限制画稿溪桫椤潜在适生区分布的主要影响因素，MaxENT模型与地理探测器分析结果具有较好的一致性。

桫椤III—V级潜在适生区大部分位于保护区内海拔约680 m以下河谷地带，源于河流区域的温热湿润小气候与桫椤生长所需的生长环境相一致[3,24-25]。地理探测器分析结果表明不同垂直梯度上桫椤潜在适生区分布差异明显；实地调查发现，桫椤沿河谷分布，其垂直分布范围距离谷底100～150 m，与不同海拔梯度桫椤潜在适生区分布具有差异性相一致。

(2)人类活动对画稿溪桫椤潜在适生区的影响 画稿溪保护区成立之前该区域有大量村民居住，致使大量原生常绿阔叶林和竹林遭受不同程度破坏，现天然常绿阔叶林仅分布于交通不便的陡坡地带；地势平缓、靠近公路和居民点的植被已被各类人工林如大叶慈竹(Dendrocalamusfarinosus)林、慈竹(Bambusaemeiensis)林、杉木(Cunninghamialanceolate)林以及旱地、水田、果园等代替[7]。2020年叙永县水尾镇竹产品(竹笋干)1 510 t，占全县竹笋干3 215 t的46.97%，竹笋干产量充分说明竹产业为水尾镇居民收入的主要来源之一。沿画稿溪、墩梓河开凿的2条公路是保护区内外联系的主要通道，河流、公路两侧分布大量竹林，桫椤亦沿该两条河流、公路分布(图6b)，在公路影响范围内，约有67.18%的适生区和高适生区分布，高适生区达81.63%，且公路及其缓冲区多处于实验区内，科学试验、旅游、驯化繁殖野生动植物以及生产活动，与桫椤高适生区在空间上高度重叠，均对桫椤生栖地产生负面影响。

画稿溪桫椤主要分布在乔木的中下层，林下具有一定数量的幼苗库[26]，沿河岸分布有人工竹林、果园等，干扰物种通过其根系对周边的桫椤产生有消极影响的化学物质来抑制桫椤根系生长，威胁桫椤正常生长发育甚至使桫椤立木度减少[27]，如：毛竹(Phyllostachysheterocyclecv.pubescens)根系在垂直方向上随土层深度的增加对桫椤总根系生物量影响逐渐减少，水平方向上，桫椤根系随离桫椤距离增加而减少，在离桫椤60 cm处最低，桫椤细根比根长随干扰增强而下降[28]；画稿溪桫椤适生区范围的高密度人工竹林在一定程度上限制了桫椤种群的生长空间，林分干扰显著影响桫椤种群幼苗以上各阶段的密度[7]，一定程度上造成桫椤群落生境狭窄。

(3)画稿溪自然保护区桫椤保护建议 MaxEnt模型基于物种对基础生态位的需求进行适生区预测，侧重于揭示物种潜在地理分布及气候特征[29]。画稿溪主要保护对象为亚热带原始山地常绿阔叶林生态系统，桫椤是其较为重要的种群之一。保护区内桫椤种群年龄结构呈金字塔型，属稳定的森林群落，适生区域与人类活动空间重叠，受保护力度有待增强。综合已有研究结果[7,26]，本文提出如下建议：①充分利用现有画稿溪管理委员会做好宣传与教育工作，在公路沿线适当增加物种保护和生态价值宣传栏；②约有189.63 hm2高适生区和适生区与竹林空间重叠、占高适生区和适生区面积60.05%，有67.18%的适生区和高适生区分布于公路周边100 m范围，宜对公路沿线人工林进行适当干预，降低人工竹林对桫椤生存空间的侵占强度；③有23.63%的边缘适生区与耕地缓冲区重合，严防居民用地、建筑用地及农用地扩张；④宜对强烈干扰桫椤种群的林木进行适度砍伐，以保障桫椤幼苗的成活率或其他龄级桫椤生长；定期对已建立的桫椤监测样地进行复查，并监测生物多样性，杜绝外来入侵物种；⑤适时提高桫椤潜在分布区对应的保护区等级，如将部分桫椤潜在适生区域的功能由缓冲区划入核心区，实验区划为缓冲区甚至核心区等。

本文考虑到桫椤物候、降水、温度、地形以及人类活动等因素对其潜在适生区分布的影响，今后拟结合画稿溪保护区植物资源特性进行植物物种调查，以期能更好地探究画稿溪桫椤适生空间，分析保护区范围与画稿溪实际生物多样性热点分布空间关系，提高保护效率。

3.2 结论

本文利用2011—2020年泸州市气象数据、桫椤分布采样信息和数字地形产品，采用MaxEnt模型预测画稿溪国家自然保护区内桫椤潜在适生区、运用地理探测器分析适生区影响因素，探讨适生区分布空间、与人类活动空间关系，获得以下结论。

(1)海拔及因地形形成的小气候控制桫椤的分布空间，海拔是限制画稿溪桫椤分布最大的因素，海拔与9月份降水变量叠加对研究区内桫椤潜在适生区空间分布的解释水平最高。

(2)画稿溪桫椤潜在适生区沿墩梓河和画稿溪两岸呈狭长分布且海拔低于680 m，约86.52%的潜在适生区分布在实验区，IV、V级潜在适生区在核心区分布仅为1.53%，对桫椤的保护存在一定空缺。

(3)保护区竹类植物根系对桫椤根系生长具有限制作用，竹林主要分布于保护区内公路沿线，现有桫椤植株亦大部分位于竹林内，竹林分布区与桫椤潜在适生区和高适生区在空间上亦相互重叠。竹林经营是影响桫椤潜在分布最激烈的人类活动，耕地垦殖则限制了桫椤边缘适生区的扩展。

(4)了解桫椤生长特性，加强人为干预以确保桫椤生存空间，提高画稿溪桫椤潜在适生区分布空间对应的保护区功能等级，如将部分桫椤潜在适生区域的功能由缓冲区划入核心区、实验区划为缓冲区甚至核心区，可能更有利于画稿溪桫椤种群的保护。