气候环境变化对冬虫夏草产量的影响与虫草产量预测模型研究

朱寿东　黄璐琦　郭兰萍　马兴田　郝庆秀　乐智勇　张小波　杨光　张燕　陈美兰

[摘要] 冬虫夏草是我国特有的名贵中药材，其道地产区是青藏高原那曲、昌都，青海玉树等地区。最近10年，虫草资源在盲目和过度的采挖中已经到达临界点，如何加快保护、发展和合理利用冬虫夏草这一宝贵资源，已提到了冬虫夏草研究和利用的重要议事日程。青藏高原生态环境和气候变化比其他地区具有超前变化趋势，这意味着冬虫夏草种群的分布和演化与其他种群比起来更明显和剧烈，该文基于RS（remote sensing，遥感）/GIS（geographic information system，地理信息系统）技术，利用雪线高程，采收期（4，5月）的平均气温，降水量和日照时长等环境因子和虫草产量之间的关系，建立加权几何平均模型，在每年的6月上旬预测当年的虫草产量，经过检验，该模型的精度可达82.16%以上。该研究可为虫草产业的健康、可持续发展提供基础数据和信息。

[关键词] 冬虫夏草； 遥感； 环境因子； 加权几何平均； 产量预测

Climate change impacts on yield of Cordyceps sinensis and

research on yield prediction model of C. sinensis

ZHU Shoudong1， HUANG Luqi1*， GUO Lanping1， MA Xingtian2， HAO Qingxiu1，

LE Zhiyong2， ZHANG Xiaobo1， YANG Guang1， ZHANG Yan1， CHEN Meilan1

（1.State Key Laboratory Breeding Base of Daodi Herbs， National Resource Center for Chinese

Materia Medica， China Academy of Chinese Medical Sciences， Beijng 100700， China；

2.Kangmei Pharmaceutical Co.， Ltd.， Shenzhen 518000， China）

[Abstract] Cordyceps sinensis is a Chinese unique precious herbal material， its genuine producing areas covering Naqu， Changdu in Qinghai Tibet Plateau， Yushu in Qinghai province and other regions. In recent 10 years， C. sinensis resources is decreasing as a result of the blindly and excessively perennial dug. How to rationally protect， develop and utilize of the valuable resources of C. sinensis has been referred to an important field of research on C. sinensis. The ecological environment and climate change trend of Qinghai Tibet plateau happens prior to other regions， which means that the distribution and evolution of C. sinensis are more obvious and intense than those of the other populations. Based on RS （remote sensing）/GIS（geographic information system） technology， this paper utilized the relationship between the snowline elevation， the average temperature， precipitation and sunshine hours in harvest period （April and may） of C. sinensis and the actual production of C. sinensis to establish a weighted geometric mean model. The model′s prediction accuracy can reach 82.16% at least in forecasting C. sinensis year yield in Naqu area in every early June. This study can provide basic datum and information for supporting the C. sinensis industry healthful， sustainable development.

[Key words] Cordyceps sinensis； remote sensing； environmental factor； weighted geometric mean； yield prediction

目前世界上已知的廣义虫草属有500多种，其中冬虫夏草是我国特有的名贵中药资源，是由冬虫夏草菌Cordyceps sinensis（Berk.） Sacc.侵染鳞翅目蝙蝠蛾幼虫后发育而成的由子座和菌核组成的菌体。冬草夏草一般生长在海拔3 000 m以上的山地阴坡、半阴坡的灌丛和草甸中，在我国主要分布于青藏高原及边缘地区[1]，零星分布于境外的尼泊尔、不丹的喜马拉雅山脉高寒草甸，我国冬虫夏草的采挖量占世界总采挖量的98%以上[2]。

由于天然冬虫夏草的生长有严格的寄生性及要求特殊的生态地理环境，野生资源十分有限[3]。另一方面，目前冬虫夏草不能进行规模化人工培育[3]，只能自然生长，于是其蕴藏量在日益趋高的价格逼迫下迅速下降。据统计，20世纪50年代以前全国产量曾达到1 000 t以上，60年代初为5 080 t，自70年代以来，因全球气候变化、生态环境破坏，以及掠夺式的采挖，加之其自身生长缓慢，自然资源更新能力低下，更使天然蕴藏量迅速减少[4]。据80年代中期第三次全国中草药资源普查结果显示其野生蕴藏量仅约400 t，90年代就仅有300 t[5]，总体还呈逐年下降趋势。目前冬虫夏草天然资源已濒于灭绝，被列为国家二级重点保护野生植物。

由于各种原因从第三次中药资源普查至今，没有对虫草进行过全面的调查，特别是虫草在西藏的具体分布和蕴藏量，有关部门只能从各乡群众采挖范围上获取一些基本情况，基本上是经验数据。因此，准确掌握虫草的产量信息是可持续开发冬虫夏草资源的基础，也是一个亟待解决的问题。

杨大荣等[6]通过研究表明，冬虫夏草的分布受到地形、地貌、海拔、气候等众多因子的影响。吴庆贵等[7]通过样方法指出植被盖度、植物数量、坡度、天敌因子、温度和湿度是影响冬虫夏草种群分布的主导因子。李芬等[8]通过文献调研和野外调查相结合的方法选择海拔高度、植被和土壤作为影响冬虫夏草生长的环境因子。李晖等[9]指出西藏地区冬虫夏草分布于年均温在[-1.9 ℃， 4.7 ℃]，年降水量在[370 mm， 850 mm]。综合分析发现，冬虫夏草的分布和气候环境变化关系较大，因此本文选择雪线高程，采收期（4月和5月）的平均气温，降水量和日照时长等环境因子，分别建立这些因子与虫草产量的线性关系，利用独立性权数法确定系数（相关系数比值作为权重），建立加权几何平均模型，来预测当年的虫草产量。

1 材料

1.1 实地调查与采样

本文在2013—2014年通过实地考察、问卷调查等方法对众多环境因子进行筛选，结合冬虫夏草的生长周期，根据经验初步选出虫草所在地区的5月份雪线高程，4，5，6月的降水量，4，5，6月的平均气温，4，5，6月的日照时长共10个因子作为影响虫草产量的因子，供后期建立预测模型使用。

1.2 降水量、平均气温、日照时长等相关因子

根据虫草所在地区的经纬度，从中国气象科学数据共享服务网（http：//cdc.nmic.cn）找到该区域内的气象台站，利用这些气象台站的平均值作为环境因子的值。按照这样的方法得到4，5，6月的降水量，平均气温，日照时长，以及由此衍生的4月5月的降水量和，平均气温和，日照时长和。

1.3 研究区介绍

那曲地区地处西藏北部，位于青藏高原腹地，是长江、怒江、拉萨河、易贡河等大江大河的源头。与新疆维吾尔自治区和青海省交界，东邻昌都市，南接拉萨市、日喀则市、林芝市，西与阿里地区相连。境域面积45万km2，地理位置约E83°55′—E95°5′， N29°55′—N36°30′。那曲地区除东部少量的半农半牧区外，基本上是一个纯牧区，拥有草地面积5亿多亩。那曲地区是冬虫夏草的道地产区，并且是虫草分布核心带。

2 方法

2.1 构建雪线高程

本文利用NDSI指数（normalized difference snow index， 归一化雪被指数）提取雪线[10]。NDSI是观测冰雪定量指标，积雪有很强的可见光反射和很强的短波红外吸收特性。积雪在0.5 μm附近有高反射率，在1.6，2.1 μm附近反射率较低。

NDSI=（Ref0.555 μm-Ref1.640 μm）/（Ref0.555 μm+Ref1.640 μm） （1）

公式（1）中Ref0.555 μm，Ref1.640 μm分别为0.555，1.64 μm处的反射率，NDSI≥0.4的区域即为冰雪覆盖区。

2.1.1 把DN值图像转化成反射率图像 根据遥感影像的来源，有2种情况，一种是针对Landsat TM5/Landsat TM7影像，另一种是针对Landsat TM8影像。

对于Landsat TM5和Landsat TM7影像：0.555，1.64 μm波段分别对应着影像的第2和第5波段。

第一步，分别计算各个波段每个像元的辐射亮度L值。

L=Lmax-LminQCALmax-QCALmin×（QCAL-QCALmax）+Lmin （2）

公式（2）中，QCAL为某一像元的DN值（digital number， 像元亮度值），即QCAL=DN。QCALmax为像元可以取的最大值，QCALmin为像元可以取的最小值。Lmax为该波段的最大辐射亮度，Lmin为该波段的最小辐射亮度，可以从遥感影像的头文件中获取以上参数的信息。

第二步，计算各波段反射率（反照率、反射率）ρ。

ρi=π×L×D2ESUN×cosθ（i为第i波段） （3）

公式（3）中，D为日地距離，可从影像的头文件中获取该参数信息。ESUN（solar exoatmospheric spectral irradiances，或ESUNI）为大气顶层太阳辐照度可从遥感权威单位定期测定并公布的信息中获取，见表1。

2.1.2 在生成的反射率图像上计算NDSI 对于Landsat TM5/7，NDSI=（B2-B5）/（B2+B5）；对于Landsat TM8，NDSI=（B3-B6）/（B3+B6）。式中B2， B5， B3， B6分别代表各波段反射率图像。

2.1.3 提取NDSI上大于等于0.4的区域 根据MODIS产品说明，认为NDSI大于等于0.4的区域被积雪覆盖。利用ENVI软件里的band math工具来提取积雪覆盖的区域。

2.1.4 提取积雪区域的边界 由于提取的积雪覆盖的区域是栅格图形，利用ENVI软件转化为矢量后是多边形，接着利用Arctools工具把多边形转化为线类型的矢量图形，至此雪线已提取完成。

2.1.5 利用空间分析提取雪线处最低高程 利用Arcmap 10软件，加载DEM（digital elevation model，数字高程模型）数据，与雪线分布数据进行空间分析，提取雪线处的高程信息，最后去除图像边缘的异常值。

2.2 气候环境因子与虫草产量的关系

利用相关分析工具，分别确定各因子与虫草产量之间的相关系数，按照相关系数大小排序来确定入选模型的因子，各因子与虫草产量的相关系数，就是该因子的权重。确定了因子后，基于时间序列的虫草产量数据和各因子数据建立虫草产量的单因子线性回归模型。

2.3 虫草产量预测模型的建立

目前主流的赋权方法主要分为2类，一类是主观赋权法，主要是以德尔菲法（Delphi），层次分析法（analytic hierarchy process， AHP），专家会议法为代表。而另一类客观赋权法则是以主成分分析法，因子分析法，熵值法，离差、均差法等纯数学分析为代表[11]。本文通过对比试验，采用客观赋权法中的独立性权系数法求出各因子的权重系数，然后采用加权几何平均模型来预测虫草产产量。

2.3.1 独立性权系数法原理 独立性权系数法[3]是根据各指标与其他指标之间的共线性强弱来确定指标权重的。

设有指标项，若指标与其他指标的复相关系数越大，则说明与其他指标之间的共线性关系越强，越容易由其他指标的线性组合表示，重复信息越多，因此该指标的权重也就应该越小。

即若指标与其他指标的复相关系数R越大，该指标的权重越小。

R=（y-）（-）（y-）2（-）2

由于本文中强调各相关因子与虫草产量的关系，因此用各因子与虫草产量的相关系数替代了各因子间的复相关系数。

2.3.2 加权几何平均模型 如果预测因子与虫草产量的线性函数为fi，i=1， 2， …， N，N为最终入选模型的因子数目。各个因子权重比值为wi，则虫草产量Y为公式（6）。

Y=（∏Nf wii）1Ni=1wi（6）

2.4 虫草产量数据调研和采集

在2013，2014年6月，通过实地走访那曲地区聂荣县、比如县、索县和那曲虫草交易市场，和采挖虫草的藏民进行交流，获取了虫草的采挖情况。结合查阅资料获取了近6年的那曲虫草产量[1521]，见表2。

2.5 生态因子与虫草产量的关系

2.5.1 雪线高程 根据虫草的采挖时间，选择收集那曲地区的5月份的landsat卫星影像数据，如果5月份数据含云量特别大，可选择4月份的数据。那曲地区所在的影像条带号是：WRS_PATH=137，WRS_ROW=038。最终选择的数据情况见表3。

按照2.1的方法对遥感进行处理，得到雪线处的高程信息见表4，提取雪線的过程中虽然删除了影像边缘的异常值，但是影像零散的冰雪覆盖区对雪线的边界有影响，根据正态分布原理[12]，把提取出来的雪线高程的“平均值-2×标准差”作为雪线的最低高程。按照2012年、2014年、2011年、2009年、2010年的顺序，雪线的最低高程依次上升，见图1。

由于2013年的影像含云量较大，建模过程中舍去这一年的数据。利用剩余的5期数据，计算得雪线高程和虫草产量的相关系数r=0.40，线性回归关函数为y=0.030 1x-99.367。

2.5.2 4月和5月降水量之和 选择那曲地区附近的5个气象台站（编号为55294，55299，55493，56106，56202）的降水量数据的平均值作为该因子的值。通过比较单月降水量、各月降水量之和与虫草

产量的相关系数发现，4月和5月降水量之和比单独月降水量与虫草产量的相关性更大，于是选择4月和5月降水量之和作为预测因子。那曲地区各年4月和5月降水量之和见表5。

计算得到4月和5月降水量的和与虫草产量的相关系数r=0.95，线性回归关函数为y=0.13x-49.44。

2.5.3 5月平均气温 同降水量一样，选择5个气象台站的数据平均值作为该因子的值。通过比较单月平均气温、各月平均气温之和与虫草产量的相关系数发现，5月份的平均气温比其他月、月份之和与虫草产量的相关性更大，于是选择5月平均气温作为预测因子，那曲地区各年5月平均气温见表5。计算得到5月平均气温与虫草产量的相关系数r=-0.69，线性回归关函数为y=-1.477 1x+114.09。

2.5.4 4月和5月日照时长之和 通过比较单月日照时长、各月日照时长之和与虫草产量的相关系数发现，4月和5月日照时长之和比其他月、月份的和与虫草产量的相关性更大，于是选择4月和5月日照时长之和作为预测因子，那曲地区各年4月和5月日照时长之和见表5。计算得到5月平均气温与虫草产量的相关系数r=-0.50，线性回归关函数为y=-0.022 1x+144.4。

2.6 构建虫草产量预测模型

根据各个因子与虫草产量的相关系数，得到相关系数的比值是r1∶r2∶r3∶r4=0.4∶0.95∶|-0.69|∶|-0.50|=1∶3∶2∶1.于是各因子权重比值w1∶w

上式中，X1是雪线高程，X2是4月和5月降水量的和，X3是5月平均气温，X4是4月和5月日照时长的和，Y是预测的虫草产量（单位t）。

3 结果

3.1 模型精度验证

把2009—2014年的环境因子数据带入得到的模型中，得到的预测值和实际值见表6。特别地，2013年的遥感影像含云量较大，无法提取雪线高程，于是对预测模型进行修正，去掉f1，只利用剩余的3个因子进行预测，得到2013年的虫草产量预测值。综合过去6年的预测结果和真实值，计算得到模型预测的平均偏差是7.98%，最大偏差是17.84%。

3.2 2015年虫草产量预测

在2015年6月上旬，根据生态因子的构建方法，得到X1雪线高程是4 574 m，X2 4月和5月降水量之和是80.72 mm，X3 5月平均气温是5.08 ℃，X4 4月和5月日照时长之和由于数据暂缺无法获取，所以去掉f4，只利用前3个因子进行预测，得到2015年那曲地区的虫草产量是46.42 t，比去年增产约60%，是2013年产量的85%左右。综合分析来看，今年高产与今年虫草生长季降雨量丰沛有较大关系。

4 讨论

本文通过RS/GIS技术来构建气候环境因子，把气候环境因子与虫草产量进行线性拟合，配合几何加权平均模型，得到了虫草产量的预测模型。经过检验，发现模型预测偏差较小。而且，如果有非重要的生态因子因为实际情况（天气或者时间关系）暂缺时，模型的精度不会受到较大影响，具有较高的稳定性和灵活性。

[参考文献]

[1] 章力建，李兵，胡育骄，等.中国冬虫夏草资源管理概况[J]. 中国草地学报， 2010， 32（增刊）： 1.

[2] 刘天平，卓嘎，郭健斌. 西藏冬虫夏草资源可持续利用探讨[J]. 西藏研究， 2010（3）： 114.

[3] 蔡佩云，孙世英. 青海省冬虫夏草资源持续开发利用浅析[J]. 中国草地学报， 2010， 32（增刊）： 6.

[4] 张古忍，吴光国，古德祥，等. 冬虫夏草资源可持续利用面临的问题及其对策[J]. 中国草地学报[J]， 2010， 32（增刊）： 51.

[5] 陈仕江，钟国跃，马开森. 冬虫夏草资源可持续利用的思考[J]. 中国草地学报， 2010， 32（增刊）： 44.

[6] 杨大荣，彭艳琼，陈吉岳，等. 中国冬虫夏草分布格局与环境变化对其分布的影响[J]. 中国草地学报， 2010，32（增刊）： 22.

[7] 吴庆贵，苏智先，苏瑞军，等. 冬虫夏草生境选择主导因子[J]. 广西植物， 2009，29（3）： 331.

[8] 李芬，吴志丰，徐翠，等. 三江源区冬虫夏草资源适宜性空间分布[J]. 生态学报， 2014，34（5）： 1318.

[9] 李晖，王立辉. 西藏冬虫夏草生境的气候[J]. 西藏科技， 2011（6）： 68.

[10] 赵英时. 遥感应用分析原理与方法[M]. 北京： 科学出版社， 2003.

[11] 李亮. 評价中权系数理论与方法比较[M]. 上海： 上海交通大学， 2009.

[12] 刘自远，刘成福. 综合评价中指标权重系数确定方法探讨[J]. 中国卫生质量管理， 2006， 13（2）： 44.

[13] 蔡辉，沈婉兰. 医院综合效益评价中的权数（二）——独立性权数[J]. 中国医院统计， 1997， 4（3）：131.

[14] 杜子芳. 抽样技术及其应用[M]. 北京：清华大学出版社， 2005.

[15] 李芬，吴志丰，徐翠，等. 三江源区冬虫夏草资源适宜性空间分布[J]. 生态学报， 2014，34（5）： 1318.

[16] 新华网. 西藏冬虫夏草价格进入下行通道[EB/OL]. （20100705） [20160209].http：//news.xinhuanet.com/fortune/201007/05/c_12300368.htm .

[17] 西藏商城藏善堂. 新产冬虫夏草价格2010年5月走势分析[EB/OL]. （20100423） [20160209 ].http：//www.xizangshop.com/article138.html.

[18] 冬虫夏草屋. 同仁堂虫草价格调价涨幅达三成[EB/OL]. （20120602） [20160209].http：//www.dcxc5.com/article488.html.

[19] 冬虫夏草屋. 虫草产量连年下降，冬虫夏草资源日趋稀缺[EB/OL]. （20120602）[20160209].http：//www.dcxc5.com/article585.html.

[20] 中国经济网. 青海计划今年采集冬虫夏草35吨[EB/OL]. （20140529） [20160209].http：//city.ce.cn/ztpd/ysgjccj/ysxwlb/yshdxc/201405/29/t20140529_1653483.shtml.

[21] 中国经济网. 2014年新草动态[EB/OL]. （20140529） [20160209].http：//city.ce.cn/ztpd/ysgjccj/ysxwlb/yshdxc/201405/29/t20140529_1653566.shtml.

[责任编辑 吕冬梅]