基于MaxEnt的無翼坡壘在三亞潛在適生區(qū)的研究

（中南林業(yè)科技大學，湖南 長沙 410004）

無翼坡壘Hopea exalata，又名鐵凌，是在1978年才發(fā)現(xiàn)的龍腦香科新種[1]，生于海拔400 m 左右的丘陵、坡地、山嶺的森林中[2]。它分布的地域范圍很小，僅在海南三亞東北部山區(qū)甘什嶺一帶見其蹤跡[3]，屬于珍稀瀕危種，國家Ⅱ級重點保護野生植物（國務院1999年8月4日批準）。目前，對無翼坡壘的研究報道較少且多為早期研究成果。楊小波等[4-5]通過最小面積取樣調(diào)查法得到了無翼坡壘的種群結(jié)構(gòu)呈現(xiàn)集群性的特征且隨著無翼坡壘個體越大，集群的程度不斷減小。隨后楊小波等[6]繼續(xù)從種群高度與種群密度的關(guān)系與種群生物量和種群稀疏規(guī)律解釋了無翼坡壘種群調(diào)節(jié)的規(guī)律。胡榮桂等[7]研究了無翼坡壘不同立木級下葉、樹干皮與根中的N、P、K 營養(yǎng)元素的含量變化。伍月花等[8]總結(jié)出無翼坡壘有葉尖干枯病、灰斑病、褐斑病、枯斑病4 種病害，并分析了4 種病害的癥狀及發(fā)生環(huán)境。羊金殿等[9]提出在無翼坡壘的遷地保護中適當?shù)恼陉幱欣跓o翼坡壘的生長發(fā)育。后人也有對無翼坡壘裸根壘幼苗移植技術(shù)[10]、無翼坡壘的生長與光合生理特性的差異[11]、無翼坡壘于土壤營養(yǎng)元素及化學計量學特征[12]方面的研究。但目前對無翼坡壘的潛在適生區(qū)的研究或影響無翼坡壘分布的環(huán)境條件的研究都暫無報道。

關(guān)于物種分布格局分析及潛在分布預測的模型眾多，研究表明在BIO_CLIM、DOMAIN、ENFA、GARP、Mahalanobis、MaxEnt 這常用的6 種生態(tài)位模型中，MaxEnt 的預測效果最好[13]。MaxEnt 是以最大熵理論為基礎(chǔ)，將物種分布點數(shù)據(jù)與對應的環(huán)境變量相結(jié)合，從符合條件的分布中選擇熵最大的分布作為最優(yōu)分布，建立物種最大可能預測模型。MaxEnt 模型被廣泛地應用于珍稀動物[14]、病蟲害[15]、入侵植物[16]甚至是病 菌[17]的潛在適生區(qū)的預測。而且MaxEnt 模型的準確度對樣本量并不太敏感，是目前針對可用分布數(shù)據(jù)較少情況下最有效的模型[18]。綜上所述，MaxEnt 模型在物種潛在分布預測方面的應用已十分成熟且不僅可以處理復雜的非線性關(guān)系，還能預測模型結(jié)果并檢驗，現(xiàn)根據(jù)無翼坡壘分布范圍小、可用分布點樣本少的特點，本研究選擇MaxEnt 模型對無翼坡壘的潛在適生區(qū)進行研究。研究得出的適生區(qū)可為對無翼坡壘的保護區(qū)邊界的劃定提供參考依據(jù)，分析得到的主要影響無翼坡壘分布的環(huán)境變量，可為無翼坡壘的遷地保護提供理論支持。

1 材料與方法

1.1 無翼坡壘分布點數(shù)據(jù)

無翼坡壘獨特的分布特點決定了獲取無翼坡壘分布點的難度較大，且分布點樣本點較少。一般常用的物種分布點數(shù)據(jù)大多來自各平臺物種標本館，在物種分布點數(shù)據(jù)量上占較大優(yōu)勢。而無翼坡壘的分布點數(shù)據(jù)在各平臺物種標本館皆不可獲得，根據(jù)無翼坡壘目前可考文獻得知無翼坡壘的分布主要集中分布于海南三亞甘什嶺一帶，因此構(gòu)建MaxEnt 模型所需的無翼坡壘分布點數(shù)據(jù)主要來源于海南甘什嶺省級自然保區(qū)科學考察報告（2011年）以及外業(yè)實地考察作為驗證（2018年，胸徑起測徑階為5 cm）。外業(yè)考察中于海南甘什嶺省級保護區(qū)實驗區(qū)及其南部范圍共設置兩條樣線，調(diào)查20 m×20 m 的樣地共41 塊。

調(diào)查過程中主要記錄的數(shù)據(jù)為樣地中無翼坡壘分布的經(jīng)緯度、高程、坡度和坡向，便于后續(xù)數(shù)據(jù)處理與驗證。將科考報告和外業(yè)樣方調(diào)查中的無翼坡壘分布點數(shù)據(jù)整理成經(jīng)緯度數(shù)據(jù)，剔除地理位置重復與異常數(shù)據(jù)后共得到無翼坡壘有效分布樣點26 個。無翼坡壘分布點分布見圖1。

1.2 環(huán)境數(shù)據(jù)

參考常見預測物種分布模型所需環(huán)境變量，并根據(jù)無翼坡壘的生長特性，現(xiàn)選取氣候數(shù)據(jù)、地形數(shù)據(jù)、徑流累積量數(shù)據(jù)、土壤數(shù)據(jù)以及太陽輻射數(shù)據(jù)（從高程數(shù)據(jù)計算獲得的太陽柵格圖）共5 套環(huán)境數(shù)據(jù)?？紤]到無翼坡壘胸徑年平均生長量為0.13～0.18 cm，而無翼坡壘在幼年時期的發(fā)病率較高存活率較低[8-9]，在隨著苗齡的增加，對環(huán)境的適應能力將不斷增強[2]，故在綜合考慮數(shù)據(jù)的可獲得性的情況下，本研究確定的氣候區(qū)間為1970—2000年。數(shù)據(jù)來源于Worldclim 氣候數(shù)據(jù)庫（http://worldclim.org/），數(shù)據(jù)庫內(nèi)包括19個氣候變量（BIO_1～BIO_19）的1970—2000年間共30年間的平均值，分辨率為30″×30″（約 1 km×1 km）。該套氣候變量是世界各地氣象站的記錄數(shù)據(jù)，通過空間差值生成全球氣候柵格數(shù)據(jù)，這被廣泛地應用于潛在適生區(qū)的研究。生物氣候變量信息見表1。土壤數(shù)據(jù)來源于“黑河計劃數(shù)據(jù)管理中心”(http://westdc.westgis.ac.cn)，采用的土壤分類系統(tǒng)主要為FAO-90，數(shù)據(jù)格式為投影成WGS84 坐標的grid 柵格數(shù)據(jù)。土壤變量中以T_開頭屬性字段表示上層土壤屬性（0～30 cm），若以S_開頭屬性字段表示下層土壤屬性（30～100 cm）。土壤變量信息見表2。

圖1 無翼坡壘分布點Fig.1 Distribution point of Hopea exalata

表1 生物氣候變量及其信息Table 1 Information of BIO_climatic variables

地形數(shù)據(jù)中高程數(shù)據(jù)來源于地理空間數(shù)據(jù)云GDEMDEM 30M 分辨率數(shù)字高程數(shù)據(jù)（http://www.gscloud.cn/search），數(shù)據(jù)格式為投影成WGS84 坐標的tiff柵格數(shù)據(jù)。其中坡度、坡向、匯流累積量和太陽輻射數(shù)據(jù)均可由ArcGIS 空間分析工具分析高程數(shù)據(jù)得到，得到的數(shù)據(jù)的坐標和分辨率保持與原文件相一致。

表2 土壤變量及其信息Table 2 Information of soil variables

1 .3 研究方法

本研究基于MaxEnt 模型，將無翼坡壘分布點和環(huán)境數(shù)據(jù)作為源數(shù)據(jù)預測無翼坡壘在甘什嶺一帶的潛在適生區(qū)。MaxEnt 軟件的運行需要環(huán)境變量具有相同的空間范圍、坐標和分辨率。現(xiàn)通過ArcGIS 將以上環(huán)境數(shù)據(jù)統(tǒng)一定義為WGS_1984 坐標，并利用掩膜提取工具統(tǒng)一環(huán)境數(shù)據(jù)的空間范圍和分辨率一致，分辨率設置統(tǒng)一為30 m。將預處理好之后的各環(huán)境柵格變量通過轉(zhuǎn)換工具轉(zhuǎn)換成asc 格式文件，無翼坡壘分布點數(shù)據(jù)轉(zhuǎn)換成經(jīng)緯度數(shù)據(jù)導入Excel 表，另存為csv 格式。

由于參與建模的環(huán)境變量的數(shù)據(jù)量較大，為排除環(huán)境變量之間的數(shù)據(jù)聯(lián)系同時使得參與建模的環(huán)境變量具有代表性，現(xiàn)對19 個氣候變量和32個土壤變量標準化后分別進行主成分分析，選取累計方差百分比達90%的主成分。氣候變量與土壤變量總方差解釋達90%以上數(shù)據(jù)分別見表3與表4。在氣候變量篩選出的3 個主成分中，選取各主成分中成分得分系數(shù)絕對值最大的氣候變量，得到BIO_2、BIO_8、BIO_9、BIO_13、BIO_14、BIO_18、BIO_19 共7 個氣候變量，代替原19 個氣候變量參與建模。在土壤變量中選取各主成分中成分得分系數(shù)絕對值最大的土壤變量，得到 s_oc、s_ref_bulk、t_bs、t_gravel、t_oc、t_ph_h2o、t_ref_bulk、t_sand 共8 個土壤變量代替原32 個土壤變量參與建模。氣候變量與土壤變量主成分分析得分系數(shù)矩陣見表3和表4。最后將MaxEnt 軟件生成的asc 格式柵格導入ArcGIS 分級顯示：適宜性指數(shù)范圍在（0，0.3]為不適宜區(qū)；適宜性指數(shù)范圍在（0.3，0.45]為低適宜區(qū)；適宜性指數(shù)范圍在（0.45，0.6]，為一般適宜區(qū)；適宜性指數(shù)范圍在（0.6，0.75]，為高適宜區(qū)；適宜性指數(shù)范圍在（0.75，1]，為核心適宜區(qū)。

表3 氣候變量總方差解釋Table 3 Interpretation of total variance of climate variables

表4 土壤變量總方差解釋Table 4 Interpretation of total variance of soil variables

表5 氣候變量成分系數(shù)矩陣Table 5 Component coefficient matrix of climate variables

表6 土壤變量成分系數(shù)矩陣Table 6 Component coefficient matrix of soil variables

本次研究以氣候數(shù)據(jù)、土壤數(shù)據(jù)、地形數(shù)據(jù)、徑流累積量和太陽輻射量之間相關(guān)性低為前提，因此只對氣候數(shù)據(jù)與土壤數(shù)據(jù)分別進行主成分分析。將主成分分析篩選得出的7 個氣候變量、8 個土壤變量、高程、坡度、坡向、徑流累積量和太陽輻射量共20 個變量輸入MaxEnt 模型進行初次建模。在篩選初次建模后對模型貢獻率低于1 的環(huán)境變量，得到了BIO__14.tif、高程、BIO__19.tif、BIO__02.tif、BIO__13.tif、t_gravel、BIO__18.tif 7 個環(huán)境變量。將這7 個環(huán)境變量作為二次建模環(huán)境變量導入MaxEnt 模型得到最終結(jié)果。

2 結(jié)果與分析

2.1 模型精度檢驗

MaxEnt 模型利用ROC 曲線檢驗模型精度，結(jié)果顯示訓練樣本AUC 值為0.991，測試樣本的AUC 值為0.992。AUC 值非常接近1，說明篩選后環(huán)境變量對無翼坡壘潛在適生區(qū)分布的關(guān)系大，該次模型預測精度高，本次模型運行取得很好的結(jié)果。無翼坡壘的ROC 曲線驗證見圖2。

圖2 MaxEnt 模型對預測結(jié)果的ROC 曲線驗證Fig.2 Verification of ROC curve of prediction results by MaxEnt model

2.2 無翼坡壘的潛在適生區(qū)

在對無翼坡壘在三亞的適生區(qū)進行預測后，得到無翼坡壘在三亞的潛在分布（圖3）。無翼坡壘的核心適宜區(qū)主要分布于甘什嶺森林公園核心區(qū)范圍內(nèi)，在大茅村北部與羅蓬村交界處也有少許分布，兩處核心適宜區(qū)的總面積約632.09 hm2；高適宜區(qū)的總面積約為825.12 hm2，其分布除了甘什嶺保護范圍，也分布于大茅村與三亞北部邊界三郎嶺林地處。預測結(jié)果中南田農(nóng)場愛國分廠的東北部和東風分廠東北部赤田水庫北面林地也有極少高適宜區(qū)分布；一般適宜區(qū)的總面積約為699.93 hm2，其主要緊密圍繞于高適宜區(qū)周圍分布，也有少部分分布于坡灣村中甘什嶺森林公園范圍內(nèi)；三亞市其他區(qū)縣及島嶼都多為適生指數(shù)較低的不適宜區(qū)與低適宜區(qū)。

發(fā)現(xiàn)預測結(jié)果與現(xiàn)有資料對無翼坡壘的分布描述相同，主要在甘什嶺一帶有分布，因此將甘什嶺一帶作為后續(xù)分析的研究區(qū)。研究區(qū)的劃定在保證連續(xù)性的同時，排除水域、村莊和交通運輸用地后，并使得研究區(qū)內(nèi)的地類基本為林地，最終確定無翼坡壘潛在適生區(qū)的研究區(qū)。研究區(qū)橫跨三亞市吉陽區(qū)與海棠區(qū)，貼近三亞市的西北面邊界，總面積為5 596.69 hm2，無翼坡壘在研究區(qū)范圍內(nèi)的潛在適生區(qū)分布見圖4。MaxEnt 模型的輸出結(jié)果統(tǒng)計如下：不適宜區(qū)面積為3 078.72 hm2，約占研究區(qū)面積的55%；低適宜區(qū)面積為 585.29 hm2，約占研究區(qū)面積的10%；一般適宜區(qū)面積為529.46 hm2，約占研究區(qū)面積的9%；高適宜區(qū)面積為711.13 hm2，約占研究區(qū)面積的14%；核心適宜區(qū)面積為632.09 hm2，約占研究區(qū)面積的11%，主要分布于甘什嶺自然保護區(qū)的核心區(qū)。

2.3 環(huán)境變量結(jié)果分析

在訓練算法的每次迭代中，將正則增益的增加量加到相應變量的貢獻中，若λ的絕對值的變化是負的，則從其減去，定義此為第一估計值即貢獻率。若對每個環(huán)境變量的值隨機置換，在置換后的數(shù)據(jù)基礎(chǔ)上重新評估模型，將訓練AUC 的下降結(jié)果數(shù)據(jù)歸一化為百分比定義為第二估計值即排列重要性。隔離環(huán)境變量單獨使用時增益最高的環(huán)境變量為BIO_14，因此，BIO_14 本身似乎具有最有用的信息。環(huán)境變量的刀切試驗結(jié)果見圖5。根據(jù)模型輸出的結(jié)果反饋可以直觀地看到，對模型的貢獻率從大到小依次是：BIO_14、高 程、BIO_19、BIO_2、BIO_13、t_gravel、BIO_18。其中BIO_14 貢獻率為41.4%，排列重要性為42.3%，對模型貢獻率與排列重要性均為第一；高程對模型的貢獻率為17.9%，排列重要性為34.8%；BIO_19 對模型的貢獻率為12.8%，排列重要性為7.7%。綜上所述BIO_14、高程和BIO_19三個環(huán)境變量對模型的貢獻率最大，貢獻率總和為72.1%，重要性總和為84.8%；t_gravel 對模型的貢獻率只有5.2%，但排列重要性排第三達9.5%。環(huán)境變量對模型貢獻率見表5。

圖3 無翼坡壘在三亞的潛在適生區(qū)分布Fig.3 Distribution of potential suitable growing area of Hopea exalata in sanya

圖4 無翼坡壘在甘什嶺的潛在適生區(qū)分布Fig.4 Distribution of potential suitable growing area of Hopea exalata in Ganshenling

MaxEnt 模型通過迭代法不斷修改環(huán)境變量的參數(shù)來計算不同環(huán)境變量對于模型的貢獻率，即在其他環(huán)境變量保持在其平均樣本值條件下，分別對BIO_14、高程、BIO_19、t_gravel 四個環(huán)境變量與預測概率變化的關(guān)系進行分析，并以響應曲線的形式進行表現(xiàn)。同時通過單因子響應曲線與變量在所有環(huán)境變量參與下的響應曲線進行對比，來確定環(huán)境變量的閾值。同一變量的兩個響應曲線的一致性越高，也進一步說明了該環(huán)境變量對模型的貢獻率高，是確定主導變量的依據(jù)之一。響應曲線結(jié)果顯示BIO_14 的閾值為 19～25 mm，在BIO_14 達到25 mm 左右時，適應性最高；高程的閾值為0～1 200 m，當高程在300 m 左右時，適宜性最高；BIO_19 閾值為 70～105 mm，當BIO_19 在90 左右時，適應性最高；t_gravel 的閾值為0%～30%，當t_gravel 為30% 時，適宜性最低。BIO_14、高程、BIO_19、t_gravel 的響應曲線見圖6。

表7 環(huán)境變量對模型貢獻率Table 7 Contribution rate of environmental variables to model

圖5 環(huán)境變量對模型貢獻率Fig.5 Contribution rate of environmental variables to model

圖6 主導環(huán)境變量響應曲線Fig.6 Response curve of dominant environmental variables

3 結(jié)論與討論

3.1 結(jié) 論

通過整理無翼坡壘物種分布數(shù)據(jù)、氣候變量、地形變量和土壤變量，用MaxEnt 模型軟件對無翼坡壘在海南甘什嶺一帶的潛在分布進行模擬預測，同時對主導無翼坡壘潛在分布的環(huán)境變量進行分析，主要結(jié)論如下：

在總共56 個環(huán)境變量中，BIO_14、高程、BIO_19 和t_gravel 對無翼坡壘的適生區(qū)有重要影響。其中BIO_14 對無翼坡壘分布的影響最大,貢獻率為41.4%，高程的貢獻率為17.9%，BIO_19的貢獻率為12.8%，t_gravel 的貢獻率為5.2%。環(huán)境變量的響應曲線表明：無翼坡壘在干燥和寒冷季對降水量的需求較低；上層土壤中較低的碎石比例利于無翼坡壘的生長，上層碎石比例達到15%時存在無翼坡壘的概率幾乎為0；無翼坡壘最適生長高程范圍為200～400 m。

根據(jù)無翼坡壘潛在適生區(qū)分布圖可得：無翼坡壘的核心適宜區(qū)和高適宜區(qū)主要分布于甘什嶺保護區(qū)的核心區(qū)范圍內(nèi)，還有部分核心適宜區(qū)分布于甘什嶺保護區(qū)外西北面；一般適宜區(qū)的分布除了分布于甘什嶺保護區(qū)范圍內(nèi)，還在保護區(qū)外南邊有零星分布，無翼坡壘的潛在適生區(qū)有往南擴張的趨勢，可能是無翼坡壘群落仍處于演替階段，尚未達到穩(wěn)定[19]。

3.2 討 論

本研究包含的環(huán)境變量為氣候數(shù)據(jù)、土壤數(shù)據(jù)、地形數(shù)據(jù)、徑流累積量和根據(jù)地形模擬的太陽輻射總量。其中氣候數(shù)據(jù)、土壤數(shù)據(jù)、地形數(shù)據(jù)為一般預測植物潛在分布的常用變量[20-23]，在此基礎(chǔ)上本研究加入了徑流累積量和根據(jù)地形模擬的太陽輻射兩組變量。在對環(huán)境變量的選擇時，由于甘什嶺設有省級自然保護區(qū)，人為活動痕跡較不明顯，因此并未引入人為活動的干擾作為環(huán)境變量。且無翼坡壘與其他樹種間的競爭關(guān)系較難數(shù)據(jù)化為柵格變量，故未考慮種間競爭對無翼坡壘適生的影響。在獲取無翼坡壘分布點方面，由于現(xiàn)有野外無翼坡壘的記錄僅在甘什嶺一帶，這使得分布點可能沒有完全代表物種的全部實際生境[24]，這對可預測的范圍也是極大的限制，過大的研究范圍可能會導致極低的預測準確性。

無翼坡壘在三亞的適生區(qū)的預測結(jié)果為甘什嶺一帶，范圍包括海南甘什嶺省級自然保護區(qū)與海南甘什嶺省級森林公園。但是仍有部分適生區(qū)存在于保護區(qū)與森林公園的范圍外，后續(xù)可在本研究的基礎(chǔ)上對保護區(qū)進行空缺分析，這將更好地對無翼坡壘的生境進行保護。研究結(jié)果顯示BIO_14、高程、BIO_19 和t_gravel 對無翼坡壘的適生區(qū)選擇有重要影響。無翼坡壘在最干季與最冷季對降水量的需求不高可能是由于甘什嶺一帶干濕季明顯，每年5—10月為雨季，11月至翌年4月為旱季，而最冷季與旱季持續(xù)時間基本一致且降水量都較低?？赡苁怯捎诟适矌X的季節(jié)性降水特點使得無翼坡壘的分布主要集中分布于那一帶。在對無翼坡壘的遷地保護的過程中，可為無翼坡壘提供含有較低碎石比例的表層土壤來保證較好的適生環(huán)境，同時需在最干季與最冷季控制對無翼坡壘的澆水量。