王冶 ，薛忠財 ，楊越 ，常佳寧 ，趙春穎 ，李云峰 ，曾凡然

1.河北民族師范學(xué)院資源與環(huán)境科學(xué)學(xué)院，河北 承德 067000；2.河北省山區(qū)地質(zhì)環(huán)境重點實驗室，河北 承德 067000；3.冀北山區(qū)土壤與特色植物營養(yǎng)創(chuàng)新研究中心，河北 承德 067000；4.承德醫(yī)學(xué)院河北省中藥研究與開發(fā)重點實驗室，河北 承德 067000；5.承德筆架山生態(tài)農(nóng)業(yè)科技發(fā)展有限公司，河北 承德 067000

物種分布區(qū)特征是物種重要的生態(tài)學(xué)和進(jìn)化學(xué)特征[1]，對制定物種保育策略具有重要參考價值[2]，諸多因子均可能影響物種分布區(qū)特征[3-5]，其中氣候變化已成為影響生態(tài)系統(tǒng)和生物類群的最主要因素之一[6]。氣候變化導(dǎo)致我國一些野生藥用植物資源面臨巨大威脅[7-8]，預(yù)測未來氣候變化對物種地理分布的影響，有利于指導(dǎo)瀕危物種保育策略的制定。目前用于預(yù)測物種地理分布的常見模型已較為成熟[9]，其中最大熵（MaxEnt）模型基于最大熵原理，對物種的潛在分布進(jìn)行預(yù)測，因其預(yù)測效果好、預(yù)測覆蓋面積大等特點[10]，被廣泛用于分析未來氣候?qū)ξ锓N分布的影響、野生動植物保護(hù)、物種適宜性等[11]。北蒼術(shù)Atractylodes chinensis(DC.) Koidz.是菊科蒼術(shù)屬多年生草本植物[12]，其根莖是中藥蒼術(shù)的主要來源[13]。燕山地區(qū)所產(chǎn)北蒼術(shù)“朱砂點密，香氣濃郁”，因此作為該藥材的道地產(chǎn)區(qū)[14]。燕山地區(qū)地形地貌、生態(tài)系統(tǒng)相對復(fù)雜，氣候變化對北蒼術(shù)適生區(qū)的潛在影響尚不明確。因此，研究氣候變化影響下的北蒼術(shù)適生區(qū)分布規(guī)律十分必要。

本研究基于實地調(diào)查和查閱文獻(xiàn)獲取燕山地區(qū)北蒼術(shù)地理分布信息，應(yīng)用MaxEnt模型和地理信息系統(tǒng)軟件ArcGIS，分析當(dāng)前燕山地區(qū)北蒼術(shù)適生區(qū)分布，預(yù)測不同氣候情景下未來不同時期潛在適生區(qū)變化，為北蒼術(shù)引種馴化、資源保護(hù)和產(chǎn)業(yè)發(fā)展提供理論依據(jù)。

1 資料與方法

1.1 研究區(qū)概況與樣點分布

本研究區(qū)域為廣義燕山地區(qū)，位于壩上高原以南、河北平原以北、白河谷地以東、遼河平原以西的山地（93°41＇～43°19＇N、116°60＇～112°13＇E），調(diào)查區(qū)域總面積1.3×105km2，包含北京、張家口、唐山、承德、秦皇島、赤峰、葫蘆島、朝陽、錦州等市的全部或部分區(qū)域。燕山地區(qū)處于溫帶大陸性季風(fēng)氣候區(qū)，年均溫8.7～12.3 ℃，年降水量487.6 mm左右，作為京津冀地區(qū)的天然屏障和水源涵養(yǎng)地，海拔梯度較大，涵蓋森林、灌叢、草地、農(nóng)田、濕地等多樣的生境類型。

野生北蒼術(shù)地理分布信息收集：①本團(tuán)隊于2020年8月和2021年7月赴燕山地區(qū)進(jìn)行實地考察，基于主要環(huán)境因子在燕山地區(qū)的空間分布特征進(jìn)行布點，收集北蒼術(shù)地理分布信息共36個；②查閱中國數(shù)字植物標(biāo)本館（https://www.cvh.ac.cn/）、中國國家標(biāo)本資源平臺（http://www.nsii.org.cn/）及相關(guān)文獻(xiàn)[15-17]，經(jīng)過數(shù)據(jù)篩選與核查，收集記錄信息河北12 個、內(nèi)蒙古5個、遼寧4個，利用ArcGIS10.0軟件，基于燕山地區(qū)范圍進(jìn)行分布點信息校對，除去無效記錄，剩余13個。最終確定燕山地區(qū)野生北蒼術(shù)地理分布信息共49 個（見圖1）。本文涉及的地圖來自自然資源部標(biāo)準(zhǔn)地圖服務(wù)系統(tǒng)（http://bzdt.ch.mnr.gov.cn/）。

圖1 燕山地區(qū)北蒼術(shù)樣點分布

1.2 環(huán)境因子數(shù)據(jù)

選擇生物氣候變量、地形地貌變量及土壤屬性變量等30個環(huán)境因子數(shù)據(jù)（見表1），19個溫度、降水相關(guān)的生物氣候變量數(shù)據(jù)[18]來自世界氣候數(shù)據(jù)庫（WorldClim，Version2.0，http://www.worldclim.org），包括當(dāng)前時期（1970-2000 年）及未來3 個時期（2021-2040年、2041-2060年、2061-2080年），并基于聯(lián)合國政府間氣候變化專門委員會第六次評估報告，選擇全球氣候模型預(yù)測的溫室氣體低排放情景（SSP1-2.6）和高排放情景（SSP5-8.5）氣候數(shù)據(jù)。地形地貌變量中數(shù)字高程模型（DEM）數(shù)據(jù)來源于地理空間數(shù)據(jù)云平臺（http://www.gscloud.cn），坡度和坡向數(shù)據(jù)基于DEM，采用ArcGIS10.0軟件表面分析工具計算獲取。土壤屬性變量（0～30 cm）數(shù)據(jù)來自世界土壤屬性數(shù)據(jù)庫（https://geodata.pku.edu.cn/）。以上數(shù)據(jù)空間分辨率均統(tǒng)一為2.5 km。

表1 30個環(huán)境因子

1.3 模型構(gòu)建與精度評價

采用MaxEnt3.4.1 軟件將北蒼術(shù)樣點分布信息（CSV格式）與不同時期的環(huán)境因子（ASCⅡ格式）數(shù)據(jù)輸入軟件中進(jìn)行建模運(yùn)算，刀切法計算各環(huán)境因子的貢獻(xiàn)率，隨機(jī)選擇75%數(shù)據(jù)作為訓(xùn)練集，25%數(shù)據(jù)作為測試集，設(shè)置10次重復(fù)，同時設(shè)定最大迭代次數(shù)為10 000，其他參數(shù)為默認(rèn)值[19]。利用ArcGIS10.0軟件進(jìn)行作圖。采用受試者工作特征（ROC）曲線下面積（AUC）對預(yù)測結(jié)果進(jìn)行評價[20]：AUC≤0.6表明預(yù)測結(jié)果失敗，0.6＜AUC≤0.7表明預(yù)測結(jié)果較差，0.7＜AUC≤0.8表明預(yù)測結(jié)果一般，0.8＜AUC≤0.9表明預(yù)測結(jié)果良好，0.9＜AUC≤1.0表明預(yù)測結(jié)果優(yōu)秀。AUC取值范圍為0～1，越接近1表明模型預(yù)測越準(zhǔn)確[21]。

1.4 適生區(qū)劃分

基于模型結(jié)果，采用ArcGIS10.0軟件重分類工具對適生區(qū)進(jìn)行等級劃分。根據(jù)概率P將適生區(qū)劃分為4個等級：高適生區(qū)（P≥0.8）、中適生區(qū)（0.6≤P＜0.8）、低適生區(qū)（0.4≤P＜0.6）和非適生區(qū)（P＜0.4）[22]。

1.5 時空變化分析

采用重心遷移模型分析不同時期環(huán)境因子變化影響下，燕山地區(qū)北蒼術(shù)適生區(qū)重心遷移軌跡，從而反映其時空格局變化趨勢。重心遷移模型廣泛應(yīng)用在物種遷移、景觀格局演變等空間過程與格局的研究中[7，23-24]，是反映對象空間演變過程的有效方法，其計算公式為：

式中，zi為第i個平面空間單元的AUC值，(xi,yi)為第i個平面空間單元的投影坐標(biāo)值，點(xˉ,yˉ)為北蒼術(shù)適生區(qū)地理分布重心的投影坐標(biāo)值。重心移動的距離計算公式為：

式中，dm表示地理分布重心移動的距離；yt+m和yt分別表示在時期t+m和t的地理分布重心投影y坐標(biāo)；xt+m和xt分別表示在時期t+m和t的地理分布重心投影x坐標(biāo)。

2 結(jié)果與分析

2.1 模型評價

當(dāng)前時期北蒼術(shù)MaxEnt模型AUC平均值為0.913；未來3個時期（2021-2040年、2041-2060年、2061-2080年）不同溫室氣體排放情景下，MaxEnt模型AUC平均值均大于0.9，表明模型預(yù)測結(jié)果優(yōu)秀，可用于北蒼術(shù)適生區(qū)分布的預(yù)測研究。見圖2。

圖2 不同時期北蒼術(shù)生態(tài)分布預(yù)測模型AUC比較

2.2 環(huán)境因子篩選

影響北蒼術(shù)生長的主要環(huán)境因子為溫度季節(jié)性變化（46.6%），其次為高程、降水量季節(jié)性變化、最濕季度降水量、年均降水量，累積貢獻(xiàn)率達(dá)73.4%。主要環(huán)境因子適宜值和貢獻(xiàn)率見表2。

表2 影響北蒼術(shù)分布主要環(huán)境因子適宜值與貢獻(xiàn)率

2.3 適生區(qū)及其動態(tài)變化

2.3.1 當(dāng)前時期適生區(qū)分布

當(dāng)前時期燕山地區(qū)北蒼術(shù)適生區(qū)分布見圖3?？傔m生區(qū)面積約7.54×104km2，占燕山地區(qū)總面積的59.10%，其中高適生區(qū)主要集中在承德西南部、張家口東部、秦皇島西北部和唐山北部，總面積約1.06×104km2，占燕山地區(qū)總面積的8.29%；中適生區(qū)主要集中在承德東部、秦皇島東部和赤峰西南部，總面積約2.12×104km2，占燕山地區(qū)總面積的16.61%；低適生區(qū)主要集中在朝陽西部、葫蘆島西部和北京北部，總面積約4.36×104km2，占燕山地區(qū)總面積的34.19%。綜上表明，燕山地區(qū)北蒼術(shù)適生區(qū)主要以承德、秦皇島和張家口為核心，高、中、低適生區(qū)相連并呈發(fā)散狀向東北方向分布。本研究北蒼術(shù)地理分布信息雖無位于燕山地區(qū)西北部的分布點，但根據(jù)模型預(yù)測結(jié)果可以看出，該區(qū)域也適合北蒼術(shù)生長。

圖3 當(dāng)前時期燕山地區(qū)北蒼術(shù)適生區(qū)分布

2.3.2 未來氣候變化對適生區(qū)分布的影響

未來不同時期2種氣候情景下燕山地區(qū)北蒼術(shù)適生區(qū)分布見圖4，適生區(qū)面積變化見表3。結(jié)果表明，在SSP1-2.6情景下，隨著時間變化，適生區(qū)面積增加的區(qū)域主要集中在朝陽中、西部和葫蘆島南部，而承德北部與赤峰南部適生區(qū)面積減少。2021-2040年適生區(qū)面積減少2 358.40 km2，占4.16%；2041-2060年適生區(qū)面積減少3 849.21 km2，占6.79%；2061-2080年適生區(qū)面積減少5 599.63 km2，占9.88%。在SSP5-8.5情景下，承德中北部與赤峰南部適生區(qū)面積隨時間增加，而朝陽、葫蘆島、秦皇島、唐山地區(qū)適生區(qū)面積均減少。2021-2040年適生區(qū)面積增加1 056.04 km2，占1.86%；2041-2060年適生區(qū)面積增加2 887.32 km2，占5.09%；2061-2080年適生區(qū)面積增加5 610.13 km2，占9.90%。綜上表明，SSP1-2.6情景導(dǎo)致的氣候變化使燕山地區(qū)北蒼術(shù)適生區(qū)面積呈退化趨勢，主要體現(xiàn)為高、中適生區(qū)向低、非適生區(qū)轉(zhuǎn)化，而SSP5-8.5情景導(dǎo)致的氣候變化使燕山地區(qū)北蒼術(shù)適生區(qū)呈擴(kuò)張趨勢，主要體現(xiàn)為低、非適生區(qū)向高、中適生區(qū)轉(zhuǎn)化。

表3 未來不同氣候情景下燕山地區(qū)北蒼術(shù)適生區(qū)面積轉(zhuǎn)移矩陣

圖4 未來不同氣候情景下燕山地區(qū)北蒼術(shù)適生區(qū)變化

2.4 氣候變化對適生區(qū)空間遷移的影響

2.4.1 適生區(qū)水平分布變化

未來2種氣候情景下，燕山地區(qū)北蒼術(shù)重心遷移路線見圖5。結(jié)果表明，當(dāng)前適生區(qū)重心位于承德中部，SSP1-2.6情景下，適生區(qū)重心總體向燕山地區(qū)中、東部遷移，該區(qū)域為內(nèi)蒙古高原與遼河平原的過渡帶，最終遷至承德與朝陽交界處，各時期遷移距離分別為21.02、26.34、36.16 km，總體遷移距離為83.52 km。SSP5-8.5情景下，適生區(qū)重心總體向燕山地區(qū)北部遷移，該區(qū)域為山區(qū)與內(nèi)蒙古高原，最終遷至赤峰南部地區(qū)，各時期遷移距離分別為26.27、27.27、46.94 km，總體遷移距離為100.44 km。

圖5 不同氣候情景下燕山地區(qū)北蒼術(shù)適生區(qū)重心遷移路線

2.4.2 適生區(qū)垂向分布變化

提取2種氣候情景下，不同時期的燕山地區(qū)北蒼術(shù)適生區(qū)海拔數(shù)據(jù)，進(jìn)行高度為300 m的頻率分析。結(jié)果表明，各時期適生區(qū)面積垂向分布規(guī)律一致，隨海拔升高，適生區(qū)面積呈先增加后減少的特征，主要分布海拔區(qū)間為300～900 m。見圖6。

圖6 不同氣候情景下燕山地區(qū)北蒼術(shù)適生區(qū)垂向分布變化

對燕山地區(qū)北蒼術(shù)適生區(qū)垂向分布變化統(tǒng)計后發(fā)現(xiàn)：到2061-2080年，SSP1-2.6情景下，適生區(qū)面積在海拔區(qū)間0～300 m 累積增加30.12%，在海拔區(qū)間300～600、600～900 m，適生區(qū)面積分別累積減少19.06%、7.56%，其他海拔區(qū)間變化不顯著；SSP5-8.5情景下，在海拔區(qū)間0～300 m，適生區(qū)面積累積減少24.96%，在海拔區(qū)間300～600、600～900 m，適生區(qū)面積累積增加10.14%、14.57%，其他海拔區(qū)間變化不顯著?？傮w表明，在SSP1-2.6情景下，燕山地區(qū)北蒼術(shù)適生區(qū)垂向分布呈現(xiàn)由高海拔向低海拔區(qū)域遷移的趨勢，平均海拔下降88 m；在SSP5-8.5情景下，呈現(xiàn)由低海拔向高海拔區(qū)域遷移的趨勢，平均海拔上升119 m。

3 討論

根據(jù)環(huán)境因子貢獻(xiàn)率結(jié)果，對燕山地區(qū)北蒼術(shù)影響最大的氣候變量為溫度季節(jié)性變化，地形地貌變量主要為高程，二者累積貢獻(xiàn)率達(dá)58.6%。溫度季節(jié)性變化在1 100～1 350為最適宜條件，即全年4個季度溫度標(biāo)準(zhǔn)偏差在11～13.5 ℃最有利于北蒼術(shù)生長，與標(biāo)本資料、實地調(diào)查獲取及相關(guān)文獻(xiàn)[25-26]基本一致；北蒼術(shù)最適海拔為300～900 m，此結(jié)果與中國數(shù)字植物標(biāo)本館統(tǒng)計數(shù)據(jù)及《中國植物志》[12]記載較為一致。燕山地區(qū)屬溫帶大陸性季風(fēng)氣候區(qū)，具有鮮明的季度氣溫差異，且由于坡向、海拔等地理因素[27]，氣溫呈現(xiàn)顯著的東西差異，西部地區(qū)溫度季節(jié)性變化為1 050～1 457，東部地區(qū)溫度季節(jié)性變化為890～980，表明西部地區(qū)具有更大的季節(jié)性氣溫變化。且燕山地區(qū)的東、西部地形地貌差異顯著，西部地區(qū)以高原山地為主，海拔在200～2 100 m，東部以丘陵平原為主，海拔30～180 m，總體呈現(xiàn)西高東低的地形特征。北蒼術(shù)適生區(qū)集中在燕山西部地區(qū)，主要環(huán)境因子特征與燕山氣溫變化及地形特征高度吻合，表明模型計算的主要環(huán)境因子結(jié)果合理。

當(dāng)前條件下，燕山地區(qū)北蒼術(shù)適生區(qū)主要集中在承德、秦皇島和張家口東部地區(qū)，與文獻(xiàn)記載[15，28]一致。SSP1-2.6情景下，適生區(qū)分布有所退化，地理分布重心向東遷移，總體東向遷移至海拔較低的內(nèi)蒙古高原與遼河平原過渡帶；SSP5-8.5情景下，適生區(qū)分布有所擴(kuò)張，地理分布重心向東北遷移，總體向更高海拔的內(nèi)蒙古高原方向遷移?？赡苁菧厥覛怏w低排放情景下，低海拔地區(qū)溫度下降，達(dá)到了適生條件，導(dǎo)致北蒼術(shù)向海拔較低的平原方向遷移；溫室氣體高排放情景下，氣溫升高使北蒼術(shù)原適生區(qū)物種豐富度增加，出現(xiàn)生態(tài)位重疊，進(jìn)而產(chǎn)生種間競爭，且高海拔山區(qū)溫度有所升高，導(dǎo)致北蒼術(shù)向海拔更高的高原方向遷移。趙文龍等[7]研究表明，氣候變暖使紅花綠絨蒿低海拔適生區(qū)出現(xiàn)退化，高海拔適生區(qū)面積擴(kuò)大；郭彥龍等[29]研究發(fā)現(xiàn)，氣候變暖使桃兒七適生區(qū)的平均海拔逐漸升高。

本研究顯示，溫度季節(jié)性變化與高程是影響燕山地區(qū)北蒼術(shù)生長的主要環(huán)境因子，當(dāng)前總適生區(qū)面積約7.54×104km2，適生區(qū)較為廣泛。未來隨著氣候變化，北蒼術(shù)分布范圍將出現(xiàn)向不同海拔地區(qū)遷移的現(xiàn)象，因此在北蒼術(shù)人工馴化和引種區(qū)域的海拔選擇上，應(yīng)注意氣候變化，防止適生區(qū)退化。本研究結(jié)果有助于了解氣候變化對燕山地區(qū)北蒼術(shù)適生區(qū)分布的影響，為北蒼術(shù)資源可持續(xù)利用及產(chǎn)業(yè)發(fā)展提供參考。