陳禹錦 羅喻才 于 芬 楊光耀 張文根

(江西省竹子種質(zhì)資源與利用重點實驗室 江西農(nóng)業(yè)大學(xué) 南昌 330045)

氣候影響植物生長，決定植被類型及其地理分布，同時植被是反映氣候變化最鮮明的標(biāo)志[1]?；诘乩硇畔⑾到y(tǒng)(GIS)技術(shù)和最大熵原理構(gòu)建的MaxEnt模型，可預(yù)測物種在未來氣候環(huán)境下的生境適宜分布區(qū)及其擴張或收縮情景[2-3]。近年來，該模型在群落生態(tài)學(xué)與生物多樣性研究中得到了廣泛關(guān)注與應(yīng)用[4]。例如：氣候變化對物種分布的影響研究，如海南鳳仙花(Impatienshainanensis)[5]、南酸棗(Choerospondiasaxillaris)[6]、小黃花茶(Camellialuteoflora)[7]等；物種入侵監(jiān)測研究，如預(yù)測紅火蟻(Solenopsisinvicta)在中國的適生區(qū)分布[8]；自然保護(hù)區(qū)的規(guī)劃設(shè)計，如白冠長尾雉(Syrmaticusreevesii)保護(hù)區(qū)功能區(qū)的重新劃定[9]等。

毛竹[Phyllostachysedulis(Carriere) J. Houzeau][10]是中國栽培歷史久、種植面積廣、利用價值高的竹種，在材料、食品、康養(yǎng)、綠色環(huán)保等領(lǐng)域均有著巨大的應(yīng)用潛力[11-13]。然而，毛竹入侵降低森林群落中喬木和灌木的豐富度，破壞森林植被，改變土壤養(yǎng)分循環(huán)模式，降低碳儲量和改變土壤微生物結(jié)構(gòu)，在中國南方已成為突出的生態(tài)學(xué)問題[14-16]。目前，關(guān)于毛竹環(huán)境適宜性分布、適生區(qū)預(yù)測的研究報道較少[17-18]。本文首次利用Maxent模型并結(jié)合GIS技術(shù)對氣候變化情景下毛竹的潛在適生區(qū)進(jìn)行了預(yù)測，建立了毛竹生境的變化情景，可為毛竹種植、引種以及入侵防治提供參考。

1 材料與方法

1.1 毛竹地理分布點數(shù)據(jù)

基于中國現(xiàn)有分布區(qū)內(nèi)有代表性的34個毛竹居群[19]，結(jié)合中國數(shù)字植物標(biāo)本館(CVH: http://www.cvh.ac.cn/)、國家標(biāo)本資源共享平臺(http://www.nsii.org.cn/2017/home.php)和記錄有毛竹分布地點的文獻(xiàn)資料，經(jīng)篩選和剔除，獲得毛竹分布點數(shù)據(jù)316條，利用谷歌Earth衛(wèi)星圖拾取經(jīng)緯度坐標(biāo)并導(dǎo)入Excel中，保存為.csv格式。

1.2 環(huán)境變量的采集與處理

1.2.1 環(huán)境變量采集

從Worldclim(www.worldclim.org)獲取當(dāng)前時段(1981—2000年)及未來時段(2021—2040年、2061—2080年)的19個氣候生物變量，分別為：年平均溫(bio1)、晝夜溫差月均值(bio2)、等溫性(bio3)、溫度季節(jié)變化(bio4)、最暖月最高溫度(bio5)、最冷月最低溫度(bio6)、年均溫變化范圍(bio7)、最濕季度平均溫度(bio8)、最干季度平均溫度(bio9)、最暖季度平均溫度(bio10)、最冷季度平均溫度(bio11)、年均降水量(bio12)、最濕月降水量(bio13)、最干月降水量(bio14)、降水量變異系數(shù)(bio15)、最濕季度降水量(bio16)、最干季度降水量(bio17)、最暖季度降水量(bio18)和最冷季度降水量(bio19)。

未來氣候變量是基于聯(lián)合國政府間氣候變化專門委員會(IPCC)發(fā)表的第6次評估報告對全球未來氣候變化所進(jìn)行的預(yù)估而得出的。大氣環(huán)流模型采用中國國家氣候中心開發(fā)的BCC_CSM模型[20]。未來氣候情景包括了4種典型濃度路徑(SSPs)，分別為：SSPs126、SSPs245、SSPs370、SSPs585。本文選用SSPs126、SSPs245、SSPs585分別代表低、中、高溫室氣體排放情景。

土壤因子來自于寒旱區(qū)科學(xué)大數(shù)據(jù)中心(http://data.casnw.net/portal/)的中國土壤數(shù)據(jù)集v1.1(China Soil Map Based Harmonized World Soil Database, v1.1)。土壤質(zhì)地和結(jié)構(gòu)、容重、水分含量、酸堿度、鈣離子含量等因素[21-23]均為限制毛竹生長的因子，因此選用土壤質(zhì)地分類(bio20)、土壤淤泥含量(bio21)、土壤容重(bio22)、土壤基本飽和度(bio23)、土壤酸堿度(bio24)和土壤碳酸鹽含量(bio25)數(shù)據(jù)。

地形因子選用海拔(bio26)、坡度(bio27)和坡向(bio28)數(shù)據(jù)。DEM高程數(shù)據(jù)來自中國科學(xué)院地理空間數(shù)據(jù)云(http://www.gscloud.cn)。高程圖像通過Arcgis 10.2中spatial analyze工具提取，坡度、坡向圖層由高程圖層通過Arcgis 10.2中3D Analyst工具提取。

1.2.2 環(huán)境變量處理

所有環(huán)境變量使用WGS-1984地理坐標(biāo)系，按照來源于國家基礎(chǔ)地理信息系統(tǒng)(http：//nfgis.nsdi.gov.cn/)的中國行政區(qū)劃圖剪裁，空間分辨率調(diào)整為2.5′(22.5 km2)。為避免環(huán)境變量之間的過度擬合[24]，保證模型運行的準(zhǔn)確性，將28個環(huán)境變量及毛竹分布點數(shù)據(jù)導(dǎo)入MaxEnt 3.4.1中進(jìn)行預(yù)模擬，初步得出各環(huán)境變量的貢獻(xiàn)率。使用Arcgis 10.2中的重采樣工具將28個環(huán)境變量與毛竹分布點耦合并導(dǎo)出Excel表，在Spss中采用雙變量相關(guān)并勾選Pearson檢驗進(jìn)行相關(guān)性分析。若2個環(huán)境變量的相關(guān)性≥0.9，則保留預(yù)模擬中貢獻(xiàn)率較大的環(huán)境變量。最終，確定17個環(huán)境變量為主要影響毛竹生長的環(huán)境因子，并進(jìn)行最終預(yù)測得出各環(huán)境變量的貢獻(xiàn)率及置換重要值(表1)。環(huán)境變量統(tǒng)一采用ASCⅡ格式。

表1 環(huán)境變量對毛竹分布預(yù)測的貢獻(xiàn)率及置換重要值Tab.1 The contribution and permutation importance of environmental variables to the predicted distribution of P. edulis.

1.3 毛竹分布區(qū)預(yù)測方法

將毛竹分布點數(shù)據(jù)和篩選出的17個環(huán)境變量分別導(dǎo)入MaxEnt 3.4.1的物種分布(Samples)模塊和環(huán)境變量(Environmental layers)模塊。選用75%的毛竹分布點作為訓(xùn)練集用以建立模型，25%作為測試集用以測試模型，分布值以邏輯斯蒂(Logistic)格式輸出。重復(fù)模擬10次，最終輸出的ASCⅡ格式圖層為10次重復(fù)的平均值。Liu等[25]對現(xiàn)有的13種選擇閾值的方法進(jìn)行過深入研究，本文選用最大特殊敏感性加特殊性(Maximum training sensitivity plus specificity)為該模型的閾值。

1.4 模型精確度評價

使用刀切法(Jackknife)計算各環(huán)境變量對毛竹分布區(qū)域影響的貢獻(xiàn)率。以假陽性(1-特異度)為橫坐標(biāo)，真陽性(靈敏度)為縱坐標(biāo)繪制接受者操作特征曲線(ROC曲線)。模型的優(yōu)劣由ROC曲線下方的面積AUC值決定。AUC的數(shù)值范圍為0～1，數(shù)值越大表示模型準(zhǔn)確度越高。根據(jù)Swets的標(biāo)準(zhǔn)[26]，AUC值為0.5～0.6時模型不存在預(yù)測能力，AUC值為0.6～0.7時模型預(yù)測能力一般，AUC值為0.7～0.8時模型預(yù)測能力良好，而AUC值>0.9時模型準(zhǔn)確性非常高。

1.5 適生等級劃分

選取MaxEnt重復(fù)運算10次所得到的最大特殊敏感性加特殊性平均值0.198為毛竹適生區(qū)分布概率的閾值p，p<0.198為非適生區(qū)，p>0.198為潛在適生區(qū)。將預(yù)測結(jié)果導(dǎo)入Arcmap 10.2中，使用重采樣工具以自然斷點分級法將圖層劃分為非適生區(qū)(p<0.198)、低適生區(qū)(0.1980.6)4個等級，并制作毛竹當(dāng)前及未來適生區(qū)分布圖。

1.6 未來不同氣候背景下毛竹分布區(qū)域的動態(tài)變化

為對比未來氣候情景下毛竹適生區(qū)面積與當(dāng)前適生區(qū)面積的變化趨勢，以0.198作為閾值，指示適生(p>0.198)與非適生(p<0.198)2種狀態(tài)，繪制二進(jìn)制圖層。將當(dāng)前與未來毛竹適生區(qū)的二進(jìn)制圖層導(dǎo)入Arcmap中，進(jìn)行疊加、運算，獲得毛竹適生區(qū)的動態(tài)變化圖。

2 結(jié)果與分析

2.1 模型精確度檢驗

采用ROC曲線下方面積AUC值檢驗?zāi)Ｐ途?。?jīng)過10次運算，當(dāng)前適生區(qū)模型的平均AUC值為0.912(圖1)。根據(jù)Swets的標(biāo)準(zhǔn)可知，AUC值>0.9時模型準(zhǔn)確性非常高，說明模型能夠準(zhǔn)確預(yù)測毛竹的適生區(qū)域。

圖1 毛竹MaxEnt模型的ROC曲線Fig.1 ROC curve of the MaxEnt model in P. edulis

2.2 影響毛竹分布的氣候因子

基于貢獻(xiàn)率、置換重要值和刀切法檢驗并選擇限制毛竹在中國潛在地理分布的環(huán)境變量因子。貢獻(xiàn)率較大的前4個環(huán)境因子變量分別為最干月降水量(bio14，73.3%)、年均降水量(bio12，5.9%)、最冷月最低溫(bio6，5.2%)和溫度年較差(bio7，3.5%)。置換重要值的大小表示模型對該變量依賴性的強弱，置換重要值較大的前4個環(huán)境因子變量為海拔(bio26，29.9%)、最干月降水量(bio14，12.9%)、最冷月最低溫(bio6，11.5%)和年平均氣溫(bio1，9.5%)。根據(jù)刀切法所得數(shù)據(jù)(圖2)可知，訓(xùn)練模型得分較高的前4個環(huán)境因子變量分別為最冷月最低溫(bio6)、最干月降水量(bio14)、年均降水量(bio12)和年平均氣溫(bio1)。

圖2 各環(huán)境變量對毛竹當(dāng)前潛在分布模型的Jackknife檢驗得分Fig.2 Jackknife test scores of environmental variables on the current potential distribution model of P. edulis.

由此可知，主導(dǎo)毛竹潛在地理分布的是降水變量(最干月降水量、年均降水量)、溫度變量(最冷月最低溫、溫度年較差、年平均溫)和地形變量(海拔)。土壤因子雖然貢獻(xiàn)率、置換重要值較小，但對毛竹的生長有不可替代的作用。繪制以上6個環(huán)境變量的響應(yīng)曲線(圖3)以表示毛竹分布概率與環(huán)境變量之間的關(guān)系。當(dāng)毛竹適生區(qū)分布概率大于閾值(0.198)時，限制毛竹分布的主導(dǎo)環(huán)境因子所在區(qū)間分別為最干月降水量(10.0～183.1 mm)、年均降水量(852.4～3 128.6 mm)、最冷月最低溫(-4.7～15.7 ℃)、年平均溫(13.0～23.8 ℃)、溫度年較差(<34.5 ℃)和海拔(<1 315.0 m)。

圖3 毛竹對主要環(huán)境因子的響應(yīng)曲線Fig.3 Response curve of P. edulis to major environmental factors.

2.3 當(dāng)前毛竹潛在分布區(qū)

在當(dāng)前氣候條件下，毛竹實際地理分布點均位于潛在分布區(qū)范圍內(nèi)(圖4)，預(yù)測結(jié)果與實際分布情況相吻合。毛竹潛在適生區(qū)的總面積為191.61萬km2，主要位于亞熱帶季風(fēng)氣候區(qū)。分布范圍為18.31°—37.47° N，93.46°—122.47° E，南起于海南島的南部，北至陜西南部(越過秦嶺，分布于秦嶺北麓)、河南南部及山東沿海一帶；西起于西藏東南部、四川的東部和云南中南部，東至江蘇、浙江、福建沿海一帶及臺灣地區(qū)。其中，高適生區(qū)面積為24.75萬km2，主要位于總適生區(qū)的東部和西北部，即江西中北部、湖北西部、安徽南部、江蘇沿海區(qū)域、重慶西部和四川東南部等地，呈間斷分布；中適生區(qū)面積為89.43萬km2，主要位于總適生區(qū)的東部，即江西南部、福建大部、浙江、安徽中南部，呈連續(xù)分布；低適生區(qū)面積為77.43萬km2，主要位于總適生區(qū)北部、西部和南部的邊緣區(qū)域，呈連續(xù)分布。

圖4 毛竹地理分布點及當(dāng)前潛在適生區(qū)預(yù)測Fig.4 Geographical distribution and prediction of current potential suitable area for P. edulis.

2.4 未來毛竹潛在適生區(qū)及其動態(tài)變化

2.4.1 未來毛竹潛在適生區(qū)面積變化

在BCC_CSM大氣環(huán)流模型、3種溫室氣體排放情景、未來2個時間段條件下，模擬得到6個毛竹未來分布模型(圖5)。與當(dāng)前潛在適生區(qū)面積相比(表2)，未來潛在適生區(qū)面積均有不同程度地縮減，不同等級適生區(qū)面積的變化幅度較大。高適生區(qū)面積在2021—2040年時間段內(nèi)將大幅縮減，到2061—2080年時間段內(nèi)將進(jìn)一步縮減(縮減幅度>90%)，并且高適生區(qū)面積隨溫室氣體排放情景由低到高變化而依次增大。中適生區(qū)面積將在2021—2040年時間段內(nèi)擴張，到2061—2080年時間段內(nèi)中適生區(qū)面積明顯隨溫室氣體排放情景由低到高變化而減小。低適生區(qū)面積將在2021—2040年時間段內(nèi)呈現(xiàn)先擴張后收縮的變化趨勢，在2061—2080年時間段隨溫室氣體排放情景由低到高變化明顯呈現(xiàn)擴張趨勢(擴張幅度>30%)。

圖5 未來氣候情景下的毛竹潛在適生區(qū)Fig.5 Potential suitable area of P. edulis under future climate scenarios

表2 不同氣候變化情景下毛竹適生區(qū)面積變化Tab.2 Suitable area change of P. edulis under different climate change scenarios

2.4.2 毛竹未來潛在適生區(qū)范圍變化

從毛竹潛在適生區(qū)分布圖得知，在2021—2040年時間段內(nèi)，在3種溫室氣體排放情景下，毛竹的中、低適生區(qū)范圍無明顯變化；高適生區(qū)變化范圍隨溫室氣體排放情景由低到高變化依次擴大，無明顯遷移趨勢。而在2061—2080年時間段內(nèi)，廣東、廣西境內(nèi)的低適生區(qū)范圍向北移動，且隨溫室氣體排放情景由低到高變化，低適生區(qū)南邊界向北移動范圍越廣。四川東部、重慶西部、湖南南部、廣東北部、廣西北部的中適生區(qū)范圍隨溫室氣體排放情景由低到高變化而依次減小，中適生區(qū)呈破碎化分布。高適生區(qū)范圍大幅度縮減：在SSPs126模式下，僅在浙江沿海區(qū)域有小范圍分布；在SSPs245模式下，僅在江蘇南部、浙江北部及沿海地區(qū)有小范圍分布；在SSPs585模式下，僅在浙江北部有小范圍分布。從當(dāng)前到2061—2080年時間段，毛竹高適生區(qū)有向東遷移的變化趨勢。

2.4.3 未來毛竹潛在適生區(qū)動態(tài)分布

毛竹潛在適生區(qū)在未來氣候情景下呈現(xiàn)向高緯度地區(qū)轉(zhuǎn)移、整體北移的變化趨勢。由動態(tài)分布圖(圖6)可知，毛竹適生區(qū)的南部除海南島存在面積擴張現(xiàn)象以外，其余地區(qū)面積均縮減，南界整體呈現(xiàn)北移趨勢。適生區(qū)的北部面積向北擴張，北界北移。其中，西北部有極小程度的擴張，東北部(山東、河南)適生區(qū)面積明顯擴張，東、西部適生區(qū)范圍無明顯變化趨勢。

圖6 未來氣候情景下毛竹潛在適生區(qū)動態(tài)變化Fig.6 Dynamic change of potential suitable area for P. edulis under future climate scenarios

2021—2040年時間段內(nèi)SSPs126排放路徑下，以及2061—2080年時間段內(nèi)SSPs245、585排放路徑下，毛竹潛在適生區(qū)在新疆伊寧市西部地區(qū)出現(xiàn)小范圍的擴張。

3 討論

本文基于MaxEnt預(yù)測模型，利用當(dāng)前氣候、土壤、地形數(shù)據(jù)及毛竹分布地點預(yù)測毛竹當(dāng)前潛在適生區(qū)范圍，并在3種溫室氣體排放情景下預(yù)測毛竹在未來2個時間段內(nèi)(即2021—2040年，2061—2080年)潛在適生區(qū)范圍及動態(tài)變化過程。由上文所得結(jié)果可知，毛竹的地理分布主要受到水熱條件的制約，限制毛竹分布最重要的環(huán)境變量為最干月降水量、年降水量、最冷月最低溫、年平均溫、溫度年較差和海拔。黃韜等[27]研究表明，越冬期間的降水天數(shù)及溫度是毛竹育筍的重要限制因子，若越冬期氣溫低、降水量小，則對竹筍萌長不利，從而導(dǎo)致來年毛竹成株數(shù)量減少。該研究與本文預(yù)測結(jié)果相一致，可驗證最干月降水量、最冷月最低溫是限制毛竹地理分布的重要環(huán)境因素。汪陽東等[28]研究認(rèn)為，年降水量對毛竹生長有顯著作用，年平均溫對毛竹生長有極顯著作用，與本文結(jié)果也相一致。海拔通過影響土壤理化性質(zhì)、土壤酶活性、光照條件和氣溫進(jìn)而影響毛竹的生長，從而限制毛竹的垂直地理分布。毛竹在山地地貌條件下分布的海拔上限是1 200～1 500 m[29]，與本文海拔響應(yīng)曲線所得出的上限1315 m大致相同。

毛竹當(dāng)前潛在適生區(qū)范圍為18.31°—37.47° N，93.46°—122.47° E，面積為191.61萬km2。在未來氣候情景下，毛竹適生區(qū)面積部分縮減，高適生區(qū)面積縮減幅度大，適生區(qū)范圍整體呈現(xiàn)北移的趨勢。王佳鑫等[17]所得的當(dāng)前毛竹適生區(qū)面積為246萬km2，差異主要在于未考慮地形變量，因此我們推斷地形變量尤其是海拔變量對毛竹地理分布有明顯的限制作用。與張雷[18]的預(yù)測結(jié)果相比，本文中當(dāng)前潛在毛竹適生區(qū)范圍更廣，在未來氣候情景下北移幅度明顯較小。造成結(jié)果差異的原因主要有：1)2001年中國植被圖所反映的是20世紀(jì)80和90年代的植被分布現(xiàn)狀，而毛竹根系發(fā)達(dá)，入侵速度快，地理分布區(qū)域在較短時間尺度內(nèi)會有較大變化，從2001年中國植被分布圖中提取的毛竹分布數(shù)據(jù)不能很好地反映當(dāng)前毛竹實際地理分布；2)未考慮土壤變量和地形變量對毛竹生長的影響。金佳鑫等[17]研究表明，土壤因素對于毛竹潛在適生區(qū)在未來氣候情景下北移的幅度具有明顯的限制作用。

預(yù)測結(jié)果對開展實際工作具有指示作用。毛竹當(dāng)前潛在地理分布區(qū)主要位于中國東南部，是典型的亞熱帶季風(fēng)氣候區(qū)，該區(qū)域雨量充沛，氣候溫暖，水熱條件良好，與毛竹的生長習(xí)性相適應(yīng)。當(dāng)前實際地理分布點范圍遠(yuǎn)小于毛竹潛在適生區(qū)范圍，并且我國自20世紀(jì)在云南、陜西、山東等地開展的毛竹引種實驗均已取得成功[30-32]，說明在當(dāng)前潛在適生區(qū)內(nèi)進(jìn)行毛竹的引種栽培具有一定的可行性。在未來氣候情景下毛竹適生區(qū)具有北移的可能性，建議根據(jù)本文所得出的毛竹適生區(qū)變化趨勢及適生區(qū)分類等級并結(jié)合具體的社會、經(jīng)濟、環(huán)境因素選擇適宜地區(qū)進(jìn)行毛竹引種栽培。

毛竹是典型的克隆植物。通過竹鞭進(jìn)行無性繁殖向其他林分?jǐn)U張[33]，由闊葉林過渡為竹闊混交林最終成為毛竹純林，威脅其他物種的生存，降低了其生態(tài)功能[34]。目前，毛竹入侵產(chǎn)生的生態(tài)危害已成為熱點研究問題。本研究獲得毛竹當(dāng)前及未來的潛在分布范圍及適生程度，結(jié)合毛竹分布范圍的時間尺度與空間尺度可知：在全球氣候變暖的大前提下，毛竹適生區(qū)范圍有隨時間推移而北移程度加劇的變化趨勢，可為未來適生區(qū)面積擴張的地理區(qū)域提供預(yù)警。應(yīng)在毛竹當(dāng)前潛在適生區(qū)內(nèi)及未來可能擴張的區(qū)域內(nèi)密切監(jiān)測毛竹擴張動向，結(jié)合該區(qū)域的森林植被類型、物種組成、群落結(jié)構(gòu)等特性開展有針對性的管理、清除工作，保護(hù)該地區(qū)生態(tài)功能免受破壞?？傊瑧?yīng)權(quán)衡毛竹引種帶來的經(jīng)濟效益與毛竹擴張造成的生態(tài)危害，豐富毛竹入侵管控手段，控制竹林邊界，局限毛竹生長范圍，實現(xiàn)竹林的精細(xì)化管理，在實現(xiàn)經(jīng)濟價值的同時保持生態(tài)環(huán)境良好、生態(tài)功能穩(wěn)定。