氣候環(huán)境變化對冬蟲夏草產(chǎn)量的影響與蟲草產(chǎn)量預測模型研究

朱壽東　黃璐琦　郭蘭萍　馬興田　郝慶秀　樂智勇　張小波　楊光　張燕　陳美蘭

[摘要] 冬蟲夏草是我國特有的名貴中藥材，其道地產(chǎn)區(qū)是青藏高原那曲、昌都，青海玉樹等地區(qū)。最近10年，蟲草資源在盲目和過度的采挖中已經(jīng)到達臨界點，如何加快保護、發(fā)展和合理利用冬蟲夏草這一寶貴資源，已提到了冬蟲夏草研究和利用的重要議事日程。青藏高原生態(tài)環(huán)境和氣候變化比其他地區(qū)具有超前變化趨勢，這意味著冬蟲夏草種群的分布和演化與其他種群比起來更明顯和劇烈，該文基于RS（remote sensing，遙感）/GIS（geographic information system，地理信息系統(tǒng)）技術(shù)，利用雪線高程，采收期（4，5月）的平均氣溫，降水量和日照時長等環(huán)境因子和蟲草產(chǎn)量之間的關(guān)系，建立加權(quán)幾何平均模型，在每年的6月上旬預測當年的蟲草產(chǎn)量，經(jīng)過檢驗，該模型的精度可達82.16%以上。該研究可為蟲草產(chǎn)業(yè)的健康、可持續(xù)發(fā)展提供基礎(chǔ)數(shù)據(jù)和信息。

[關(guān)鍵詞] 冬蟲夏草； 遙感； 環(huán)境因子； 加權(quán)幾何平均； 產(chǎn)量預測

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.侵染鱗翅目蝙蝠蛾幼蟲后發(fā)育而成的由子座和菌核組成的菌體。冬草夏草一般生長在海拔3 000 m以上的山地陰坡、半陰坡的灌叢和草甸中，在我國主要分布于青藏高原及邊緣地區(qū)[1]，零星分布于境外的尼泊爾、不丹的喜馬拉雅山脈高寒草甸，我國冬蟲夏草的采挖量占世界總采挖量的98%以上[2]。

由于天然冬蟲夏草的生長有嚴格的寄生性及要求特殊的生態(tài)地理環(huán)境，野生資源十分有限[3]。另一方面，目前冬蟲夏草不能進行規(guī)?；斯づ嘤齕3]，只能自然生長，于是其蘊藏量在日益趨高的價格逼迫下迅速下降。據(jù)統(tǒng)計，20世紀50年代以前全國產(chǎn)量曾達到1 000 t以上，60年代初為5 080 t，自70年代以來，因全球氣候變化、生態(tài)環(huán)境破壞，以及掠奪式的采挖，加之其自身生長緩慢，自然資源更新能力低下，更使天然蘊藏量迅速減少[4]。據(jù)80年代中期第三次全國中草藥資源普查結(jié)果顯示其野生蘊藏量僅約400 t，90年代就僅有300 t[5]，總體還呈逐年下降趨勢。目前冬蟲夏草天然資源已瀕于滅絕，被列為國家二級重點保護野生植物。

由于各種原因從第三次中藥資源普查至今，沒有對蟲草進行過全面的調(diào)查，特別是蟲草在西藏的具體分布和蘊藏量，有關(guān)部門只能從各鄉(xiāng)群眾采挖范圍上獲取一些基本情況，基本上是經(jīng)驗數(shù)據(jù)。因此，準確掌握蟲草的產(chǎn)量信息是可持續(xù)開發(fā)冬蟲夏草資源的基礎(chǔ)，也是一個亟待解決的問題。

楊大榮等[6]通過研究表明，冬蟲夏草的分布受到地形、地貌、海拔、氣候等眾多因子的影響。吳慶貴等[7]通過樣方法指出植被蓋度、植物數(shù)量、坡度、天敵因子、溫度和濕度是影響冬蟲夏草種群分布的主導因子。李芬等[8]通過文獻調(diào)研和野外調(diào)查相結(jié)合的方法選擇海拔高度、植被和土壤作為影響冬蟲夏草生長的環(huán)境因子。李暉等[9]指出西藏地區(qū)冬蟲夏草分布于年均溫在[-1.9 ℃， 4.7 ℃]，年降水量在[370 mm， 850 mm]。綜合分析發(fā)現(xiàn)，冬蟲夏草的分布和氣候環(huán)境變化關(guān)系較大，因此本文選擇雪線高程，采收期（4月和5月）的平均氣溫，降水量和日照時長等環(huán)境因子，分別建立這些因子與蟲草產(chǎn)量的線性關(guān)系，利用獨立性權(quán)數(shù)法確定系數(shù)（相關(guān)系數(shù)比值作為權(quán)重），建立加權(quán)幾何平均模型，來預測當年的蟲草產(chǎn)量。

1 材料

1.1 實地調(diào)查與采樣

本文在2013—2014年通過實地考察、問卷調(diào)查等方法對眾多環(huán)境因子進行篩選，結(jié)合冬蟲夏草的生長周期，根據(jù)經(jīng)驗初步選出蟲草所在地區(qū)的5月份雪線高程，4，5，6月的降水量，4，5，6月的平均氣溫，4，5，6月的日照時長共10個因子作為影響蟲草產(chǎn)量的因子，供后期建立預測模型使用。

1.2 降水量、平均氣溫、日照時長等相關(guān)因子

根據(jù)蟲草所在地區(qū)的經(jīng)緯度，從中國氣象科學數(shù)據(jù)共享服務網(wǎng)（http：//cdc.nmic.cn）找到該區(qū)域內(nèi)的氣象臺站，利用這些氣象臺站的平均值作為環(huán)境因子的值。按照這樣的方法得到4，5，6月的降水量，平均氣溫，日照時長，以及由此衍生的4月5月的降水量和，平均氣溫和，日照時長和。

1.3 研究區(qū)介紹

那曲地區(qū)地處西藏北部，位于青藏高原腹地，是長江、怒江、拉薩河、易貢河等大江大河的源頭。與新疆維吾爾自治區(qū)和青海省交界，東鄰昌都市，南接拉薩市、日喀則市、林芝市，西與阿里地區(qū)相連。境域面積45萬km2，地理位置約E83°55′—E95°5′， N29°55′—N36°30′。那曲地區(qū)除東部少量的半農(nóng)半牧區(qū)外，基本上是一個純牧區(qū)，擁有草地面積5億多畝。那曲地區(qū)是冬蟲夏草的道地產(chǎn)區(qū)，并且是蟲草分布核心帶。

2 方法

2.1 構(gòu)建雪線高程

本文利用NDSI指數(shù)（normalized difference snow index， 歸一化雪被指數(shù)）提取雪線[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的區(qū)域即為冰雪覆蓋區(qū)。

2.1.1 把DN值圖像轉(zhuǎn)化成反射率圖像 根據(jù)遙感影像的來源，有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為該波段的最小輻射亮度，可以從遙感影像的頭文件中獲取以上參數(shù)的信息。

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

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

公式（3）中，D為日地距離，可從影像的頭文件中獲取該參數(shù)信息。ESUN（solar exoatmospheric spectral irradiances，或ESUNI）為大氣頂層太陽輻照度可從遙感權(quán)威單位定期測定并公布的信息中獲取，見表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的區(qū)域 根據(jù)MODIS產(chǎn)品說明，認為NDSI大于等于0.4的區(qū)域被積雪覆蓋。利用ENVI軟件里的band math工具來提取積雪覆蓋的區(qū)域。

2.1.4 提取積雪區(qū)域的邊界 由于提取的積雪覆蓋的區(qū)域是柵格圖形，利用ENVI軟件轉(zhuǎn)化為矢量后是多邊形，接著利用Arctools工具把多邊形轉(zhuǎn)化為線類型的矢量圖形，至此雪線已提取完成。

2.1.5 利用空間分析提取雪線處最低高程 利用Arcmap 10軟件，加載DEM（digital elevation model，數(shù)字高程模型）數(shù)據(jù)，與雪線分布數(shù)據(jù)進行空間分析，提取雪線處的高程信息，最后去除圖像邊緣的異常值。

2.2 氣候環(huán)境因子與蟲草產(chǎn)量的關(guān)系

利用相關(guān)分析工具，分別確定各因子與蟲草產(chǎn)量之間的相關(guān)系數(shù)，按照相關(guān)系數(shù)大小排序來確定入選模型的因子，各因子與蟲草產(chǎn)量的相關(guān)系數(shù)，就是該因子的權(quán)重。確定了因子后，基于時間序列的蟲草產(chǎn)量數(shù)據(jù)和各因子數(shù)據(jù)建立蟲草產(chǎn)量的單因子線性回歸模型。

2.3 蟲草產(chǎn)量預測模型的建立

目前主流的賦權(quán)方法主要分為2類，一類是主觀賦權(quán)法，主要是以德爾菲法（Delphi），層次分析法（analytic hierarchy process， AHP），專家會議法為代表。而另一類客觀賦權(quán)法則是以主成分分析法，因子分析法，熵值法，離差、均差法等純數(shù)學分析為代表[11]。本文通過對比試驗，采用客觀賦權(quán)法中的獨立性權(quán)系數(shù)法求出各因子的權(quán)重系數(shù)，然后采用加權(quán)幾何平均模型來預測蟲草產(chǎn)產(chǎn)量。

2.3.1 獨立性權(quán)系數(shù)法原理 獨立性權(quán)系數(shù)法[3]是根據(jù)各指標與其他指標之間的共線性強弱來確定指標權(quán)重的。

設(shè)有指標項，若指標與其他指標的復相關(guān)系數(shù)越大，則說明與其他指標之間的共線性關(guān)系越強，越容易由其他指標的線性組合表示，重復信息越多，因此該指標的權(quán)重也就應該越小。

即若指標與其他指標的復相關(guān)系數(shù)R越大，該指標的權(quán)重越小。

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

由于本文中強調(diào)各相關(guān)因子與蟲草產(chǎn)量的關(guān)系，因此用各因子與蟲草產(chǎn)量的相關(guān)系數(shù)替代了各因子間的復相關(guān)系數(shù)。

2.3.2 加權(quán)幾何平均模型 如果預測因子與蟲草產(chǎn)量的線性函數(shù)為fi，i=1， 2， …， N，N為最終入選模型的因子數(shù)目。各個因子權(quán)重比值為wi，則蟲草產(chǎn)量Y為公式（6）。

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

2.4 蟲草產(chǎn)量數(shù)據(jù)調(diào)研和采集

在2013，2014年6月，通過實地走訪那曲地區(qū)聶榮縣、比如縣、索縣和那曲蟲草交易市場，和采挖蟲草的藏民進行交流，獲取了蟲草的采挖情況。結(jié)合查閱資料獲取了近6年的那曲蟲草產(chǎn)量[1521]，見表2。

2.5 生態(tài)因子與蟲草產(chǎn)量的關(guān)系

2.5.1 雪線高程 根據(jù)蟲草的采挖時間，選擇收集那曲地區(qū)的5月份的landsat衛(wèi)星影像數(shù)據(jù)，如果5月份數(shù)據(jù)含云量特別大，可選擇4月份的數(shù)據(jù)。那曲地區(qū)所在的影像條帶號是：WRS_PATH=137，WRS_ROW=038。最終選擇的數(shù)據(jù)情況見表3。

按照2.1的方法對遙感進行處理，得到雪線處的高程信息見表4，提取雪線的過程中雖然刪除了影像邊緣的異常值，但是影像零散的冰雪覆蓋區(qū)對雪線的邊界有影響，根據(jù)正態(tài)分布原理[12]，把提取出來的雪線高程的“平均值-2×標準差”作為雪線的最低高程。按照2012年、2014年、2011年、2009年、2010年的順序，雪線的最低高程依次上升，見圖1。

由于2013年的影像含云量較大，建模過程中舍去這一年的數(shù)據(jù)。利用剩余的5期數(shù)據(jù)，計算得雪線高程和蟲草產(chǎn)量的相關(guān)系數(shù)r=0.40，線性回歸關(guān)函數(shù)為y=0.030 1x-99.367。

2.5.2 4月和5月降水量之和 選擇那曲地區(qū)附近的5個氣象臺站（編號為55294，55299，55493，56106，56202）的降水量數(shù)據(jù)的平均值作為該因子的值。通過比較單月降水量、各月降水量之和與蟲草

產(chǎn)量的相關(guān)系數(shù)發(fā)現(xiàn)，4月和5月降水量之和比單獨月降水量與蟲草產(chǎn)量的相關(guān)性更大，于是選擇4月和5月降水量之和作為預測因子。那曲地區(qū)各年4月和5月降水量之和見表5。

計算得到4月和5月降水量的和與蟲草產(chǎn)量的相關(guān)系數(shù)r=0.95，線性回歸關(guān)函數(shù)為y=0.13x-49.44。

2.5.3 5月平均氣溫 同降水量一樣，選擇5個氣象臺站的數(shù)據(jù)平均值作為該因子的值。通過比較單月平均氣溫、各月平均氣溫之和與蟲草產(chǎn)量的相關(guān)系數(shù)發(fā)現(xiàn)，5月份的平均氣溫比其他月、月份之和與蟲草產(chǎn)量的相關(guān)性更大，于是選擇5月平均氣溫作為預測因子，那曲地區(qū)各年5月平均氣溫見表5。計算得到5月平均氣溫與蟲草產(chǎn)量的相關(guān)系數(shù)r=-0.69，線性回歸關(guān)函數(shù)為y=-1.477 1x+114.09。

2.5.4 4月和5月日照時長之和 通過比較單月日照時長、各月日照時長之和與蟲草產(chǎn)量的相關(guān)系數(shù)發(fā)現(xiàn)，4月和5月日照時長之和比其他月、月份的和與蟲草產(chǎn)量的相關(guān)性更大，于是選擇4月和5月日照時長之和作為預測因子，那曲地區(qū)各年4月和5月日照時長之和見表5。計算得到5月平均氣溫與蟲草產(chǎn)量的相關(guān)系數(shù)r=-0.50，線性回歸關(guān)函數(shù)為y=-0.022 1x+144.4。

2.6 構(gòu)建蟲草產(chǎn)量預測模型

根據(jù)各個因子與蟲草產(chǎn)量的相關(guān)系數(shù)，得到相關(guān)系數(shù)的比值是r1∶r2∶r3∶r4=0.4∶0.95∶|-0.69|∶|-0.50|=1∶3∶2∶1.于是各因子權(quán)重比值w1∶w

上式中，X1是雪線高程，X2是4月和5月降水量的和，X3是5月平均氣溫，X4是4月和5月日照時長的和，Y是預測的蟲草產(chǎn)量（單位t）。

3 結(jié)果

3.1 模型精度驗證

把2009—2014年的環(huán)境因子數(shù)據(jù)帶入得到的模型中，得到的預測值和實際值見表6。特別地，2013年的遙感影像含云量較大，無法提取雪線高程，于是對預測模型進行修正，去掉f1，只利用剩余的3個因子進行預測，得到2013年的蟲草產(chǎn)量預測值。綜合過去6年的預測結(jié)果和真實值，計算得到模型預測的平均偏差是7.98%，最大偏差是17.84%。

3.2 2015年蟲草產(chǎn)量預測

在2015年6月上旬，根據(jù)生態(tài)因子的構(gòu)建方法，得到X1雪線高程是4 574 m，X2 4月和5月降水量之和是80.72 mm，X3 5月平均氣溫是5.08 ℃，X4 4月和5月日照時長之和由于數(shù)據(jù)暫缺無法獲取，所以去掉f4，只利用前3個因子進行預測，得到2015年那曲地區(qū)的蟲草產(chǎn)量是46.42 t，比去年增產(chǎn)約60%，是2013年產(chǎn)量的85%左右。綜合分析來看，今年高產(chǎn)與今年蟲草生長季降雨量豐沛有較大關(guān)系。

4 討論

本文通過RS/GIS技術(shù)來構(gòu)建氣候環(huán)境因子，把氣候環(huán)境因子與蟲草產(chǎn)量進行線性擬合，配合幾何加權(quán)平均模型，得到了蟲草產(chǎn)量的預測模型。經(jīng)過檢驗，發(fā)現(xiàn)模型預測偏差較小。而且，如果有非重要的生態(tài)因子因為實際情況（天氣或者時間關(guān)系）暫缺時，模型的精度不會受到較大影響，具有較高的穩(wěn)定性和靈活性。

[參考文獻]

[1] 章力建，李兵，胡育驕，等.中國冬蟲夏草資源管理概況[J]. 中國草地學報， 2010， 32（增刊）： 1.

[2] 劉天平，卓嘎，郭健斌. 西藏冬蟲夏草資源可持續(xù)利用探討[J]. 西藏研究， 2010（3）： 114.

[3] 蔡佩云，孫世英. 青海省冬蟲夏草資源持續(xù)開發(fā)利用淺析[J]. 中國草地學報， 2010， 32（增刊）： 6.

[4] 張古忍，吳光國，古德祥，等. 冬蟲夏草資源可持續(xù)利用面臨的問題及其對策[J]. 中國草地學報[J]， 2010， 32（增刊）： 51.

[5] 陳仕江，鐘國躍，馬開森. 冬蟲夏草資源可持續(xù)利用的思考[J]. 中國草地學報， 2010， 32（增刊）： 44.

[6] 楊大榮，彭艷瓊，陳吉岳，等. 中國冬蟲夏草分布格局與環(huán)境變化對其分布的影響[J]. 中國草地學報， 2010，32（增刊）： 22.

[7] 吳慶貴，蘇智先，蘇瑞軍，等. 冬蟲夏草生境選擇主導因子[J]. 廣西植物， 2009，29（3）： 331.

[8] 李芬，吳志豐，徐翠，等. 三江源區(qū)冬蟲夏草資源適宜性空間分布[J]. 生態(tài)學報， 2014，34（5）： 1318.

[9] 李暉，王立輝. 西藏冬蟲夏草生境的氣候[J]. 西藏科技， 2011（6）： 68.

[10] 趙英時. 遙感應用分析原理與方法[M]. 北京： 科學出版社， 2003.

[11] 李亮. 評價中權(quán)系數(shù)理論與方法比較[M]. 上海： 上海交通大學， 2009.

[12] 劉自遠，劉成福. 綜合評價中指標權(quán)重系數(shù)確定方法探討[J]. 中國衛(wèi)生質(zhì)量管理， 2006， 13（2）： 44.

[13] 蔡輝，沈婉蘭. 醫(yī)院綜合效益評價中的權(quán)數(shù)（二）——獨立性權(quán)數(shù)[J]. 中國醫(yī)院統(tǒng)計， 1997， 4（3）：131.

[14] 杜子芳. 抽樣技術(shù)及其應用[M]. 北京：清華大學出版社， 2005.

[15] 李芬，吳志豐，徐翠，等. 三江源區(qū)冬蟲夏草資源適宜性空間分布[J]. 生態(tài)學報， 2014，34（5）： 1318.

[16] 新華網(wǎng). 西藏冬蟲夏草價格進入下行通道[EB/OL]. （20100705） [20160209].http：//news.xinhuanet.com/fortune/201007/05/c_12300368.htm .

[17] 西藏商城藏善堂. 新產(chǎn)冬蟲夏草價格2010年5月走勢分析[EB/OL]. （20100423） [20160209 ].http：//www.xizangshop.com/article138.html.

[18] 冬蟲夏草屋. 同仁堂蟲草價格調(diào)價漲幅達三成[EB/OL]. （20120602） [20160209].http：//www.dcxc5.com/article488.html.

[19] 冬蟲夏草屋. 蟲草產(chǎn)量連年下降，冬蟲夏草資源日趨稀缺[EB/OL]. （20120602）[20160209].http：//www.dcxc5.com/article585.html.

[20] 中國經(jīng)濟網(wǎng). 青海計劃今年采集冬蟲夏草35噸[EB/OL]. （20140529） [20160209].http：//city.ce.cn/ztpd/ysgjccj/ysxwlb/yshdxc/201405/29/t20140529_1653483.shtml.

[21] 中國經(jīng)濟網(wǎng). 2014年新草動態(tài)[EB/OL]. （20140529） [20160209].http：//city.ce.cn/ztpd/ysgjccj/ysxwlb/yshdxc/201405/29/t20140529_1653566.shtml.

[責任編輯 呂冬梅]