氣候變化情景下毛竹潛在分布及動態預測

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

(江西省竹子種質資源與利用重點實驗室 江西農業大學 南昌 330045)

氣候影響植物生長，決定植被類型及其地理分布，同時植被是反映氣候變化最鮮明的標志[1]。基于地理信息系統(GIS)技術和最大熵原理構建的MaxEnt模型，可預測物種在未來氣候環境下的生境適宜分布區及其擴張或收縮情景[2-3]。近年來，該模型在群落生態學與生物多樣性研究中得到了廣泛關注與應用[4]。例如：氣候變化對物種分布的影響研究，如海南鳳仙花(Impatienshainanensis)[5]、南酸棗(Choerospondiasaxillaris)[6]、小黃花茶(Camellialuteoflora)[7]等；物種入侵監測研究，如預測紅火蟻(Solenopsisinvicta)在中國的適生區分布[8]；自然保護區的規劃設計，如白冠長尾雉(Syrmaticusreevesii)保護區功能區的重新劃定[9]等。

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

1 材料與方法

1.1 毛竹地理分布點數據

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

1.2 環境變量的采集與處理

1.2.1 環境變量采集

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

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

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

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

1.2.2 環境變量處理

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

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

1.3 毛竹分布區預測方法

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

1.4 模型精確度評價

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

1.5 適生等級劃分

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

1.6 未來不同氣候背景下毛竹分布區域的動態變化

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

2 結果與分析

2.1 模型精確度檢驗

采用ROC曲線下方面積AUC值檢驗模型精度。經過10次運算，當前適生區模型的平均AUC值為0.912(圖1)。根據Swets的標準可知，AUC值>0.9時模型準確性非常高，說明模型能夠準確預測毛竹的適生區域。

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

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

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

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

由此可知，主導毛竹潛在地理分布的是降水變量(最干月降水量、年均降水量)、溫度變量(最冷月最低溫、溫度年較差、年平均溫)和地形變量(海拔)。土壤因子雖然貢獻率、置換重要值較小，但對毛竹的生長有不可替代的作用。繪制以上6個環境變量的響應曲線(圖3)以表示毛竹分布概率與環境變量之間的關系。當毛竹適生區分布概率大于閾值(0.198)時，限制毛竹分布的主導環境因子所在區間分別為最干月降水量(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 毛竹對主要環境因子的響應曲線Fig.3 Response curve of P. edulis to major environmental factors.

2.3 當前毛竹潛在分布區

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

圖4 毛竹地理分布點及當前潛在適生區預測Fig.4 Geographical distribution and prediction of current potential suitable area for P. edulis.

2.4 未來毛竹潛在適生區及其動態變化

2.4.1 未來毛竹潛在適生區面積變化

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

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

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

2.4.2 毛竹未來潛在適生區范圍變化

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

2.4.3 未來毛竹潛在適生區動態分布

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

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

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

3 討論

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

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

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

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