肖 祺，劉 成，2，徐夢珍，周雄冬

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

1 研究背景

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

2.1.2 模型參數(shù) MaxEnt模型是一種基于機器學(xué)習(xí)的復(fù)雜模型，準(zhǔn)確率高但可轉(zhuǎn)移性差，大多數(shù)研究采用默認參數(shù)來構(gòu)建模型，默認參數(shù)的設(shè)置來源于模型開發(fā)者對部分地區(qū)的266個物種數(shù)據(jù)的測試，且以模擬其現(xiàn)實分布為目的[13]，所以在模擬物種未來潛在適生區(qū)時，繼續(xù)使用默認參數(shù)會使MaxEnt模型產(chǎn)生過擬合的問題，使預(yù)測結(jié)果不可靠[14]，需調(diào)用ENMeval數(shù)據(jù)包[15]對MaxEnt模型參數(shù)中的調(diào)控倍頻（regularization multiplier，RM）和特征組合（feature combination，F(xiàn)C）進行調(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、LQ（同時使用線型與二次型）。利用ENMeval所輸出的各類數(shù)據(jù)選擇最適生的RM與FC。

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

2.1.4 適生區(qū)劃分 根據(jù)MaxEnt模型的結(jié)果，按照0～1的取值范圍輸出適生程度，0為最不適生，1為最適生。參考聯(lián)合國政府間氣候變化專門委員會第5次評估報告中的描述[19]，將適生區(qū)劃分為4個等級，分別為不適生區(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ù)來源于兩方面，一是全球生物多樣性信息網(wǎng)絡(luò)（https：//www.gbif.org/）、中國數(shù)字植物標(biāo)本館平臺（http：//www.cvh.ac.cn/）網(wǎng)站；二是相關(guān)文獻，根據(jù)怒江金絲猴相應(yīng)的出沒地點信息，借助 GeoNames（http：//www.geonames.org/）查詢其經(jīng)緯度數(shù)據(jù)。剔除冗余數(shù)據(jù)點、坐標(biāo)缺失的點及坐標(biāo)無意義的點，獲得44個怒江金絲猴分布點。

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

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

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

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

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

2.4 模型參數(shù)選擇在小樣本情況下AIC轉(zhuǎn)變成AICc，隨著樣本量增加，AICc收斂成AIC。所以AICc可以應(yīng)用于任何樣本大小的情況下。得到各個模型的AICc值后，每個AICc值與最小的AICc值相減，得到ΔAICc。通過ENMeval計算 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ù)雜度和擬合度均較低，具有較高的可靠性。進一步對 ENMeval生成的其他評估指標(biāo)進行比較，當(dāng) RM=3，F(xiàn)C選取 LQ時，其 AUCdiff、orMTP和or10值均為最低，表明RM=3，F(xiàn)C選取LQ時的擬合度更好，并且其AUCtest的值更高，表明其能更好地區(qū)分測試分布點和背景點，因此，以RM=3，F(xiàn)C選取LQ作為最終的設(shè)置參數(shù)。在該設(shè)置參數(shù)條件下，MaxEnt模型模擬怒江金絲猴適生區(qū)的平均AUC值為0.9995，表明模型的預(yù)測效果良好。

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

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

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

3 結(jié)果與討論

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

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

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

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

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

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

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

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

4 結(jié)論

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