基于最大熵生態(tài)位模型的中華穿山甲潛在適宜生境預(yù)測(cè)

塔 旗,李言闊,*,范文青,單繼紅,涂曉斌,應(yīng) 欽,繆瀘君,邵瑞清,申 錦

1 江西師范大學(xué)生命科學(xué)學(xué)院, 南昌 330022 2 江西省野生動(dòng)植物保護(hù)管理局, 南昌 330046 3 江西農(nóng)業(yè)大學(xué), 南昌 330045 4 江西省林業(yè)科學(xué)院, 南昌 330032

穿山甲是哺乳動(dòng)物中一個(gè)較小的類群,隸屬于哺乳綱(Mammalia) 、鱗甲目(Pholidota)、鯪鯉科 (Manidae)、穿山甲屬(Manis)[1]。目前,全世界現(xiàn)有8種穿山甲：中華穿山甲(Manispentadactyla)、印度穿山甲(M.crassicaudata)、馬來穿山甲(M.javanica)、大穿山甲(M.gigantea)、樹穿山甲(M.tricuspis)、南非穿山甲(M.temminckii)、長(zhǎng)尾穿山甲(M.teyradactyla)和巴拉望穿山甲(M.culionensis)。我國(guó)分布的穿山甲主要是中華穿山甲[2-4]。該物種2014年被IUCN紅色名錄(Red List)列為“極?！?CR)[5]；2017年被列入CITES附錄I(www.cites.org)；2020年6月5日,我國(guó)將其從國(guó)家二級(jí)保護(hù)動(dòng)物提升為國(guó)家一級(jí)保護(hù)動(dòng)物。在過去幾十年間,濫捕濫獵和非法貿(mào)易使中華穿山甲種群數(shù)量急速下降；據(jù)估計(jì),其國(guó)內(nèi)種群數(shù)量已減少90%以上[6- 7]。盡管中華穿山甲的瀕危程度不斷提高,但近年來我國(guó)一直沒有開展過全國(guó)性的野生中華穿山甲專項(xiàng)調(diào)查,也未見相關(guān)報(bào)道,其種群現(xiàn)狀尚不明顯,分布位置尚不清晰,在很大程度上限制了對(duì)該物種的有針對(duì)性的保護(hù)管理。因此,基于已知中華穿山甲的分布位點(diǎn)預(yù)測(cè)其潛在適宜區(qū)對(duì)于其種群保護(hù)和棲息地管理具有十分重要的意義。

掌握瀕危物種的種群狀況和生態(tài)需求是開展物種保護(hù)的基礎(chǔ)[8]。掌握目標(biāo)物種的空間分布和生境需求是瀕危物種保護(hù)的一項(xiàng)重要工作,但是有限的野外數(shù)據(jù)有時(shí)會(huì)限制對(duì)瀕危物種分布區(qū)的掌握。生態(tài)位模型通過目標(biāo)物種的地理分布數(shù)據(jù)和物種所在地的環(huán)境數(shù)據(jù),在物種已知的地理分布區(qū)域內(nèi)進(jìn)行測(cè)量或估算,進(jìn)而預(yù)測(cè)該區(qū)域內(nèi)哪些區(qū)域滿足該物種生態(tài)位的要求,從而成為該物種潛在分布區(qū)的一部分[9]。主流的生態(tài)位模型有GARP、ENFA、DOMAIN、BIOCLIM和MaxEnt。其中,最大熵模型(the maximum entropy model software,MaxEnt)預(yù)測(cè)的結(jié)果精確度更高[10]；即使在目標(biāo)物種分布點(diǎn)比較少的情況下,MaxEnt模型也可以獲得較準(zhǔn)確的預(yù)測(cè)結(jié)果[11-14]。自Phillips提出該模型以來,MaxEnt在外來物種入侵[15- 17]、動(dòng)植物保護(hù)[18-20]和全球氣候變化對(duì)物種分布的影響等方面得到廣泛應(yīng)用[21-23]。

在此背景下,本研究對(duì)近年來我國(guó)各省市有關(guān)中華穿山甲的救護(hù)記錄、救助新聞進(jìn)行了整理,收集了67個(gè)中華穿山甲分布數(shù)據(jù),選取地形地貌、地表類型、土地覆被、氣候變量、生物學(xué)變量和人類影響變量等25個(gè)全球尺度的環(huán)境變量,結(jié)合地理信息系統(tǒng)(GIS)、統(tǒng)計(jì)分析軟件SPSS和最大熵模型等軟件預(yù)測(cè)了中華穿山甲在我國(guó)的潛在適宜分布區(qū),以期能夠?yàn)槿珖?guó)中華穿山甲種群專項(xiàng)調(diào)查和種群保護(hù)管理提供科技支撐。

1 材料與方法

1.1 數(shù)據(jù)收集與處理

利用MaxEnt模型預(yù)測(cè)物種適宜分布區(qū)需要物種的分布數(shù)據(jù)和環(huán)境數(shù)據(jù)。中華穿山甲的分布數(shù)據(jù)主要來源于互聯(lián)網(wǎng)信息與各地的救助記錄。我們通過百度搜索引擎查找關(guān)鍵詞“中華穿山甲+省份”搜尋近十年來我國(guó)關(guān)于中華穿山甲的救護(hù)新聞,篩選出有關(guān)中華穿山甲的救護(hù)報(bào)道,提取位置信息,借助百度地圖輸入分布點(diǎn)地理位置信息,確定其十進(jìn)制格式的經(jīng)緯度。將上述途徑所得到的分布點(diǎn)進(jìn)行仔細(xì)篩選、整理,獲取分布點(diǎn)共67個(gè)(表1)。

表1 中華穿山甲地理分布位點(diǎn)

圖1 中華穿山甲救護(hù)記錄的空間分布圖 Fig.1 The spatial distribution map of rescued Chinese pangolin in China

使用ArcGIS軟件鄰域分析中的緩沖區(qū)模塊對(duì)分布數(shù)據(jù)點(diǎn)進(jìn)行剔除,保證每個(gè)柵格僅有一個(gè)分布點(diǎn),避免人為因素導(dǎo)致分布點(diǎn)過多而出現(xiàn)過度擬合現(xiàn)象。本研究所用環(huán)境數(shù)據(jù)空間分辨率為2.5 arc-minutes(約4.5 km),設(shè)置緩沖直徑為3 km,當(dāng)兩個(gè)分布點(diǎn)之間的距離小于3 km時(shí),緩沖區(qū)面積會(huì)重疊,此時(shí)只保留其中一個(gè)分布點(diǎn)[24],最后保留有效分布點(diǎn)65個(gè)(圖1)。分布數(shù)據(jù)按species、longitude和latitude格式依次輸入Excel單元格,并最終保存為MaxEnt模型支持的CSV格式文件用于后續(xù)分析。

本研究選取的環(huán)境數(shù)據(jù)包括1.地形地貌變量：海拔、坡度和坡向；2.地表類型變量：土地覆被類型；3.氣候變量：19個(gè)氣溫和降水因子,如年平均氣溫、最冷月最低溫度和最冷季度降水量等；4.生物學(xué)變量：歸一化植被指數(shù)(NDVI)；5.人類影響變量：人口密度。這些數(shù)據(jù)來源及生物學(xué)意義詳見表2。使用ArcGIS軟件將上述25個(gè)全球尺度的環(huán)境變量掩膜提取出中國(guó)區(qū)域的環(huán)境數(shù)據(jù),并重采樣統(tǒng)一精度2.5 arc-minutes,通過柵格轉(zhuǎn)換工具將柵格數(shù)據(jù)轉(zhuǎn)換為*.ASCII格式數(shù)據(jù)以進(jìn)行MaxEnt模型運(yùn)算。

19個(gè)氣候變量是以溫度和降水量為基礎(chǔ)數(shù)據(jù)計(jì)算得出,因此它們之間不可避免地存在自相關(guān)和多重共線性等問題,在模型運(yùn)行時(shí)帶入冗余信息,導(dǎo)致輸出的AUC值偏高[24]。為此,我們對(duì)所選的25個(gè)環(huán)境變量進(jìn)行了主成分分析,以剔除貢獻(xiàn)率低的環(huán)境變量,僅保留貢獻(xiàn)率大于1.0%的環(huán)境變量[25]。同時(shí),使用SPSS軟件對(duì)所有環(huán)境變量進(jìn)行Spearman相關(guān)性分析(附錄),如果兩個(gè)環(huán)境變量之間的相關(guān)系數(shù)絕對(duì)值大于0.8,說明兩個(gè)變量間存在顯著相關(guān)性,剔除兩者中對(duì)物種影響較小的環(huán)境變量[26-28]。最終,我們選擇了6個(gè)環(huán)境變量對(duì)中華穿山甲進(jìn)行適宜分布區(qū)預(yù)測(cè),包括最冷季度降水量(Bio_19)、人口密度(PD)、年降水量(Bio_12)、坡度(Slope)、坡面(Aspect)、海拔(Alt)。

1.2 模型運(yùn)算與評(píng)價(jià)

將處理好的分布數(shù)據(jù)與環(huán)境數(shù)據(jù)導(dǎo)入MaxEnt模型軟件,創(chuàng)建響應(yīng)曲線并設(shè)置75%的分布點(diǎn)作為訓(xùn)練集(training data),25%為隨機(jī)測(cè)試集(Random test percentage)；勾選“Random seed”選項(xiàng),重復(fù)次數(shù)(Replicates)設(shè)為10次[29-30]；高級(jí)設(shè)置中勾選“Write plot data”選項(xiàng),方便后期制圖與數(shù)據(jù)統(tǒng)計(jì),其他選項(xiàng)保持默認(rèn),進(jìn)行模型運(yùn)算。根據(jù)最大訓(xùn)練敏感度和特異度法[31- 33](Maximum Training Sensitivity Plus Specificity, MaxSS)確定閾值為0.166,然后通過重分類,將物種適宜區(qū)劃分為4個(gè)等級(jí)：適生值<0.166為非適宜區(qū)、0.166—0.332 為低適宜區(qū)、0.332—0.498為中適宜區(qū)、適生值>0.498為高適宜區(qū),接著計(jì)算4個(gè)適宜區(qū)在中國(guó)的分布面積。

采用刀切法(Do jackknife)來評(píng)價(jià)各變量對(duì)模型的重要性,使用受試者工作特征曲線(ROC, Receiver Operating Characteristic Curve)與橫坐標(biāo)圍成的面積即(Area Under Curve,AUC)來評(píng)價(jià)模型預(yù)測(cè)結(jié)果的準(zhǔn)確性。AUC值在0—1之間：0.5—0.6表現(xiàn)為失敗,0.6—0.7表現(xiàn)為較差,0.7—0.8表現(xiàn)為一般,0.8—0.9表現(xiàn)為好,0.9—1表現(xiàn)為極好；值越接近1表示模型預(yù)測(cè)準(zhǔn)確度越高。

2 結(jié)果與分析

2.1 主導(dǎo)環(huán)境因子貢獻(xiàn)率

通過MaxEnt模型分析了6種環(huán)境變量對(duì)中華穿山甲分布的貢獻(xiàn)率,結(jié)果表明氣候變量中的降水因子是限制中華穿山甲分布的主導(dǎo)因素,累計(jì)貢獻(xiàn)率達(dá)85.4%；其中,最冷季度降水量單因子貢獻(xiàn)率高達(dá)68.9%。人類影響因子和海拔因子的累計(jì)貢獻(xiàn)率達(dá)14.2%,坡面和坡度的貢獻(xiàn)率相對(duì)較小,兩者累計(jì)僅為0.4%(圖2)。

表2 中華穿山甲分布預(yù)測(cè)所使用的環(huán)境變量

圖2 影響中華穿山甲分布的環(huán)境變量貢獻(xiàn)率與累計(jì)貢獻(xiàn)率 Fig.2 The accumulated contribution of each environmental variable to the potential distribution of Chinese pangolin in China

2.2 刀切法分析單一環(huán)境變量的重要程度

我們用刀切法進(jìn)一步分析了單一環(huán)境變量對(duì)中華穿山甲分布結(jié)果的影響率(圖3)。藍(lán)色表示當(dāng)僅用此環(huán)境變量用于MaxEnt模型建模時(shí),該變量對(duì)模型的增益效果；綠色表示除此變量,其他變量對(duì)模型的總體增益效果；紅色表示所有環(huán)境變量應(yīng)用模型所產(chǎn)生的增益效果。結(jié)果表明最冷季度降水量和年降水量對(duì)中華穿山甲的分布起主導(dǎo)作用,訓(xùn)練增益分別為1.85和1.83；人口密度和海拔是影響中華穿山甲分布的次要因素,訓(xùn)練增益分別為0.86和0.71；影響最小的是坡面和坡度,訓(xùn)練增益不超過0.10。這與相對(duì)貢獻(xiàn)率分析的結(jié)果是一致的,表明所選的環(huán)境變量能夠較好地反映中華穿山甲真實(shí)的生活環(huán)境[25]。

圖3 刀切法分析主要環(huán)境變量對(duì)中華穿山甲在中國(guó)的分布影響程度Fig.3 Analysis of the impact of major environmental variables on the distribution of Chinese pangolin in China by Jackknife test

圖4 中華穿山甲對(duì)4個(gè)主要環(huán)境因子的響應(yīng)曲線Fig.4 Response of the Chinese pangolin to 4 major environmental factors Bio_19: 最冷季度降水量Precipitation of coldest quarter； PD: 人口密度 Population density； Bio_12: 年降水量Annual precipitation

2.3 主要環(huán)境變量的響應(yīng)分析

中華穿山甲的分布與主要的環(huán)境變量之間具有較強(qiáng)的響應(yīng)關(guān)系(圖4),當(dāng)存在概率大于0.498時(shí),可以認(rèn)為此環(huán)境變量值最適合中華穿山甲生存。其中,最冷季度降水量、年降水量、人口密度、海拔四個(gè)主要環(huán)境變量對(duì)中華穿山甲分布的最適宜范圍分別為：最冷季度降水量141.22—439.46 mm、年降水量1471.67—2386.56 mm、人口密度390—28894 人/km2、海拔<316.98 m。

最冷季度降水量小于29.80 mm時(shí),中華穿山甲的存在概率幾乎為零；小于93.85 mm時(shí),中華穿山甲的存在概率小于0.166,為非適宜區(qū)；93.85—120.54 mm時(shí),存在概率大于0.166且小于0.332,為低適宜區(qū)；120.54—141.22 mm或大于439.46 mm時(shí),存在概率大于0.332且小于0.498,為中適宜區(qū)；141.22—439.46 mm時(shí),存在概率大于0.498,為高適宜區(qū)。

年降水量小于1218.53 mm時(shí),為非適宜區(qū)；1218.53—1365.09 mm或大于3234.83 mm時(shí),為低適宜區(qū)；1365.09—1471.67 mm或2386.56—3234.83 mm時(shí),為中適宜區(qū)；1471.67—2386.56 mm時(shí),為高適宜區(qū),且在1560.50—1675.97 mm之間時(shí),中華穿山甲的存在概率達(dá)到峰值。

圖5 中華穿山甲在中國(guó)的適宜區(qū)分布Fig.5 Potential distribution of the Chinese pangolin in China

人口密度在30—6795 人/km2時(shí),存在概率隨人口密度的增加而增加,在6795—36742 人/km2時(shí),存在概率隨人口密度增加而減少。低適宜區(qū)范圍：0—30 人/km2；中適宜區(qū)范圍：30—390 人/km2或28894—32457 人/km2；高適宜區(qū)范圍：390—28894 人/km2。

海拔范圍大于3351.13 m時(shí),中華穿山甲存在概率為零。隨著海拔的降低,中華穿山甲的分布概率增加。海拔范圍在602.92—1166.86 m時(shí),為低適宜區(qū)；在324.92—563.20 m,為中適宜區(qū)；海拔<316.98 m,為高適宜區(qū)。

2.4 中華穿山甲在中國(guó)的適宜區(qū)預(yù)測(cè)

MaxEnt模型運(yùn)行十次后,篩選AUC值最高的模型作為最優(yōu)模型。按照重分類標(biāo)準(zhǔn)將中華穿山甲適宜性指數(shù)重分類成四類(圖5)。結(jié)果表明,中華穿山甲適宜分布區(qū)總面積約為74.27×104km2,占中國(guó)國(guó)土面積的7.73%,其中低、中、高適宜區(qū)面積分別為32.53×104km2、26.86×104km2和14.88×104km2,占總適宜區(qū)面積的43.80%、36.16%和20.04%。中華穿山甲在江西省和湖南省的高適宜區(qū)面積最大,分別為6.87×104km2和3.13×104km2,占該省面積的44.96%和16.14%；其次是廣東省和廣西省,面積分別為2.96×104km2和2.07×104km2,占該省面積的19.15%和9.91%,具體各適宜區(qū)在各省(市、自治區(qū))所占的面積及比例見表3。

表3 中華穿山甲適宜分布區(qū)主要覆蓋地區(qū)數(shù)據(jù)

2.5 ROC曲線和AUC值對(duì)模型預(yù)測(cè)的準(zhǔn)確性檢驗(yàn)

基于所有環(huán)境變量與主導(dǎo)環(huán)境變量構(gòu)建MaxEnt模型獲得ROC曲線和AUC值(圖6),獲得的訓(xùn)練集AUC值為別為0.982和0.961(誤差范圍±0.014),模型預(yù)測(cè)結(jié)果達(dá)到極好標(biāo)準(zhǔn),表明該模型對(duì)中華穿山甲在中國(guó)的潛在適宜區(qū)的預(yù)測(cè)結(jié)果非常好,可信度比較高。

圖6 MaxEnt模型預(yù)測(cè)結(jié)果的ROC曲線Fig.6 ROC curve of distribution of the Chinese pangolin predicted by MaxEnt model

3 討論

為了更有效地保護(hù)瀕危野生動(dòng)物,首先要明確的是該物種的分布情況,還要確定該物種的生境需求與所處環(huán)境之間的關(guān)系,從而為該物種提供最合理、最有效、最具體的保護(hù)方案[34]。本研究基于近年來我國(guó)中華穿山甲救護(hù)記錄,利用MaxEnt模型對(duì)其潛在適宜分布區(qū)進(jìn)行了預(yù)測(cè)。研究結(jié)果揭示了中華穿山甲在中國(guó)的適宜分布區(qū)主要在我國(guó)長(zhǎng)江以南地區(qū),主要集中在江西、廣東、湖南和浙江省,預(yù)測(cè)值與觀測(cè)值高度吻合,預(yù)測(cè)結(jié)果與中國(guó)獸類野外手冊(cè)[35]比較吻合,預(yù)測(cè)結(jié)果比較準(zhǔn)確,預(yù)測(cè)結(jié)果可為我國(guó)中華穿山甲的種群保護(hù)和棲息地管理提供有力的科技支撐。

分布數(shù)據(jù)和環(huán)境數(shù)據(jù)的過度擬合通常會(huì)導(dǎo)致模型輸出的AUC值偏高(圖6)。為了最大限度地避免數(shù)據(jù)點(diǎn)的偏差和模型的過度擬合,我們對(duì)已獲取的分布數(shù)據(jù)點(diǎn)進(jìn)行了篩選和剔除,保證每個(gè)柵格里只保留一個(gè)數(shù)據(jù)分布點(diǎn)。同時(shí),對(duì)所有的環(huán)境變量進(jìn)行了相關(guān)性分析來確保選取的環(huán)境數(shù)據(jù)與中華穿山甲的分布密切相關(guān),以此減少冗余信息對(duì)模型結(jié)果的影響,使模型輸出后的AUC值更加真實(shí)。MaxEnt模型輸出的AUC平均值0.961,表明該模型的預(yù)測(cè)結(jié)果具有較高的可信度。

我們也分析了影響中華穿山甲空間分布的關(guān)鍵環(huán)境因子。結(jié)果表明,氣候變量中的降水因子是限制中華穿山甲分布的主要因素,原因可能是降水因子與中華穿山甲的食物——白蟻的地理分布密切相關(guān)。1只中華穿山甲每年能吃700萬只螞蟻和白蟻,這些蟻類在我國(guó)的活動(dòng)分布主要在淮河以南的廣大地區(qū)(熱帶、亞熱帶等溫暖多濕地區(qū)),向北逐漸稀少,向南逐漸遞增。在這些地區(qū),白蟻的分布與雨量關(guān)系極為密切,尤其是江西及江西周圍省份溫度適宜、雨量充沛且白蟻種類繁多,這為中華穿山甲提供了很好的食物來源。同時(shí),白蟻的分布區(qū)也與本研究預(yù)測(cè)的中華穿山甲潛在適宜區(qū)高度吻合[36]。

物種與環(huán)境的關(guān)系是研究物種生態(tài)學(xué)需求和空間分布的重要方面。本研究認(rèn)為,雖然中華穿山甲的棲息地和食物來源會(huì)受到氣候變化的影響,但物種變化的過程是一個(gè)極其緩慢的過程[37- 38],物種的生態(tài)位壓縮主要還是由人類影響所驅(qū)動(dòng)。欒曉峰等研究發(fā)現(xiàn)短時(shí)間內(nèi)的氣候變化對(duì)于物種分布呈現(xiàn)中性的效果,或在某種情況下甚至可能是積極的,人類活動(dòng)干擾是影響中華穿山甲分布變化的最大因素[39]。國(guó)家林業(yè)局1998年前后開展的中國(guó)重點(diǎn)陸生野生動(dòng)物資源調(diào)查結(jié)果表明,我國(guó)野生穿山甲的數(shù)量約為64000只[3]。到2008年,張立等人認(rèn)為中國(guó)野生穿山甲的數(shù)量在25100—49450只之間,即10年內(nèi)其種群就減少了約1/2的數(shù)量[3]。雖然中華穿山甲的瀕危等級(jí)不斷提高,但近十年來我國(guó)一直沒有開展全國(guó)性的野生穿山甲專項(xiàng)調(diào)查,也未見相關(guān)報(bào)道,在很大程度上限制了對(duì)該物種有效保護(hù)。穿山甲是地球上最古老的物種之一,在地球上至少已延續(xù)了5000年。因其進(jìn)化程度低、對(duì)新的環(huán)境適應(yīng)能力差、繁殖力低下、人工難于馴養(yǎng)等原因,一旦遭到大量捕殺導(dǎo)致種群數(shù)量下降后就難以恢復(fù)[6]。如若不及時(shí)開展對(duì)中華穿山甲的生物學(xué)研究,潛在的滅絕危險(xiǎn)將進(jìn)一步增加。

MaxEnt模型在物種分布數(shù)據(jù)很少的情況下預(yù)測(cè)結(jié)果也比較準(zhǔn)確,但是不充分的分布信息也會(huì)導(dǎo)致模型預(yù)測(cè)產(chǎn)生偏差[40]。在避免過擬合的情況下,物種的分布數(shù)據(jù)越多越好。因此,基于更多野外活動(dòng)點(diǎn)的中華穿山甲分布區(qū)研究仍然有待深入。本研究所獲得的部分物種分布數(shù)據(jù)是通過百度地圖提取其救助地點(diǎn)所獲取的,少數(shù)新聞救助點(diǎn)的具體位置只能精確到鎮(zhèn)或鄉(xiāng)。因此,該救助地點(diǎn)的環(huán)境狀況無法完全反應(yīng)物種實(shí)際分布點(diǎn)的環(huán)境狀況,導(dǎo)致結(jié)果會(huì)有一定偏差。譬如，救護(hù)地點(diǎn)大多接近人類聚集區(qū)，周邊人口密度較高，可能會(huì)導(dǎo)致模型分析過高估計(jì)穿山甲適宜生境中的人口密度；本研究中MaxEnt模型預(yù)測(cè)穿山甲適宜生境的人口密度為390—28894 人/km2，雖然在一定程度上反映出穿山甲適宜生境距離人類活動(dòng)區(qū)較近，但可能過高地估計(jì)了穿山甲適宜生境中的人口密度。另外,雖然本研究選取的土地覆被類型和植被歸一化指數(shù)有助于精確模擬物種適宜生長(zhǎng)的微環(huán)境,但是兩者對(duì)于模型的貢獻(xiàn)度過小,使得我們無法從區(qū)域或微生境尺度分析環(huán)境因子對(duì)其分布的作用關(guān)系,只能從較大的空間尺度預(yù)測(cè)其潛在分布區(qū)。因此,多尺度和高精度上的中華穿山甲生境利用和空間分布適宜性預(yù)測(cè)有待進(jìn)一步開展。