龙溪-虹口国家级自然保护区川金丝猴栖息地现状及震后影响
2018-04-09雷雨，青菁，何可
雷 雨，青 菁，何 可
(1.甘孜卫生学校，四川 甘孜 626001；2.西华师范大学 生命科学学院，四川 南充 637009)
在物种保护工作中，物种的分布区信息对于制定保护措施有着极为重要的价值，在保护决策过程中有着关键的作用[1-2]。同样，在控制入侵物种的生态管理工作中，预测物种的潜在分布区可以为管理决策提供重要依据, 也有着非常重要的价值[3-4]。但是，由于自然条件、人力、物力的限制，某一区域的物种分布区域数据往往难以获得，或者只能有少数分布点的数据。在这种情况下，利用仅有的信息推断物种的分布区域就显得尤其重要。
川金丝猴(Rhinopithecusroxellanae)是我国特有濒危珍惜野生动物。龙溪-虹口国家级自然保护区内有一定种群的川金丝猴分布，是岷山-邛崃山区川金丝猴分布最接近平原地区的一个种群，但分布点数据缺乏，不足以构建需要分布点数据的栖息地评价模型。2008年地震及随后的次生地质灾害对保护区的生态环境造成了巨大破坏，此时评估区内川金丝猴栖息地受损状况对保护工作开展有极为重要的价值。为此，本文基于已有川金丝猴的栖息地选择文献，确定川金丝猴栖息地生态因子环境包线，结合保护区内各项栖息地生态因子评估龙溪-虹口自然保护区内川金丝猴栖息地分布。在此基础上，根据滑坡体及泥石流分布，分析川金丝猴栖息地受灾害破坏情况。
1 研究地区与方法
1.1 研究区域
四川龙溪-虹口国家级自然保护区位于四川省都江堰市北部，总面积310 km2。保护区与其邻近地区处于中亚热带，北亚热带和青藏高原区三个气候区的多种地理要素交汇的区域，是中国中部的一个大尺度、复合性的生态过渡带。植被垂直分布上，除基带植被亚热带常绿阔叶林受人为活动影响，只在局部地段残遗外，常绿落叶阔叶混交林带、落叶阔叶针叶混交林带、亚高山暗针叶林带、亚高山灌丛带、高山草甸带和高山流石滩稀疏植被均以原生性或原始状态在保护区内连续分布，呈现较完整的植被垂直带谱(图1)。
1.2 研究方法
采用直线包线法(Rectilinear envelope method)评估龙溪-虹口自然保护区内川金丝猴栖息地现状，分析原理是，首先确定研究区域中各项生态因子的适宜范围，通过GIS系统划定研究区域内单个生态因子的适宜区分布，然后叠加多个生态因子，确定单个栅格是否处于各项生态因子适宜范围内。由于龙溪-虹口保护区内的川金丝猴种群栖息地选择信息缺乏，因此我们分析文献资料，基于其他区域的川金丝猴栖息地选择研究和栖息地描述信息建立各项栖息地生态因子的选择范围。在此基础上建立栖息地质量评估模型，并结合保护区现有的少量川金丝猴分布记录，通过ROC判断模型诊断效果并曲线筛选模型评估参数。保护区植被数据来自全国大熊猫第三次调查(1999—2003年)的数据，震后泥石流及滑坡体分布数据也来自全国大熊猫第四次调查(2010—2013年)的数据。保护区内川金丝猴分布记录来自2005年以来保护区开展的各次监测、调查的数据积累。
1.2.1生态因子适宜范围确定
选取海拔、植被、水源距离、坡度、坡位、坡向、人类干扰7项生态因子进行栖息地质量评估。在综合分析四川唐家河自然保护区[5]、白河自然保护区[6]、甘肃白水江自然保护区[7]、四川卧龙保护区[7]、陕西周至保护区[7]和湖北神龙架保护区[7-8]等地的川金丝猴栖息地选择研究的基础上，确定了川金丝猴适宜栖息地各项生态因子的包线范围：
1)适宜海拔： 1 500～3 000 m之间。
2)适宜植被类型：常绿落叶混交林、落叶阔叶林、针阔混交林、针叶林。根据龙溪-虹口保护区内川金丝猴种群观察，川金丝猴虽然也会利用一些中海拔地带的灌丛斑块生境，但是仅仅是偶尔取食活动，其主要的活动还是在高大乔木上。大面积的灌丛并不是川金丝猴的适宜栖息地，因此未纳入适宜植被型范围。
3)水源距离：设置川金丝猴距离水源的适宜范围为1 000 m内，介于郑维超等[5]和李阳[8]的研究结论之间。
4)坡 度：根据唐家河的研究，川金丝猴冬季栖息地偏好坡度为6°～20°，不选择平地[5]。而四川白河自然保护区的研究则表明川金丝猴秋季偏好25°以上陡坡，而不好缓坡[6]。李阳[8]的分析也表明，神龙架的川金丝猴主要偏好坡度15°甚至45°以上的斜坡和陡坡，而不好平缓地带。由于保护区内坡度基本都在35°以下，本文设定坡度适宜范围为大于5°的坡地。
5)坡向：依据郑维超等[5]、顾志宏[6]、李阳[8]的研究，本文设定川金丝猴适宜栖息地的坡向特征为非阴坡。
6)坡位：依据郑维超等[5]、顾志宏[6]、李阳[8]的研究，本文设定川金丝猴的对坡位的偏好为脊部、中部、和下部。
7)人类干扰因子：一般认为，在公路对山区生态环境的直接影响不会超过1 km[9]。龙溪-虹口自然保护区内无居民及各项生产经营活动。保护区外南面有通乡公路一条。保护区外西南侧原有龙池景区，由于2008年地震及其后的山洪、泥石流破坏，景区及公路均已停止运营，居民完全撤出，各种生产活动停止，人为干扰可以忽略。保护区外东南侧有虹口景区，目前正常运营，各项旅游设施和活动区域基本沿公路分布，故与公路影响合并统计。
1.2.2模型评估与参数确定
当物种的家域面积大于评估的单个栅格面积时，物种事实上会同时利用多个栅格的栖息地。这意味着，栖息地适宜度不仅仅由本区域生态因子决定，还受周边栖息地生态因子的影响。因此，我们将目标区域周边一定范围内的包线法识别的适宜栖息地比例作为该区域的适生指数(habitat suitability index)，用于指示该区域的栖息地质量，这事实上是一种邻域平均(Neighbor Averaging)的处理方法。根据陕西周至国家级保护区的研究川金丝猴的全年家域面积10～20 km2左右[10-11]，而起一天之内的移动距离在0.3～3 km之间[10-12]。本研究中我们分别以0.1、0.5、1.0、1.5、2.0、2.5 km和3.0 km为半径(面积为0.03～28.26 km2)进行领域平均估算栖息地的适生指数。然后将自龙溪-虹口自然保护区野外调查及监测过程中收集的川金丝猴分布点作为验证数据，通过受试者工作特征曲线(receiver operating characteristic curve,ROC曲线)计算曲线下面积(Area Under Curve,AUC)值，确定最佳半径，在此基础上划定保护区内的川金丝猴适宜栖息地。一般认为，当0.5
2 结 果
通过包线法直接识别保护区内有川金丝猴适宜栖息地面积69.9 km2，占保护区总面积的22.52%。保护区内川金丝猴的分布记录基本与直接识别的栖息地重合或者非常临近(图2A)。保护区内有川金丝猴分布记录26处(图2)，以此为验证数据，ROC曲线分析表明，基于0.1 ～ 3 km半径估算的适生指数之中，半径2 km的评估得到的适生指数诊断效力最高，AUC为0.7330， 敏感度为0.95，特异性为0.58(图3)。
ROC曲线分析确定的最佳截断点(cutoff point)适生指数为48。基于此，我们对保护区内川金丝猴适宜栖息地进行了判定(图2C)，结果显示保护区有125.8 km2的适宜栖息地，占保护区总面积的310 km2的40.6%。仅有1个已知分布点落在判定的适宜栖息地以外。
5·12地震及其后的滑坡、泥石流对保护区生态环境造成了较大的影响。根据卫片解译的受滑坡体分布，保护区内川金丝猴被破坏的适生指数48以上栖息地面积达45.5 km2，占总适宜栖息地面积的36.1%。将保护区内滑坡体遗迹作为不宜栖息地后，再次计算保护区内栖息地的适生指数，发现地震及其次生灾害对保护区川金丝猴栖息地造成了很大破坏，适生指数均值由震前的41.3±24.8(SD) 降为震后14.8±11.5(SD)。震后保护区内最大适生指数仅为47，全面低于震前适宜栖息地的划分标准(图2D、图4)。
3 讨 论
基于文献资料的包线法评估栖息地质量的诊断效力并非很高，但是这一方法的意义在于，当对研究区域的物种分布资料缺乏，但是被研究物种又有着较为丰富的栖息地选择研究资料时，我们可以依靠现有资料对物种可能分布区作出具有一定可靠性的快速评估判断。保护区内川金丝猴记录点仅有26处，对于ENFA等需要分布点数据的评估模型而言，这远远低于训练模型所需的数据样本量，跟不用说还需预留样本作为检验。
龙溪-虹口自然保护区历年监测调查纪录到的川金丝猴分布点基本上都落在包线法确定的适宜栖息地内部或者周边非常邻近的区域。通过2 km半径邻域平均得到的适生指数的AUC值达到0.7330，具有一定的实用价值。这一半径也处于各地区确定的川金丝猴研究家域大小范围之内。事实上，由于我们只有分布点(present)数据，而没有非分布点(absence)数据，ROC曲线分析中只能使用随机抽取的伪非分布点(false absence)进行估算，因而造成ROC曲线右移, 相应地AUC值会偏小。所以包线法和领域平均的真实诊断效力应该还会更高一点[4,14]。
基于最佳截断点判定的川金丝猴适宜栖息地主要分为两大块。除了保护区中心及东部区域的主要适宜栖息地外，在保护区西南角狭长地带还有几块面积相对较小的适宜栖息地。本分析仅仅针对保护区内的区域，如果将保护区外的区域也纳入分析，西南角的几块栖息地极有可能是完整的一大块。这也提示我们，如果将保护区西南角狭长地带的边界扩大，将更有利于保护该川金丝猴种群。两大块栖息地之间的距离很近，仅仅1 km左右，此外，对照包线法直接确定的适宜栖息地，这两大块适宜栖息地之间事实上存在着大量的适宜栖息地斑块碎片。因此这两块栖息地之间的种群交流应该不存在困难。
地震及其次生地质灾害造成的滑坡体在保护区内分布广泛，导致川金丝猴栖息地质量全面下降，保护区内所有区域的适生指数全部低于震前的适宜栖息地标准。其中下降最为严重的是保护区西南角狭长地带和保护区东南部。但是这并不意味着保护区内川金丝猴栖息地已经丧失，低于震前标准的栖息地仍然可以支持保护区内的川金丝猴种群。在2012年保护区的调查、监测中也记录到不少的川金丝猴活动痕迹。
当前预测物种分布区的方法很多，可以分为两大类：一类是根据已知物种分布数据预测物种分布区，而另一类则是根据物种对各种生态因子的偏好选择预测物种分布区。
第一类是是根据已知物种分布数据建立栖息地适宜度(habitat suitability)评估模型，推断物种整个区域的分布状况。早期的预测方法基于逻辑斯蒂的广义线性模型(generalized linear models,GLMs)或者广义加性模型(and generalized additive models,GAMs)[15]。这些方法不仅仅要求有物种分布点的生态因子数据，还需要作为对照的非分布点生态因子数据。由于大量物种分布资料仅仅记载了分布点的数据，而没有采集对照点的数据，因此这一类方法的应用受到了很大的限制。为回避这一限制，基于生态位假说，目前已经发展出一系列不需要对照点、只要有分布点(present-only)的栖息地适宜度评估模型，比如BIOCLIM[16]，DOMAIN[17]，ENFA[18]，Maxent[19]等；这一系列方法。但是，当我们需要对一个完全缺乏物种分布信息的的区域进行栖息地适宜度评估时，这一类方法就难以实施了。虽然可以将一个区域建立的栖息地适宜度评估模型移用到其他区域，但是由于地理特征的变化，移用的模型的准确性往往会受到较大影响[20]。
本文采用的是第二类方法，这类方法不需要研究范围内物种分布点数据,而是基于物种的对环境生态因子的要求,通过环境包线(environmental envelope)确定物种分布范围[21-22]。CLIMEX模型即是一种基于已有的物种生理生态参数知识,根据气候条件推断物种潜在分布区的方法,但是这种方法也可以利用物种原有分布区气候资料推断物种分布气候包线,并非完全依据既有知识[23]。在国内，野生动物，尤其是濒危珍惜野生动物的栖息地选择研究则开展的非常广泛[24]。通过文献资料往往可以获得研究物种在各个分布区域的定量或者定性描述的栖息地选择研究资料。如果充分利用这些资料，构建环境包线，则可以快速推断缺少研究区域的物种分布和潜在分布信息。
致谢：本文研究过程中，中科院成都生物所戴强副研究员；杨彪博士、龙溪-虹口国家级自然保护区全体员工提供了大力支持和帮助。在此一并致谢！
参考文献：
[1]GUISAN A,ZIMMERMANN N E.Predictive habitat distribution models in ecology[J].Ecological Modelling,2000,135(2/3):147-186.
[2]SANTOS X,BRITO J C,CARO J,et al.Habitat suitability,threats and conservation of isolated populations of the smooth snake (Coronella austriaca) in the southern Iberian Peninsula[J].Biological Conservation,2009,142(2):344-352.
[3]BIDINGER K,LOTTERS S,RODDER D,et al.Species distribution models for the alien invasive Asian Harlequin ladybird (Harmonia axyridis)[J].Journal of Applied Entomology,2012,136(1/2):109-123.
[4]王运生,谢丙炎,万方浩，等.ROC曲线分析在评价入侵物种分布模型中的应用[J].生物多样性,2007,15(4):365-372.
[5]郑维超,黎大勇,谌利民，等.唐家河国家级自然保护区川金丝猴冬季栖息地选择[J].四川动物,2012,31(2):208-211.
[6]顾志宏,金崑,刘世荣，等.四川省白河自然保护区川金丝猴生境评价[J].林业科学,2007,43(8):96-99.
[7]刘广超.川金丝猴栖息地质量评价和保护对策研究[D].北京：北京林业大学,2007.
[8]李阳.“3S”技术在神农架川金丝猴栖息地选择及动态迁徙中的应用[D].北京：北京林业大学,2010.
[9]FORMAN R T T and ALEXANDER L E.Roads and their major ecological effect [J].Annual Review of Ecology & Systematics,1998,29(4):207-231.
[10]LI B,CHEN C,JI W,et al.Seasonal Home Range Changes of the Sichuan Snub-Nosed Monkey(Rhinopithecus roxellana) in the Qinling Mountains of China[J].Folia primatologica;international journal of primatology,2000,71(6):375.
[11]TAN C L,GUO S,Li B.Population Structure and Ranging Patterns of Rhinopithecus roxellana,in Zhouzhi National Nature Reserve,Shaanxi,China[J].International Journal of Primatology,2007,28(3):577-591.
[12]李义明,廖明尧,喻杰等.社群大小的年变化、气候和人类活动对神农架自然保护区川金丝猴日移动距离的影响[J].生物多样性,2005,13(5):432-438.
[13]ARAUJO M B,PEARSON R G,THUILLER W,et al.Validation of species-climate impact models under climate change.Global Change Biology[J].2005,11(9):1504-1513.
[14]PHILLIPS S J,ANDERSON R P,SCHAPIRE R E.Maximum entropy modeling of species geographic distributions[J].Ecological Modelling,2006,190(3/4):231-259.
[15]GUISAN A,THOMAS C EDWARDS Jr,HASTIE T.Generalized linear and generalized additive models in studies of species distributions: setting the scene[J].Ecological Modelling,2002,157(2/3):89-100.
[16]BUSBY J R.A biogeoclimatic analysis of Nothofagus cunninghamii (Hook.) Oerst.in southeastern Australia[J].Austral Ecology,1986,11(1):1-7.
[17]CARPENTER G,GILLISON A N,WINTER J.DOMAIN:a flexible modelling procedure for mapping potential distributions of plants and animals[J].Biodiversity & Conservation,1993,2(6):667-680.
[18]HIRZEL A H,HAUSSER J,CHESSEL D,et al.Ecological-Niche Factor Analysis:How to Compute Habitat-Suitability Maps without Absence Data?[J].Ecology,2002,83(7):2027-2036.
[19]PHILLIPS S J,SCHAPIRE R E.A maximum entropy approach to species distribution modeling[M]//Proceedings of the 21st International Conference on Machine Learning.ACM Press,2004：655-662.
[20]RANDIN C F,THOMAS D,STEFAN D,et al.Are niche-based species distribution models transferable in space?[J].Journal of Biogeography,2010,33(10):1689-1703.
[21]FISCHER J,LINDENMAYER D B,NIX H A,et al.Climate and animal distribution:a climatic analysis of the Australian marsupial Trichosurus caninus[J].Journal of Biogeography.2001,28(3):293-304.
[22]邓浩,纪力强.物种潜在分布地预测系统PSDS的设计与实现[J].生物多样性,2008,16(1):96-102.
[23]SUTHERST R W,MAYWALD G F,RUSSELL B L.Estimating vulnerability under global change:modular modelling of pests[J].Agriculture Ecosystems & Environment,2000,82(1/3):303-319.
[24]蒋志刚,李春旺,曾岩.中国野生动物研究进展[C]//第五届全国生物多样性保护与持续利用研讨会.杭州，2002.
