基于MaxEnt模型的怒江金絲猴潛在適生區(qū)預(yù)測(cè)

肖 祺，劉 成，2，徐夢(mèng)珍，周雄冬

（1.中國(guó)水利水電科學(xué)研究院，北京 100038；2.國(guó)際泥沙研究培訓(xùn)中心，北京 100048；3.清華大學(xué) 土木水利學(xué)院，北京 100084）

1 研究背景

怒江金絲猴（Rhinopithecus strykeri）是我國(guó)一級(jí)重點(diǎn)保護(hù)野生動(dòng)物，2010年被首次發(fā)現(xiàn)，是目前世界上記錄的第五種金絲猴，被世界自然保護(hù)聯(lián)盟（IUCN）列為極危物種。怒江金絲猴首先被發(fā)現(xiàn)于緬甸，此后由中國(guó)科學(xué)家在高黎貢山地區(qū)發(fā)現(xiàn)其蹤跡。種群數(shù)量?jī)H幾百只，生活在1700 m至3100 m的中山濕性常綠闊葉林和溫涼性針葉林中。國(guó)內(nèi)對(duì)于仰鼻猴屬的研究主要集中于早年發(fā)現(xiàn)的川金絲猴和滇金絲猴，怒江金絲猴作為新發(fā)現(xiàn)物種，其相關(guān)成果多局限于生物學(xué)描述及種群動(dòng)態(tài)研究[1-3]，而對(duì)其適生棲息地分布格局及變化規(guī)律缺乏科學(xué)研究。金絲猴生活在橫斷山區(qū)西緣，該地區(qū)地形氣候條件復(fù)雜，是全球生物多樣性熱點(diǎn)區(qū)域之一[4]，同時(shí)此地還存在的水能開(kāi)發(fā)對(duì)生態(tài)的影響問(wèn)題[5]，對(duì)其進(jìn)行適生棲息地分布研究對(duì)于保護(hù)我國(guó)生物多樣性具有重要意義。近年來(lái)使用生態(tài)位模型進(jìn)行生物多樣性保護(hù)、入侵生物潛在分布區(qū)的研究成為新的熱點(diǎn)[6-8]。常見(jiàn)的生態(tài)位模型有MaxEnt、Domain、Bioclim、ENFA和GARP等。王國(guó)崢等[9]使用4種生態(tài)位模型對(duì)我國(guó)特有孑遺植物金錢(qián)松的潛在適生區(qū)進(jìn)行預(yù)測(cè)，發(fā)現(xiàn)MaxEnt模型結(jié)果更精確；曹向鋒等[10]使用5種生態(tài)位模型對(duì)入侵生物黃頂菊的潛在適生區(qū)進(jìn)行預(yù)測(cè)，篩選后采納了MaxEnt的預(yù)測(cè)結(jié)果；段義忠等[11]在不同時(shí)段使用3種生態(tài)位模型預(yù)測(cè)，結(jié)果也表明MaxEnt模型的最適生預(yù)測(cè)與現(xiàn)實(shí)分布保持一致。本研究使用MaxEnt模型對(duì)怒江金絲猴現(xiàn)代以及未來(lái)氣候情景下的適生區(qū)進(jìn)行研究，為深入了解這一物種的時(shí)空分布格局、進(jìn)行生態(tài)適生區(qū)保護(hù)提供理論依據(jù)。

2 研究方法

2.1 模型分析

2.1.1 模型含義 MaxEnt最大熵模型（maximum entropy model）是以最大熵理論（the theory of maximum entropy）為基礎(chǔ)的物種地理尺度空間分布模型，是目前應(yīng)用最廣、表現(xiàn)最好的生態(tài)位模型之一[12]，主要用于預(yù)測(cè)物種分布區(qū)。本研究中使用開(kāi)發(fā)者Steven Phillips提供的3.4.1版本。

2.1.2 模型參數(shù) MaxEnt模型是一種基于機(jī)器學(xué)習(xí)的復(fù)雜模型，準(zhǔn)確率高但可轉(zhuǎn)移性差，大多數(shù)研究采用默認(rèn)參數(shù)來(lái)構(gòu)建模型，默認(rèn)參數(shù)的設(shè)置來(lái)源于模型開(kāi)發(fā)者對(duì)部分地區(qū)的266個(gè)物種數(shù)據(jù)的測(cè)試，且以模擬其現(xiàn)實(shí)分布為目的[13]，所以在模擬物種未來(lái)潛在適生區(qū)時(shí)，繼續(xù)使用默認(rèn)參數(shù)會(huì)使MaxEnt模型產(chǎn)生過(guò)擬合的問(wèn)題，使預(yù)測(cè)結(jié)果不可靠[14]，需調(diào)用ENMeval數(shù)據(jù)包[15]對(duì)MaxEnt模型參數(shù)中的調(diào)控倍頻（regularization multiplier，RM）和特征組合（feature combination，F(xiàn)C）進(jìn)行調(diào)整，分析各種參數(shù)條件下模型的復(fù)雜度，選取最低復(fù)雜度的模型參數(shù)。特征類型分為：線型（linear）-L，二次型（quadratic）-Q，鉸鏈型（hinge）-H，乘積型（product）-P和閾值型（threshold）-T。RM的調(diào)整區(qū)間為0.5～4.0，每次增加0.5；由于樣本量較少，F(xiàn)C選取2種組合方式，分別為L(zhǎng)、LQ（同時(shí)使用線型與二次型）。利用ENMeval所輸出的各類數(shù)據(jù)選擇最適生的RM與FC。

2.1.3 模型評(píng)價(jià) 運(yùn)用赤池信息量準(zhǔn)則（AIC）評(píng)價(jià)模型的復(fù)雜度和擬合優(yōu)度，AIC越小，模型越好，通常選擇AIC最小的模型；運(yùn)用AUC（ROC曲線下方的面積大?。UCtest（訓(xùn)練AUC）、AUCdiff（測(cè)試AUC值之差）、orMTP（最小訓(xùn)練集遺漏率）和or10（10%訓(xùn)練遺漏率）[16-18]評(píng)估不同參數(shù)組合的擬合度，綜合以上數(shù)據(jù)選擇出最適合怒江金絲猴的參數(shù)組合。利用最優(yōu)參數(shù)組合在MaxEnt模型中進(jìn)行建模和預(yù)測(cè)適生區(qū)分布，AUC評(píng)價(jià)模型的性能好壞時(shí)取值范圍為0.5～1，AUC值越大表示模型性能越好，當(dāng)AUC值≥0.9時(shí)預(yù)測(cè)結(jié)果精度高且較為可靠。

2.1.4 適生區(qū)劃分 根據(jù)MaxEnt模型的結(jié)果，按照0～1的取值范圍輸出適生程度，0為最不適生，1為最適生。參考聯(lián)合國(guó)政府間氣候變化專門(mén)委員會(huì)第5次評(píng)估報(bào)告中的描述[19]，將適生區(qū)劃分為4個(gè)等級(jí)，分別為不適生區(qū)（0～0.3）、低適生區(qū)（0.3～0.5）、中適生區(qū)（0.5～0.7）和高適生區(qū)（0.7～1.0）。

2.2 物種分布數(shù)據(jù)物種分布數(shù)據(jù)來(lái)源于兩方面，一是全球生物多樣性信息網(wǎng)絡(luò)（https：//www.gbif.org/）、中國(guó)數(shù)字植物標(biāo)本館平臺(tái)（http：//www.cvh.ac.cn/）網(wǎng)站；二是相關(guān)文獻(xiàn)，根據(jù)怒江金絲猴相應(yīng)的出沒(méi)地點(diǎn)信息，借助 GeoNames（http：//www.geonames.org/）查詢其經(jīng)緯度數(shù)據(jù)。剔除冗余數(shù)據(jù)點(diǎn)、坐標(biāo)缺失的點(diǎn)及坐標(biāo)無(wú)意義的點(diǎn)，獲得44個(gè)怒江金絲猴分布點(diǎn)。

2.3 環(huán)境變量獲取及分析現(xiàn)代及未來(lái)氣候變量數(shù)據(jù)來(lái)源于WorldClim（https：//www.worldclim.org/），包含19個(gè)氣候變量，最冷季度和最干季度等氣候變量來(lái)源于每月的溫度和降雨量值，目的是生成更具生物意義的變量，由于不同地區(qū)出現(xiàn)極端環(huán)境的月份不同，所以在生成變量時(shí)不固定月份而是以最冷季度和最干季度等來(lái)表示?？臻g分辨率為1 km×1 km，時(shí)間分辨率分別為現(xiàn)代（1970—2000年），未來(lái)（a：2021—2040年；b：2041—2060年；c：2061—2080年；d：2081—2100年）。未來(lái)氣候模式基于第六次耦合模式比較計(jì)劃（CMIP6），使用新全球氣候模式BCC-CSM 2-MR和中等未來(lái)共享社會(huì)經(jīng)濟(jì)路徑（SSP245）數(shù)據(jù)。將19個(gè)氣候變量進(jìn)行聚類和皮爾遜相關(guān)性分析，在此基礎(chǔ)上篩選生物氣候變量。聚類分析將生物氣候變量分為七組（圖1）：（1）最冷月份最低溫度、最干季度平均溫度、最冷季度平均溫度；（2）年降水量、最暖季度降水量、最濕月份降水量、最濕季度降水量；（3）最熱月份最高溫度、最濕季度平均溫度、最暖季度平均溫度；（4）溫度變化系數(shù)、年溫差；（5）等溫性、晝夜溫差月均值、降水量變化方差；（6）年均溫；（7）最冷季度降水量、最干月份降水量、最干季度降水量。在相關(guān)性分析中，各組內(nèi)變量間的相關(guān)性均大于0.9，與聚類分析的結(jié)果相符。組內(nèi)挑選與其他變量最不相關(guān)的一個(gè)變量（表1），保留最重要的生態(tài)變量，單獨(dú)計(jì)算各變量與其他變量的相關(guān)性，取其絕對(duì)值的和作為單變量總相關(guān)性，組內(nèi)單變量總相關(guān)性最低的變量即為與其他變量最不相關(guān)的變量。篩選出7個(gè)變量（年均溫、等溫性、溫度變化系數(shù)、最濕季度平均溫度、最冷季度平均溫度、最暖季度降水量、最冷季度降水量），作為預(yù)測(cè)怒江金絲猴適生區(qū)的氣候變量。

圖1 不同生物氣候變量的聚類分析結(jié)果

表1 相關(guān)性分析結(jié)果

在氣候變量的基礎(chǔ)上加入影響怒江金絲猴生長(zhǎng)發(fā)育的其他環(huán)境變量，包括：高程（DEM）、土地利用（LUCC）、歸一化植被指數(shù)（NDVI）。高程數(shù)據(jù)來(lái)源于SRTM，土地利用數(shù)據(jù)與歸一化植被指數(shù)數(shù)據(jù)來(lái)源于中科院地理科學(xué)與資源研究所（https：//www.resdc.cn/Default.aspx）。通過(guò)ArcGIS剪裁與轉(zhuǎn)換得到與環(huán)境數(shù)據(jù)相同的空間分辨率。最終選擇10個(gè)變量（年均溫、等溫性、溫度變化系數(shù)、最濕季度平均溫度、最冷季度平均溫度、最暖季度降水量、最冷季度降水量、高程、土地利用、歸一化植被指數(shù)），作為預(yù)測(cè)怒江金絲猴適生區(qū)的環(huán)境變量（表2）。

表2 環(huán)境變量及其參數(shù)

2.4 模型參數(shù)選擇在小樣本情況下AIC轉(zhuǎn)變成AICc，隨著樣本量增加，AICc收斂成AIC。所以AICc可以應(yīng)用于任何樣本大小的情況下。得到各個(gè)模型的AICc值后，每個(gè)AICc值與最小的AICc值相減，得到ΔAICc。通過(guò)ENMeval計(jì)算 16種不同參數(shù)設(shè)置組合的 ΔAICc值（表3）。在 16種參數(shù)組合中，ΔAICc＜2的參數(shù)組合共有2組，分別為：RM=3，RM=3.5，F(xiàn)C均選取 LQ。它們的 ΔAICc分別為1.8837和0.0520，2組參數(shù)組合的ΔAICc都在接受范圍內(nèi)，其復(fù)雜度和擬合度均較低，具有較高的可靠性。進(jìn)一步對(duì) ENMeval生成的其他評(píng)估指標(biāo)進(jìn)行比較，當(dāng) RM=3，F(xiàn)C選取 LQ時(shí)，其 AUCdiff、orMTP和or10值均為最低，表明RM=3，F(xiàn)C選取LQ時(shí)的擬合度更好，并且其AUCtest的值更高，表明其能更好地區(qū)分測(cè)試分布點(diǎn)和背景點(diǎn)，因此，以RM=3，F(xiàn)C選取LQ作為最終的設(shè)置參數(shù)。在該設(shè)置參數(shù)條件下，MaxEnt模型模擬怒江金絲猴適生區(qū)的平均AUC值為0.9995，表明模型的預(yù)測(cè)效果良好。

表3 不同參數(shù)設(shè)置下模型評(píng)估結(jié)果

2.5 主導(dǎo)環(huán)境變量分析Maxent模型在訓(xùn)練時(shí)會(huì)記錄變量為適應(yīng)模型做出的貢獻(xiàn)，最終將變量貢獻(xiàn)量轉(zhuǎn)換為百分?jǐn)?shù)輸出為貢獻(xiàn)率，置換重要性是由該變量在訓(xùn)練過(guò)程中，訓(xùn)練點(diǎn)的隨機(jī)置換和由此導(dǎo)致的AUC值的下降幅度來(lái)決定，AUC下降幅度越大，模型對(duì)該環(huán)境變量依賴越大。不同氣候變量對(duì)MaxEnt模型的貢獻(xiàn)率和置換重要性（表2）結(jié)果顯示，共有4個(gè)氣候變量對(duì)模型的貢獻(xiàn)率大于10%，分別為等溫性（11.95%）、溫度變化系數(shù)（28.60%）、最冷季度平均溫度（16.44%）和 NDVI（17.39%），其累計(jì)貢獻(xiàn)率達(dá)74.38%，其中最冷季度平均溫度和NDVI的重要性較高，最冷季度平均溫度重要性達(dá)到50.17%。為了減少氣候變量之間的相互影響，以進(jìn)一步解釋哪些變量在模型中是最重要的，對(duì)10個(gè)變量的重要性進(jìn)行刀切（Jackknife）法檢測(cè)（圖2）。當(dāng)使用刀切法時(shí)，會(huì)創(chuàng)建多種模型，包括單獨(dú)使用每個(gè)變量創(chuàng)建一個(gè)模型、依次排除每個(gè)變量并使用剩余變量創(chuàng)建模型和使用所有變量創(chuàng)建模型。結(jié)果顯示，單獨(dú)使用年均溫、最濕季度平均溫度和土地利用時(shí)，幾乎沒(méi)有任何增益，表明這些參數(shù)對(duì)預(yù)測(cè)怒江金絲猴適生區(qū)的幫助很??；當(dāng)單獨(dú)使用時(shí)，增益最高的環(huán)境變量是溫度變化系數(shù)，因此它擁有最重要的信息，這與貢獻(xiàn)率結(jié)果相符合。當(dāng)忽略單個(gè)變量時(shí)，降低增益最多的環(huán)境變量是最冷季度降水量，因此它擁有其他變量中不存在的最多信息。綜上所述，溫度變化系數(shù)、最冷季度平均溫度、最冷季度降水量和歸一化植被指數(shù)是影響怒江金絲猴分布的主導(dǎo)變量。

圖2 變量重要性的刀切法檢測(cè)結(jié)果

3 結(jié)果與討論

3.1 怒江金絲猴現(xiàn)代適生區(qū)模擬現(xiàn)代背景下，44個(gè)怒江金絲猴分布點(diǎn)中屬于高適生區(qū)、中適生區(qū)、低適生區(qū)和不適生區(qū)的點(diǎn)數(shù)分別為 17個(gè)，19個(gè)，7個(gè)，1個(gè)，相應(yīng)比例分別為38.64%、43.18%、15.91%和2.27%，模擬結(jié)果（圖3）中我國(guó)境內(nèi)現(xiàn)存點(diǎn)均在高適生范圍內(nèi)，表明了預(yù)測(cè)結(jié)果的準(zhǔn)確性和MaxEnt模型用于評(píng)價(jià)怒江金絲猴適生區(qū)的可行性。怒江金絲猴在中國(guó)的現(xiàn)代適生區(qū)總面積約為30 790 km2，其中高適生區(qū)面積為6 207 km2，占適生區(qū)總面積的20.16%；主要適生區(qū)分布于云南西部的三江并流區(qū)域。高適生區(qū)集中于怒江兩岸的高黎貢山與怒山，垂直高差較大，有著典型的山地垂直植被帶和多樣的動(dòng)植物資源，為動(dòng)植物的生存提供了自低向高的海拔條件，具有典型的亞熱帶氣候特征。

圖3 怒江金絲猴的分布點(diǎn)與模擬現(xiàn)代適生區(qū)分布

3.2 怒江金絲猴適生區(qū)預(yù)測(cè)現(xiàn)代至未來(lái)4個(gè)階段，怒江金絲猴總適生區(qū)均有增長(zhǎng)（圖4）。根據(jù)模型對(duì)未來(lái)適生區(qū)的預(yù)測(cè)（圖5），云南西部的適生區(qū)在2021—2040年面積最大并且是高適生區(qū)最多的時(shí)間段，與現(xiàn)代相比總適生區(qū)面積增加98.59%，高適生區(qū)面積增加64.91%，主要表現(xiàn)為瀾滄江東部適生區(qū)的擴(kuò)大，并且延伸至金沙江以東；2021—2080年，各類型適生區(qū)面積都呈下降趨勢(shì)，向高適生區(qū)收縮，更集中于怒江與瀾滄江之間的山區(qū)；2041—2060年的適生區(qū)面積與輪廓和2081—2100年相似，分別比現(xiàn)代增加了79.36%和79.01%。

圖4 怒江金絲猴現(xiàn)代以及未來(lái)各適生區(qū)面積

圖5 未來(lái)氣候條件下怒江金絲猴的適生區(qū)分布預(yù)測(cè)

3.3 討論本研究中，綜合變量貢獻(xiàn)率、重要性以及刀切法檢測(cè)，發(fā)現(xiàn)溫度變化系數(shù)、最冷季度平均溫度、最冷季度降水量和歸一化植被指數(shù)是影響怒江金絲猴分布的主導(dǎo)變量。怒江金絲猴為典型的森林樹(shù)棲動(dòng)物，生活在1700 m至3100 m的原始森林中，其分布受到溫度、降水及棲息地植被情況的多重影響[1]，目前已知存在點(diǎn)高黎貢山區(qū)域海拔高差較大，且存在明顯的山地植被垂直分布帶[20]，而生境中的植物對(duì)怒江金絲猴的覓食行為產(chǎn)生影響，冬季的低溫及降雪令食物變得匱乏，使其向低海拔地區(qū)移動(dòng)生活，上述分析表明與低溫期相關(guān)的變量是決定怒江金絲猴地理空間分布的重要環(huán)境因素。

相較于已知存在區(qū)怒江以西的高黎貢山，模擬現(xiàn)代高適生區(qū)多出了怒江和金沙江之間的怒山山脈。實(shí)際情況下在這兩個(gè)區(qū)域并未觀測(cè)到怒江金絲猴的存在，在生境相似的情況下模擬結(jié)果會(huì)偏向于存在可能性高[21]，怒江以東的區(qū)域由于怒江、瀾滄江的隔離作用，即使西南山區(qū)生境相似，也并沒(méi)有怒江金絲猴生活在此區(qū)域，同為仰鼻猴屬的滇金絲猴分布于瀾滄江與金沙江之間云嶺山脈兩側(cè)，主要活動(dòng)在2500 m至4500 m高度的云杉、冷杉林中，兩個(gè)物種之間沒(méi)有發(fā)現(xiàn)交流跡象。在僅考慮下墊面情況下，模擬結(jié)果顯示在臺(tái)灣地區(qū)也存在怒江金絲猴的高適生區(qū)，但由于地理上的隔離怒江金絲猴不可能出現(xiàn)在臺(tái)灣，陳之瑞等[22]總結(jié)了中國(guó)西南地區(qū)與臺(tái)灣植物的間斷分布現(xiàn)象，表示兩者植物區(qū)系上有極大的相似性，臺(tái)灣的中央山脈最高主峰近4000 m，從植被和海拔氣候都與高黎貢山有十分相近的條件，這是模擬在臺(tái)灣也有高適生區(qū)分布的客觀原因，在進(jìn)行珍稀物種保護(hù)時(shí)可考慮讓臺(tái)灣地區(qū)成為異地遷移保護(hù)的備選區(qū)。

SSP245為綜合考慮了SSP-RCP的未來(lái)氣候情景，其中共享社會(huì)經(jīng)濟(jì)路徑為 “中間路線”的情景（SSP2），即在整個(gè)21世紀(jì)都將延續(xù)歷史發(fā)展的情景，典型濃度排放路徑為溫室氣體中等排放（RCP4.5）。在此情景下，未來(lái)中國(guó)西南部氣溫有顯著的增長(zhǎng)趨勢(shì)，降水在不同時(shí)段內(nèi)保持平穩(wěn)變化[23]。受到全球變暖的影響，對(duì)低溫敏感的怒江金絲猴適生區(qū)可能會(huì)向外擴(kuò)張，由瀾滄江向東擴(kuò)大，在近現(xiàn)代經(jīng)歷一個(gè)高速擴(kuò)張的時(shí)期，適生區(qū)達(dá)到最大，接著下降到較為穩(wěn)定的2041—2060年和2081—2100年時(shí)段，該時(shí)期相較于現(xiàn)代南北方向上縮減，東西方向上擴(kuò)大。綜合未來(lái)氣候情景下的適生區(qū)預(yù)測(cè)，怒江金絲猴可以長(zhǎng)久的生活在高黎貢山區(qū)域，并且適生區(qū)有所擴(kuò)大，可以在怒江以東以尋找其蹤跡，為該種金絲猴的種群保育工作和保護(hù)區(qū)劃分提供了依據(jù)。

4 結(jié)論

本研究利用MaxEnt優(yōu)化模型模擬了怒江金絲猴現(xiàn)代以及未來(lái)氣候情景下4個(gè)時(shí)段的適生區(qū)。模型的AUC值為0.9995，模型表現(xiàn)良好。結(jié)果表明，影響怒江金絲猴生存的主導(dǎo)變量為溫度變化系數(shù)、最冷季度平均溫度、最冷季度降水量和歸一化植被指數(shù)。怒江金絲猴在當(dāng)前氣候條件下適生區(qū)總面積為30 790 km2，高適生區(qū)面積為6 207 km2，主要集中在云南的高黎貢山區(qū)域，怒江與瀾滄江之間的怒山也是最適生怒江金絲猴生存的地方。未來(lái)?xiàng)l件下其適生區(qū)增加，2021—2040年為增加最多的時(shí)間段，相較于現(xiàn)代增加了總適生區(qū)98.59%，西南山區(qū)部分向東擴(kuò)張延伸至金沙江。本研究為怒江金絲猴這一極危物種的種群研究和保護(hù)區(qū)建立提供了科學(xué)依據(jù)。