張世林，高潤紅，高明龍，韓淑敏，張文英，趙 靜

（1.內(nèi)蒙古農(nóng)業(yè)大學(xué) 林學(xué)院，內(nèi)蒙古 呼和浩特 010019；2.包頭市林業(yè)和草原工作站，內(nèi)蒙古 包頭 014030）

全球氣候持續(xù)變暖會導(dǎo)致森林產(chǎn)生不同程度的響應(yīng)?；?050年氣候變化情景，THOMAS 等[1]研究發(fā)現(xiàn)：中等氣候變化背景下，覆蓋地球表面20%樣本區(qū)域?qū)⒂?5%～37%的物種面臨滅絕。物種分布模型(SDM)常用來模擬物種潛在分布及適應(yīng)生境變化[2]。目前，主要的物種分布模型有生物氣候模型(Bioclim)、藥用植物全球產(chǎn)地生態(tài)適宜性區(qū)劃信息系統(tǒng)(GMPGIS)、氣候動態(tài)分析模型(Climex)、規(guī)則集預(yù)測的遺傳算法模型(GARP)和最大熵模型(MaxEnt)[3-5]，其中最大熵模型較其他同類模型具有預(yù)測結(jié)果準(zhǔn)確性高、靈活性強(qiáng)等優(yōu)點(diǎn)[6]。因此，最大熵模型在一些植物的潛在分布區(qū)研究中有廣泛應(yīng)用，如銀杉Cathayaargyrophylla[7]、祁連圓柏Juniperusprzewalskii[8]、沙冬青Ammopiptanthusmongolicus[9]等。

樟子松Pinussylvestrisvar.mongolica作為歐洲赤松P.sylvestris的一個地理變種，屬松科Pinaceae 松屬Pinus，因其較強(qiáng)的耐寒、耐旱、耐貧瘠土壤和較速生等優(yōu)良特性被廣泛用于三北工程等人工林體系的建設(shè)[10-11]。早在1955年，樟子松作為防風(fēng)固沙林首次成功引種到章古臺。迄今為止，樟子松已先后在中國三北工程區(qū)的13 個省(市、自治區(qū))590 余縣引種成功。然而，年平均氣溫達(dá)約7 ℃時，樟子松人工林無法自然更新；因章古臺、毛烏素沙地和塞罕壩等地受降水限制[12-17]，樟子松出現(xiàn)生產(chǎn)力水平降低和枯死等現(xiàn)象。根據(jù)《三北工程六期規(guī)劃(2021—2035年)》，樟子松仍是三北工程中最重要的常綠針葉造林樹種。為防止樟子松引種后出現(xiàn)嚴(yán)重衰退或更新困難等情況，本研究利用優(yōu)化后MaxEnt 模型預(yù)測不同氣候環(huán)境下樟子松潛在分布，以期開展未來樟子松引種和優(yōu)良種質(zhì)資源培育工作。

1 材料與研究方法

1.1 樟子松數(shù)據(jù)來源

2018—2021年內(nèi)蒙古自治區(qū)林業(yè)和草原監(jiān)測規(guī)劃院通過實(shí)地調(diào)查，在內(nèi)蒙古自治區(qū)范圍內(nèi)樟子松分布區(qū)域共獲取130 個分布數(shù)據(jù)；結(jié)合中國數(shù)字植物標(biāo)本館(CVH)、中國國家標(biāo)本植物平臺(NSII)、全球生物多樣性信息網(wǎng)(GBIF) 及其他相關(guān)文獻(xiàn)共獲取樟子松分布點(diǎn)數(shù)據(jù)115 個，經(jīng)檢索匯總后共保留245 條記錄。將文獻(xiàn)記載及實(shí)地調(diào)研數(shù)據(jù)整合，配合ArcGIS 軟件對樣點(diǎn)經(jīng)緯度進(jìn)行篩選剔除。為避免樟子松群體分布得不均勻，以經(jīng)緯度2.5′×2.5′為1 個網(wǎng)格單元，每個單元篩選1 條信息記錄，截至2021年最終保留200 個樣點(diǎn)進(jìn)行后續(xù)擬合分析。

1.2 環(huán)境因子變量及篩選

本研究選取空間分辨率均為30″的19 個生物氣候因子和1 個地形因子作為環(huán)境因子(表1)。氣候數(shù)據(jù)來源于世界氣候數(shù)據(jù)庫(https://www.worldclim.org)，地圖數(shù)據(jù)來源于國家基礎(chǔ)地理信息中心(http://www.ngcc.cn/ngcc/)的中國標(biāo)準(zhǔn)地圖。

表1 全球氣候數(shù)據(jù)庫的環(huán)境因子Table 1 Environmental factors in the Global Climate Database

未來氣候數(shù)據(jù)選取第6 次國際耦合模式計(jì)劃(CMIP6) 的國家氣候中心中等分辨率氣候系統(tǒng)模式(BCC-CSM2-MR)。此模式為中國未來的氣候模式。為更均衡合理模擬未來時期(2050s 和2100s)樟子松潛在分布，同時考慮社會經(jīng)濟(jì)共享路徑(SSP) 和代表性濃度路徑(RCP)，選擇SSP-RCP2.6 (SSP126)、SSP-RCP4.5 (SSP245)和SSP-RCP8.5 (SSP585)等3 個情景對樟子松未來分布區(qū)進(jìn)行預(yù)測。SSP126 情景使用可持續(xù)發(fā)展路徑，代表溫室氣體排放量處于較低水平；SSP245 情景使用中度發(fā)展路徑，代表溫室氣體排放量處于中等水平；SSP585 情景使用以化石燃料為主的發(fā)展路徑，代表溫室氣體排放處于高水平[18]。

為避免20 個生物氣候變量間的相似性過高導(dǎo)致模型過度擬合，通過方差膨脹因子(VIF)和Pearson 相關(guān)性檢驗(yàn)，結(jié)合MaxEnt 模型預(yù)實(shí)驗(yàn)中各個因子貢獻(xiàn)率，選擇相關(guān)系數(shù)小于0.8、VIF 小于10 且對樟子松分布具有較高貢獻(xiàn)率和生態(tài)學(xué)意義的因子。最終保留8 個氣候因子和1 個地形因子進(jìn)行建模，包括年平均氣溫、最熱月最高氣溫、最冷月最低氣溫、年降水量、最冷季度平均氣溫、氣溫季節(jié)變動系數(shù)、最干季度平均氣溫、降水量季節(jié)性變化和海拔。

1.3 模型和方法

本研究利用MaxEnt 模型預(yù)測當(dāng)前和未來時期樟子松在中國的潛在分布。將樟子松200 個分布點(diǎn)和9 個環(huán)境因子導(dǎo)入MaxEnt 3.4.1 軟件，java 環(huán)境下運(yùn)行，隨機(jī)選取25%的樣本作為測試數(shù)據(jù)集，75%的樣本作為訓(xùn)練數(shù)據(jù)集[19]，以Bootstrap 模式，重復(fù)迭代10 次，建模結(jié)果以Logistic 形式輸出。

MaxEnt 模型準(zhǔn)確度可通過受試者工作特征曲線(ROC)進(jìn)行檢驗(yàn)。用該曲線下的面積(AUC)表示所建模型的精度，AUC 通常為0.5～1.0，其值越大，表明模型預(yù)測越精準(zhǔn)。其中AUC 在0.5～0.6 時表明模型預(yù)測失敗，0.6～0.7 時表明模型預(yù)測效果較差，0.7～0.8 時表明模型預(yù)測效果一般，0.8～0.9 時表明模型預(yù)測效果良好，0.9～1.0 時表明預(yù)測效果極好[20]，可以較準(zhǔn)確地反映物種潛在分布。

1.4 模型優(yōu)化

本模型采用優(yōu)化方法為Checkerboard 2 法，調(diào)整調(diào)控倍率(RM)和特征組合(FC) 2 個參數(shù)改變模型正規(guī)化水平。通過改變特征數(shù)量、組合方式以及RM 數(shù)值，將MaxEnt 模型中的5 種特征：線性特征(L)、二次型特征(Q)、片段化特征(H)、乘積型特征(P)和閾值型特征(T)，與16 個RM 值(0.5～8.0，以0.5 為間隔梯度)排列組合出148 種組合方式。利用R 語言中的ENMeval 程序包對其進(jìn)行優(yōu)化測試，結(jié)果中各組合的ΔAICc 值和10%測試遺漏率越低，表示模型預(yù)測精準(zhǔn)度越高[21-22]。

1.5 數(shù)據(jù)的分析和處理

運(yùn)用ArcGIS 將模型運(yùn)行后的數(shù)據(jù)進(jìn)行適宜區(qū)劃分和可視化處理，根據(jù)軟件中提供的多種分類方法預(yù)測樟子松適宜性閾值。根據(jù)冉巧等[7]對孑遺植物銀杉進(jìn)行適宜區(qū)劃分的方法，本研究使用自然斷點(diǎn)法劃分樟子松適宜指數(shù)。參照張偉萍等[8]對祁連圓柏適宜區(qū)劃分的等級，將樟子松生境劃分為不適宜區(qū)、較不適宜區(qū)、一般適宜區(qū)和高度適宜區(qū)(適宜指數(shù)分別為0～0.1、＞0.1～0.2、＞0.2～0.6、＞0.6～1.0)。對比樟子松不同時期的適宜區(qū)分布差異，獲得樟子松不同適宜區(qū)空間分布圖；利用ArcGIS 10.4.1 軟件中SDMtoolbox 2.0 工具包，計(jì)算樟子松在不同時期的適生區(qū)質(zhì)心位置得到樟子松分布區(qū)的遷移方向，并計(jì)算其質(zhì)心遷移距離。

2 結(jié)果與分析

2.1 模型優(yōu)化及準(zhǔn)確性分析

根據(jù)200 個樟子松分布點(diǎn)和9 個環(huán)境因子利用MaxEnt 模型對樟子松的潛在分布區(qū)域進(jìn)行模擬預(yù)測。在MaxEnt 為默認(rèn)參數(shù)設(shè)置時，RM 為1.0，F(xiàn)C 為LQHPT，ΔAICc 為29.77。當(dāng)RM 為1.5，F(xiàn)C 為LQHPT 時，ΔAICc 為0，此時模型最優(yōu)，10% 的訓(xùn)練遺漏率降低4.63% (表2)。因此選取RM 為1.5，F(xiàn)C 為LQHPT 為模型最終參數(shù)。該參數(shù)下10 次模擬訓(xùn)練AUC 均值為0.964 (圖1)，說明預(yù)測結(jié)果精確。

圖1 MaxEnt 模型的ROC 響應(yīng)曲線Figure 1 ROC response curve under MaxEnt model

表2 不同參數(shù)MaxEnt 數(shù)模型優(yōu)化結(jié)果Table 2 Environmental results of MaxEnt model under different parameter settings

2.2 當(dāng)前中國樟子松分布

樟子松當(dāng)前氣候條件下的潛在適生區(qū)集中在53.55°N 以南，135.09°E 以西，總適生區(qū)面積(一般適生區(qū)面積與高度適生區(qū)面積之和)為645 104.2 km2，主要位于大興安嶺地區(qū)、內(nèi)蒙古中部、黑龍江西北部，在山西、河北、吉林、遼寧也有少量分布[23]，其中，高度適生區(qū)面積為147 638.9 km2，一般適生區(qū)面積為497 465.3 km2。樟子松高度適應(yīng)區(qū)絕大部分在內(nèi)蒙古東北部，黑龍江西北部呈零星分布。

2.3 未來氣候情景下樟子松分布預(yù)測

未來時期不同氣候情景下，樟子松適生區(qū)面積較當(dāng)前時期有不同程度減少。2050s，樟子松總適生區(qū)在SSP126、SSP245 和SSP585 氣候情景下降35.7%、52.5%和82.6%；2100s 與2050s 相比，總適生區(qū)面積在SSP126、SSP245 和SSP585 情景下呈下降趨勢，降幅分別為62.1%、37.4% 和100.0%。SSP126-2050s 情景下，樟子松高度適生區(qū)面積相較當(dāng)前大幅下降，面積約為當(dāng)前的6.6%，僅存在大興安嶺較高海拔地區(qū)；SSP245-2100s 氣候情景下，樟子松高度適生區(qū)將完全消失。

在3 種氣候條件下樟子松適生區(qū)面積均下降，SSP245 和SSP585 氣候情景下樟子松適生區(qū)對氣候的響應(yīng)最為明顯 (表3和表4)。SSP126 氣候情景下，2050s 樟子松適生區(qū)面積變化幅度最小，適生區(qū)面積下降230 225.7 km2，降幅為35.7%，主要變化區(qū)域位于呼倫貝爾市西南部；2100s 適生區(qū)面積下降487 968.8 km2，降幅為75.6%，主要在興安盟和錫林郭勒盟的東北部出現(xiàn)較大面積消失。SSP245 氣候情景下，2050s 樟子松適生區(qū)面積下降338 449.4 km2，降幅為52.5%，其變化區(qū)域主要位于內(nèi)蒙古自治區(qū)呼倫貝爾市西東兩側(cè)地區(qū)和赤峰市；2100s 其適生區(qū)面積下降453 211.8 km2，降幅為70.3%，內(nèi)蒙古自治區(qū)東部和黑龍江省的大部分地區(qū)樟子松適生區(qū)消失。SSP585 氣候情景下，樟子松適生區(qū)面積變化幅度最大，2050s 樟子松適生區(qū)面積下降532 882.0 km2，比當(dāng)前減少82.6%，僅在呼倫貝爾市中部有樟子松分布；2100s 樟子松適生區(qū)將在中國完全消失。

表3 中國不同氣候情境下樟子松適生區(qū)面積Table 3 Suitable growing area of P. sylvestris var.mongolica under different climates scenarios in China

表4 不同時期樟子松適生區(qū)空間格局變化Table 4 Dynamic changes in the suitable area for P. sylvestris var.mongolica under different combination of climates cenarios

2.4 質(zhì)心遷移

在空間分布格局方面，未來時期不同氣候情景下樟子松適生區(qū)的遷移位置存在差異，但總體上呈現(xiàn)出向西北和西南方向遷移(表5)。當(dāng)前時期樟子松適生區(qū)的質(zhì)心在內(nèi)蒙古興安盟阿爾山市(47.09°N，120.25°E)。氣候情景為SSP126-2100s 時，樟子松適生區(qū)向西北方向遷移，此時樟子松適生區(qū)質(zhì)心位于內(nèi)蒙古自治區(qū)呼倫貝爾市鄂溫克族自治旗(48.41°N，120.17°E)，遷移距離為161.47 km；氣候情景為SSP585-2050s 時，樟子松向西北方向遷移，此時其適生區(qū)質(zhì)心位于內(nèi)蒙古自治區(qū)呼倫貝爾市新巴爾虎左旗(47.88°N，119.06°E)，遷移距離為126.36 km；氣候情景為SSP245-2100s 時，樟子松向西南方向遷移，此時其適生區(qū)質(zhì)心位于內(nèi)蒙古自治區(qū)錫林郭勒盟東烏珠穆沁旗(45.66°N，117.39°E)，遷移距離為270.80 km。未來時期不同氣候變化情景下，全球增溫增濕使得樟子松在中國適生區(qū)的質(zhì)心整體向西北和西南遷移，質(zhì)心隨著適生區(qū)分布范圍消減可能消失。

表5 不同氣候情景下樟子松適生區(qū)質(zhì)心變化Table 5 Core distributional shifts under different climate scenario/periods for P. sylvestris var. mongolica

2.5 環(huán)境因子分析

根據(jù)Jackknife 法分析結(jié)果(表6)，樟子松潛在地理分布貢獻(xiàn)率超過5%的環(huán)境因子分別為最冷季度平均氣溫(22.1%)、降水量季節(jié)性變化(19.8%)、最冷月最低氣溫(18.1%)、氣溫季節(jié)變動系數(shù)(17.8%)、最干季度平均氣溫(9.0%)、最熱月最高氣溫(6.1%)。這6 個環(huán)境因子總貢獻(xiàn)率達(dá)92.9%。

表6 參與建模的環(huán)境因子貢獻(xiàn)率及置換重要值Table 6 Percentage contribution and permutation importance of environment variables for P. sylvestris var.mongolica in the MaxEnt model

2.6 未來分布區(qū)生態(tài)特征

樟子松分布主要環(huán)境因子與潛在分布區(qū)的關(guān)系如表7所示。SSP126 氣候情景下，2050s 和2100s 樟子松生境適宜度比當(dāng)前降低了0.38 和0.54；SSP245 氣候情景下，2050s 和2100s 樟子松適宜度比當(dāng)前降低了0.46 和0.55；SSP585 氣候情景下，樟子松生境適宜度急劇下降，2050s 下降了0.55，到了2100s 下降0.67，降幅為91%。

表7 樟子松適生區(qū)環(huán)境因子結(jié)果分析Table 7 Results analysis of major environmental variables in P. sylvestris var.mongolica suitable area

200 個樟子松分布點(diǎn)的年平均氣溫與其生境適宜度變化正好相反，2050s，SSP126、SSP245 和SSP585 情景下年平均氣溫分別增加19.5%、56.6% 和101.3%；2100s 與2050s 相比，SSP126、SSP245和SSP585 情景下年平均氣溫分別增加15.9%、48.1%和90.7%。

200 個樟子松分布點(diǎn)的年降水量在SSP126、SSP245 和SSP585 情景下均為遞增趨勢，且年降水量與年平均氣溫變化趨勢相同，都與樟子松物種適宜度變化相反。2050s，SSP126、SSP245 和SSP585 情景下年降水量分別增加5.67、20.32 和44.58 mm；2100s 與2050s 相比，SSP126、SSP245 和SSP585 情景下年降水量分別增加29.58、18.27 和45.62 mm。

3 討論與結(jié)論

3.1 討論

3.1.1 優(yōu)化后模型對樟子松當(dāng)前適生區(qū)預(yù)測 本研究基于氣候和地形因子的研究，應(yīng)用ENMeval 數(shù)據(jù)包中的Checkerboard 2 法來優(yōu)化模型[23]。作為北方主要造林樹種，樟子松集中分布在中國北方地區(qū)，具有喜光、耐寒、耐旱、耐貧瘠和抗病能力強(qiáng)的生長特性。在中國的大興安嶺北部，西起莫爾道戈、金河、根河，東到新林、呼瑪以北有連續(xù)的分布，在南伊圖里河的銀河、免渡河、阿爾山等地呈帶狀或塊狀分布著樟子松純林以及其與興安落葉松Larixgmelinii構(gòu)成的混交林[24]。自1955年樟子松被成功引種到遼寧省阜新市彰武縣章古臺以來，已推廣到科爾沁沙地、渾善達(dá)克沙地、毛烏素沙地等地。本研究基于優(yōu)化的最大熵模型預(yù)測結(jié)果顯示：當(dāng)前樟子松高度適生區(qū)也位于上述地區(qū)，預(yù)測結(jié)果與實(shí)際高度一致，表明經(jīng)過優(yōu)化后的MaxEnt 模型用于樟子松適生區(qū)的預(yù)測結(jié)果精準(zhǔn)可靠。

3.1.2 環(huán)境變量對樟子松地理分布限制 Jackknife 法分析表明：年平均氣溫、最熱月最高氣溫、最冷月最低氣溫、年降水量、最冷季度平均氣溫、氣溫季節(jié)變動系數(shù)、最干季度平均氣溫、降水量季節(jié)性變化和海拔是影響樟子松分布的主要環(huán)境因子。當(dāng)樟子松存在概率大于0.6 時，對應(yīng)的環(huán)境因子最適合植物生長。主導(dǎo)環(huán)境因子分析表明：氣溫和降水相關(guān)因子對樟子松分布貢獻(xiàn)率達(dá)到98.3%，對其分布起決定性作用。這與趙曉彬[25]研究中樟子松適生于年平均氣溫-3.7～-1.5 ℃、年降水量375 mm 的結(jié)論相符。本研究預(yù)測樟子松的最適宜海拔為800 m 左右，與戴繼先等[26]研究發(fā)現(xiàn)樟子松天然分布在300～900 m的陽坡或向陽坡的結(jié)論一致。水熱條件一定程度上決定了樟子松潛在地理分布，地形在一定程度上也會對水熱資源進(jìn)行二次分配，因此，氣溫、降水和地形對樟子松分布格局的相互作用都不能忽視。

3.1.3 氣候變化對樟子松地理分布的影響 2050s 和2100s 不同二氧化碳排放情景下，高溫高濕使樟子松適生區(qū)不斷減少。2100s 樟子松適生區(qū)顯著減少，其降幅明顯大于2050s。未來不同的氣候情景下，樟子松適生區(qū)的空間分布變化隨著氣候變暖加劇，整體向西北高緯度和西南季節(jié)性降水量增加的地區(qū)收縮。SSP126 情景下，2050s 適生區(qū)質(zhì)心遷移距離最短，說明氣候小幅變暖對適生區(qū)分布范圍影響不大。SSP245 情景下樟子松適生區(qū)質(zhì)心點(diǎn)受季節(jié)性降水量影響整體向西南方向降水量充沛地區(qū)遷移，2100s 與2050s 相比其質(zhì)心點(diǎn)遷移距離較大。SSP585 情景下溫室氣體排放增加，2050s 樟子松適生區(qū)質(zhì)心點(diǎn)明顯向西北高緯度地區(qū)遷移，2100s 適生區(qū)在中國消失。張日升等[27]研究表明：隨著溫室氣體排放等環(huán)境問題加劇，年平均氣溫達(dá)7 ℃左右時樟子松將出現(xiàn)生長不良，導(dǎo)致大面積退化。吳祥云等[28]也曾提出：樟子松人工林引種到年平均氣溫為7.1 ℃的章古臺會發(fā)生“早衰”現(xiàn)象，導(dǎo)致面積減退。未來氣候情景下，溫室氣體排放，氣溫升高導(dǎo)致樟子松適生區(qū)面積逐步北移，在中國分布減少甚至至消失。綜上所述，樟子松遷移趨勢與CHEN[29]研究中國寒溫帶森林樹種地理分布在未來全球變暖下往高緯度和降水量充沛地區(qū)遷移的趨勢相符。

3.2 結(jié)論

本研究通過對MaxEnt 模型的優(yōu)化，發(fā)現(xiàn)樟子松潛在分布最優(yōu)模型的特征組合為片段化特征(LQHPT)、乘積型特征(P)和閾值性特征(H)，調(diào)控倍率為1.5。此模型復(fù)雜度低，預(yù)測結(jié)果較準(zhǔn)確，預(yù)測樟子松適生區(qū)范圍與實(shí)際分布情況基本一致。中國樟子松當(dāng)前高度適生區(qū)集中在47.2°～53.3°N，118.8°～121.9°E 的大興安嶺北部高緯度地區(qū)。不同氣候情景下，樟子松對水熱變化最為敏感，特別是在年平均氣溫升高5.0 ℃左右時將導(dǎo)致其在中國的分布區(qū)發(fā)生顯著變化。

通過樟子松不同適生區(qū)的劃分，根據(jù)各自區(qū)劃氣候特點(diǎn)提出以下保護(hù)對策：①在不同地區(qū)選取生長良好的樟子松單株進(jìn)行種質(zhì)資源監(jiān)測、保護(hù)、采集和培育，培育最適合當(dāng)?shù)貧夂驐l件的優(yōu)良品種。②未來在不適合樟子松生長地區(qū)引種時，應(yīng)避免引種樟子松。③提前進(jìn)行系統(tǒng)的種源區(qū)劃研究，培育出不同特性的樟子松并在本地區(qū)推廣。應(yīng)該在年平均氣溫為-3.7～1.6 ℃，降水量為375～425 mm 的東北等高緯度地區(qū)引種樟子松，以發(fā)揮其最大生態(tài)效益。