基于Maxent模型對尼泊爾芒適生區(qū)時空分布的預(yù)測

李世成，易自力，廖劍鋒，肖亮

(湖南農(nóng)業(yè)大學(xué)生物科學(xué)技術(shù)學(xué)院，湖南 長沙 410128)

尼泊爾芒為多年生禾本科芒屬植物[1]，分布于中國的西藏中部、云南西北部和四川西部等海拔高度為1 900～2 800 m的山坡或河谷漫灘草地[2]。尼泊爾芒具有多種經(jīng)濟(jì)價值，可用于能源、飼料、水土保持和園林造景等[3–4]，也可為芒屬其他種的遺傳改良提供優(yōu)良遺傳資源[5]。為開發(fā)和利用尼泊爾芒，本課題組曾將其引種到湖南農(nóng)業(yè)大學(xué)芒屬植物資源圃內(nèi)，但發(fā)現(xiàn)種植于該資源圃內(nèi)的尼泊爾芒在自然條件下逐漸消亡，推測可能是由于氣候原因造成的。查閱相關(guān)文獻(xiàn)發(fā)現(xiàn)，目前對尼泊爾芒的研究僅局限于系統(tǒng)分類和地理分布[5–7]，少有利用物種分布模型對其主導(dǎo)環(huán)境因子和具體適生區(qū)分布相關(guān)方面的研究。本研究中，利用 Maxent模型對尼泊爾芒適生區(qū)分布進(jìn)行了研究，以期為尼泊爾芒的引種和資源保護(hù)提供理論依據(jù)。

Maxent模型是一種以最大熵原理和機(jī)器學(xué)習(xí)為基礎(chǔ)的物種分布模型，它可以通過機(jī)器學(xué)習(xí)的方法，根據(jù)物種分布記錄和環(huán)境變量圖層對物種的適生區(qū)進(jìn)行模擬預(yù)測[8–9]。因Maxent模型具有操作簡易、預(yù)測效果好[10–11]等優(yōu)點(diǎn)，在物種分布研究領(lǐng)域得到了廣泛應(yīng)用。另外，利用 Maxent模型所構(gòu)建的模型能夠轉(zhuǎn)移到不同時期或地理空間中，可為預(yù)測物種其他時期的適生區(qū)分布[12]和評估物種入侵風(fēng)險[13]提供工具。但當(dāng) Maxent模型轉(zhuǎn)移到非模擬情景時，其默認(rèn)設(shè)置可能會導(dǎo)致模型過擬合[14]，影響預(yù)測結(jié)果。對此，MUSCARELLA等[15]基于R語言開發(fā)了ENMeval數(shù)據(jù)包，可優(yōu)化Maxent模型默認(rèn)參數(shù)的設(shè)置，以平衡模型的擬合度和復(fù)雜度。

本研究中，通過ENMeval數(shù)據(jù)包對Maxent模型進(jìn)行了優(yōu)化，利用優(yōu)化后的 Maxent模型對現(xiàn)代尼泊爾芒的適生區(qū)進(jìn)行了預(yù)測。根據(jù)貢獻(xiàn)率、置換重要性和Jackknife檢驗(yàn)對限制尼泊爾芒適生區(qū)分布的主導(dǎo)氣候因子進(jìn)行綜合評估，以期為科學(xué)地引種尼泊爾芒提供依據(jù)，并對末次冰盛期、全新世中期和未來(2061～2080年)3個時期的尼泊爾芒適生區(qū)加以預(yù)測，與現(xiàn)代適生區(qū)的分布進(jìn)行對比，以期了解限制現(xiàn)代尼泊爾芒在不同地點(diǎn)分布的主導(dǎo)氣候因子。

1 材料與方法

1.1 物種分布數(shù)據(jù)的獲取與處理

通過查閱全球生物多樣性信息網(wǎng)(https://www.gbif.org/)、中國數(shù)字植物標(biāo)本館平臺(http://www.cvh.ac.cn/)、Kingdonia(http://kun.kingdonia.org/)、中國國家標(biāo)本平臺(http://www.nsii.org.cn/2017/home.php)等網(wǎng)站，獲得 60個尼泊爾芒分布點(diǎn)和研究范圍內(nèi)作為背景點(diǎn)的2 709個禾本科植物的分布點(diǎn)，采用區(qū)塊法將60個尼泊爾芒分布點(diǎn)劃分為3個訓(xùn)練集和1個測試集。

1.2 氣候變量的獲取與篩選

末次冰盛期、全新世中期、現(xiàn)代(1970～2000年)和未來(2061～2080年)的氣候數(shù)據(jù)均來自于Worldclim(http://www.worldclim.org/)，每個時期的氣候數(shù)據(jù)包括19個生物氣候變量，空間分辨率為5弧分；通過ArcGIS對19個生物氣候變量圖層進(jìn)行剪裁。利用19個生物氣候變量在Maxent中建模，根據(jù)各個生物氣候變量的貢獻(xiàn)率和其在SPSS軟件中的相關(guān)性分析結(jié)果[16]，最終篩選出9個生物氣候變量(表1)， 作為預(yù)測尼泊爾芒適生區(qū)的環(huán)境變量。

表1 氣候變量及其重要性參數(shù)Table 1 Environmental variables and corresponding important parameters

1.3 模型的建立與評估

利用ENMeval數(shù)據(jù)包對Maxent模型的設(shè)置參數(shù)中的調(diào)控倍頻(regularization multiplier，RM)和特征組合(feature combination，F(xiàn)C)進(jìn)行優(yōu)化，設(shè)置RM的變化范圍為0.5～4.0，每次增加0.5；FC選取6種組合，分別為L、LQ、H、LQH、LQHP和LQHPT，其中L為線性(linear)，Q為二次型(quadratic)，H為片段化(hinge)，P為乘積型(product)，T為閾值型(threshold)。運(yùn)用 Akaike信息量準(zhǔn)則(Akaike information criterion correction，AICc)評估不同參數(shù)組合的擬合度和復(fù)雜度[17]；運(yùn)用測試受試者工作特征曲線下的面積(area under receiver operating characteristic curve，AUCTEST)[18]評估不同參數(shù)組合的區(qū)分測試點(diǎn)和背景點(diǎn)的能力；利用訓(xùn)練AUC和測試AUC值之差(the difference between training and testing AUC，AUCDIFF)[18]、10%訓(xùn)練遺漏率(10%training omission rate，OR10)[19]與“最小訓(xùn)練集”遺漏率 (‘Minimum Training Presence’ omission rate ，ORMTP)[14]評估不同參數(shù)組合的擬合度，選擇出最適合尼泊爾芒的參數(shù)組合。利用該參數(shù)組合在Maxent模型中進(jìn)行建模和預(yù)測，采用受試者工作曲線(ROC曲線)評估模型的準(zhǔn)確性[18]。在輸出文件中，選擇10次重復(fù)的平均值作為本研究的預(yù)測結(jié)果，將預(yù)測結(jié)果導(dǎo)入到ArcGIS中，根據(jù)Maxent模型中的邏輯輸出結(jié)果(P值)，參考IPCC(聯(lián)合國政府間氣候變化專門委員會)第5次評估報告中對可能性的表述，將整個區(qū)域劃分為4個等級，分別為不適宜區(qū)(0≤P＜0.3)、低適宜區(qū)(0.3≤P＜0.5)、中適宜區(qū)(0.5≤P＜0.7)和高適宜區(qū)(0.7≤P≤1.0)。

1.4 多元環(huán)境相似度面和最不相似變量分析

采用多元環(huán)境相似度面和最不相似變量對不同時期的氣候差異和引起尼泊爾芒適生區(qū)發(fā)生變化的主要?dú)夂蜃兞窟M(jìn)行探討。根據(jù)多元環(huán)境相似度面計算不同情景中的環(huán)境變量與參考圖層中環(huán)境變量的相似度(similarity，S)，當(dāng)S值為正值時，S值越小，表示該點(diǎn)的環(huán)境變量與參考圖層的差異越大，S值為100時表示不存在差異；當(dāng)S值為負(fù)值時，表示該點(diǎn)至少有1個環(huán)境變量的數(shù)值超出了參考圖層的范圍[20]。本研究中，以中國境內(nèi) N16.3°～35.6°，E80.3°～122.5°的主導(dǎo)氣候變量為參考圖層，采用Maxent模型中的density.tools.Novel工具進(jìn)行分析。

2 結(jié)果與分析

2.1 模型的優(yōu)化與準(zhǔn)確性評價

通過ENMeval建立的48種不同Maxent參數(shù)設(shè)置組合的ΔAICc(delta AICc)值如圖1所示。在默認(rèn)情況下，60個尼泊爾芒分布點(diǎn)的默認(rèn)設(shè)置參數(shù)RM=1，F(xiàn)C選取LQH，其ΔAICc值為60.08。ΔAICc值遠(yuǎn)大于 2，表明基于默認(rèn)參數(shù)條件下建立的Maxent模型是不可靠的，要選擇其他的設(shè)置參數(shù)進(jìn)行建模。在 48種參數(shù)組合中，ΔAICc＜2的參數(shù)組合共有3組，分別為：RM=2，F(xiàn)C選取LQH；RM=2.5，F(xiàn)C選取LQH和RM=3，F(xiàn)C選取LQH。它們的ΔAICc分別為0、0.12和0.17，3組參數(shù)組合的ΔAICc都在接受范圍內(nèi)，其復(fù)雜度和擬合度均較低，具有較高的可靠性。進(jìn)一步對ENMeval生成的其他評估指標(biāo)進(jìn)行比較，結(jié)果(表2)發(fā)現(xiàn)：當(dāng)RM=3，F(xiàn)C選取LQH時，其 AUCDIFF、ORMTP和 OR10值均低于其他2個參數(shù)組合的值，表明RM=3，F(xiàn)C選取LQH的擬合度小于其他 2個參數(shù)組合的，并且其AUCTEST的值更高，表明其能更好地區(qū)分測試分布點(diǎn)和背景點(diǎn)，因此，以RM=3，F(xiàn)C選取LQH作為最終的設(shè)置參數(shù)。在該設(shè)置參數(shù)條件下，Maxent模型模擬尼泊爾芒適生區(qū)的訓(xùn)練AUC值為(0.966±0.003)，測試AUC值為(0.963±0.028)，表明模型預(yù)測結(jié)果準(zhǔn)確，可用于對尼泊爾芒適生區(qū)的預(yù)測。

圖1 不同設(shè)置參數(shù)下的Maxent模型的ΔAICcFig.1 ΔAICc of Maxent models under different setting parameters

表2 ΔAICc<2的參數(shù)設(shè)置下的Maxent模型評估結(jié)果Table 2 Evaluation results of Maxent model under settings which is less than 2

2.2 主導(dǎo)氣候變量及其適宜區(qū)間

根據(jù)氣候變量對Maxent模型的貢獻(xiàn)率(表1)，共有5個氣候變量對模型的貢獻(xiàn)率大于5%，分別為 bio12、bio4、bio10、bio6和 bio7，其累計貢獻(xiàn)率達(dá)96.5%，其中bio6和bio7的置換重要性較高，二者累計置換重要性達(dá)95.2%。為了減少氣候變量之間的影響，以進(jìn)一步解釋哪些氣候變量在模型中是最重要的，對 9個氣候變量的重要性進(jìn)行Jackknife法檢測(圖2)。僅單獨(dú)使用氣候變量時， 正則化訓(xùn)練增益、測試增益和AUC值較高的3個氣候變量均為 bio4、bio7和 bio12；在使用正則化訓(xùn)練增益的Jackknife中，單獨(dú)使用bio10時，幾乎沒有任何增益，表明它對預(yù)測尼泊爾芒適生區(qū)的幫助很小；在使用除bio6以外的氣候變量時，正則化訓(xùn)練增益下降的最多，表明它對訓(xùn)練數(shù)據(jù)有更好的匹配性，具有較多其他變量中不存在的信息；在使用測試增益和AUC值的Jackknife中，使用除bio7以外的氣候變量增益下降的最多，這與bio6和bio7的置換重要性結(jié)果相一致。綜上所述，bio4、bio6、bio7和bio12是影響尼泊爾芒分布的主導(dǎo)氣候變量。根據(jù) Maxent模型中生成的響應(yīng)曲線，按照邏輯值大于0.5的范圍，劃分出氣候變量的適宜區(qū)間(表1)。4個主導(dǎo)氣候變量bio4、bio6、bio7和bio12的適宜區(qū)間分別為 2.61～5.75 °C，–5.43～3.71 °C，18.72～27.00 °C 和 864.58～1 468.75 mm。

圖2 Maxent模型中氣候變量重要性的Jackknife檢測結(jié)果Fig.2 Jackknife test result of climatic variables importance in Maxent model

2.3 尼泊爾芒的現(xiàn)代適生區(qū)

60個尼泊爾芒分布點(diǎn)中屬于高適宜區(qū)、中適宜區(qū)、低適宜區(qū)和不適宜區(qū)的比例分別為 13.33%、61.67%、16.67%和8.33%(圖3)，預(yù)測結(jié)果與其分布數(shù)據(jù)基本一致，表明了預(yù)測結(jié)果的準(zhǔn)確性和Maxent模型對評價尼泊爾芒適生區(qū)的實(shí)用性。尼泊爾芒在中國的現(xiàn)代適生區(qū)總面積約為60.92萬km2，其中高適宜區(qū)的面積為6.08萬km2，占適生區(qū)總面積的9.98%。尼泊爾芒的中、高適宜區(qū)主要分布于中國的西藏東部、四川南部、貴州西部和云南大部分地區(qū)，涵蓋了東亞植物區(qū)的8個亞地區(qū)和古熱帶植物區(qū)中的2個亞地區(qū)。這些區(qū)域多以河谷以及高原和山地為主，海拔高度差異大，局部氣候比較復(fù)雜，分布有大量孑遺物種并常為孤立類型，具有熱帶、亞熱帶植物區(qū)系的特征。

2.4 尼泊爾芒適生區(qū)時空變化分析

末次冰盛期到全新世中期期間，尼泊爾芒適生區(qū)大面積消失，適生區(qū)面積減少了29.34%，表現(xiàn)為東海和南海沿海區(qū)域內(nèi)的適生區(qū)消失，內(nèi)陸適生區(qū)的向內(nèi)收縮。全新世中期到現(xiàn)代期間，兩者適生區(qū)輪廓和面積基本相似，但現(xiàn)代的高適宜區(qū)面積與全新世中期相比增加了43.51%，主要表現(xiàn)為云南省內(nèi)的高適宜區(qū)向四周擴(kuò)散。與現(xiàn)代相比，尼泊爾芒適生區(qū)在未來氣候變暖的情景(RCP8.5)下可能會整體北移，海南省和臺灣省內(nèi)的適生區(qū)消失，云南省南部的適生區(qū)破碎化，湖北省將成為新的適生區(qū)，適生區(qū)面積與現(xiàn)代的相比會略有增加，大約增加了5.50%，其中高適宜區(qū)面積約增加了3.18萬km2。

圖3 尼泊爾芒的現(xiàn)代分布點(diǎn)(黑色點(diǎn))和不同時期適生區(qū)Fig.3 Extant occurrence points (black dots) and suitable areas of M.nepalensis during different periods

2.5 多元環(huán)境相似度面和最不相似變量分析

3個時期與現(xiàn)代的多元環(huán)境相似度面和最不相似變量如圖4所示。末次冰盛期、全新世中期和未來適生區(qū)內(nèi)的多元環(huán)境相似度分別為(25.94%±19.57%)、(28.78%±20.48%)和(27.05%±15.91%)，末次冰盛期適生區(qū)內(nèi)的主導(dǎo)氣候變量與現(xiàn)代的差異最大。末次冰盛期所特有的適生區(qū)內(nèi)的相似度為(35.26%±22.60%)，其在東海和南海沿海適生區(qū)內(nèi)的相似度較低，該區(qū)域內(nèi)的最不相似變量主要為bio12，內(nèi)陸適生區(qū)內(nèi)的相似度較高，該區(qū)域內(nèi)的最不相似變量主要為bio7；全新世中期所特有的適生區(qū)內(nèi)的相似度為(35.90%±23.72%)，除 bio7外的 3個主導(dǎo)氣候變量均占據(jù)了一定區(qū)域內(nèi)的最不相似變量；未來所特有的適生區(qū)內(nèi)的相似度為(41.08%±13.54%)，該區(qū)域內(nèi)的最不相似變量主要為 bio4?，F(xiàn)代與末次冰盛期和全新世中期相比，特有適生區(qū)位置相似，相似度分別為(45.95%±19.10%)和(54.98%±13.84%)，主要的最不相似變量均為bio4，現(xiàn)代與未來相比的特有適生區(qū)范圍最大，相似度為(11.70%±8.09%)，最不相似變量主要為bio4和bio6。

3 結(jié)論與討論

Maxent模型的默認(rèn)設(shè)置參數(shù)源于對6個不同地理區(qū)域中266個物種的測試[21]，在為多個物種同時構(gòu)建模型時能夠較好地辨別物種的存在概率，但當(dāng)研究只集中于某一特定物種時，以默認(rèn)設(shè)置參數(shù)構(gòu)建的模型可能不是最優(yōu)的[22]，可能會造成模型的擬合度和復(fù)雜度過高；因此，需要根據(jù)所研究物種調(diào)整設(shè)置參數(shù)。本研究中，利用ENMeval數(shù)據(jù)包計算了多個評估指標(biāo)，發(fā)現(xiàn)默認(rèn)設(shè)置參數(shù)的ΔAICc值遠(yuǎn)遠(yuǎn)大于 2，如果繼續(xù)采用默認(rèn)參數(shù)，則會降低模型轉(zhuǎn)移能力[23]，進(jìn)而降低預(yù)測結(jié)果的準(zhǔn)確度，而設(shè)置參數(shù)組合(RM=3，F(xiàn)C選取LQH)不僅復(fù)雜度和擬合度低，且能較好地區(qū)分測試分布點(diǎn)和背景點(diǎn)；因此，選取了此參數(shù)組合，基于它所構(gòu)建的模型轉(zhuǎn)移能力強(qiáng)，能更準(zhǔn)確地預(yù)測不同時期的尼泊爾芒適生區(qū)。

圖4 尼泊爾芒不同時期的多元環(huán)境相似度面(MESS)和最不相似變量(MoD)Fig.4 Multivariate environmental similarity surface and most dissimilar variable analysis for M.nepalensis during different periods

在利用Maxent模型對物種適生區(qū)進(jìn)行預(yù)測時，常用 19個生物氣候變量作為環(huán)境因子，但是它們之間存在著高度的相關(guān)性，使得模型用于轉(zhuǎn)移時變得不穩(wěn)定[24]；因此，在建模前需要對其進(jìn)行適當(dāng)?shù)暮Y選，這不僅可以降低共線性和過擬合的風(fēng)險，而且也可以使模型簡化和更易解釋。本研究中，綜合分析了每個氣候變量的貢獻(xiàn)率、置換重要性和Jackknife后，發(fā)現(xiàn)溫度季節(jié)性標(biāo)準(zhǔn)差、最冷月最低溫度、年溫度變化范圍和年降水量是影響尼泊爾芒分布的主導(dǎo)氣候因子。利用多元環(huán)境相似度面和最不相似變量分析后發(fā)現(xiàn)，在末次冰盛期，尼泊爾芒適生區(qū)內(nèi)的4個主導(dǎo)氣候因子與現(xiàn)代的差異最大，未來次之，全新世中期與現(xiàn)代最為接近。

末次冰盛期的氣候與現(xiàn)代的反差極大，尼泊爾芒的適生區(qū)面積在該時期最大，中國這一時期的溫度比現(xiàn)代低 5～11°C，降水量是現(xiàn)代的 30%～80%[25]。尼泊爾芒作為一種適應(yīng)于極端或特殊低溫和干旱的植物，末次冰盛期的低溫和降雨量可能會增加尼泊爾芒的適生區(qū)面積。相對于現(xiàn)代，末次冰盛期時的夏季風(fēng)降水和夏季風(fēng)降水強(qiáng)度平均減幅分別為25.1%和14.3%[26]，海平面比現(xiàn)代的低100～150 m，臺灣和海南都與大陸相連[25]，這可能是尼泊爾芒在廣東省、福建省、海南省以及臺灣大面積出現(xiàn)適生區(qū)的原因。全新世中期是歷史上最近的一個暖期，Maxent模型預(yù)測結(jié)果表明，該時期的尼泊爾芒適生區(qū)面積和輪廓與現(xiàn)代相似，沒有明顯的移動趨勢，這可能是因?yàn)槿率乐衅诘臍夂蛱卣髋c現(xiàn)代十分相似[27]，所以對尼泊爾芒適生區(qū)的面積和輪廓影響不大。在未來氣候變暖的情景(RCP8.5)下，Maxent模型預(yù)測的尼泊爾芒適生區(qū)可能會整體北移，海南省和臺灣的適生區(qū)消失，云南省南部的適生區(qū)呈破碎化。綜合4個時期的適生區(qū)發(fā)現(xiàn)，尼泊爾芒在喜馬拉雅—橫斷山脈的適生區(qū)變化不大，可為尼泊爾芒提供長久的氣候避難所，這應(yīng)該與該地區(qū)第四紀(jì)以來相對穩(wěn)定的氣候環(huán)境和復(fù)雜的地形地貌有關(guān)。另外，預(yù)測結(jié)果表明，湖南省屬于尼泊爾芒的不適宜區(qū)，與之前作出的尼泊爾芒在引種湖南省后于自然條件下逐漸消亡是由于氣候原因造成的推測相一致。