李慧，張志華，薦圣淇

（1.開封市水務(wù)開發(fā)建設(shè)有限公司，河南 開封 475000；2.鄭州大學(xué) 黃河實驗室，河南 鄭州 450001）

油松Pinus tabuliformis是我國主要的用材樹種之一，種植范圍廣泛，對于水土保持起著重要的作用。構(gòu)建油松環(huán)境空間分布和棲息地適應(yīng)性模型對于人工植被恢復(fù)和油松生態(tài)保護(hù)具有重要參考價值。張雷等[1]采用8個生態(tài)位模型分析了油松當(dāng)前和未來氣候條件下3 個時間段的潛在分布，但建模數(shù)據(jù)中僅包含7 個氣候變量，缺乏重要的地形數(shù)據(jù)。肖敏等[2]選取全國氣象站點插值的10 個氣候環(huán)境因子，利用最大熵（MaxEnt）模型模擬我國油松分布區(qū)，結(jié)果表明我國油松中高度適生區(qū)域主要是我國華北、西北和東北等地[3-6]。

物種分布模型是生態(tài)學(xué)研究中最重要的工具之一，通常用于研究入侵害蟲和植物物種的潛在威脅[7]。這些模型也是預(yù)測氣候變化對物種分布影響的主要工具[8]。雖然已經(jīng)建立了各種生態(tài)模型來預(yù)測分布，相關(guān)研究發(fā)現(xiàn)，在預(yù)測精度方面，MaxEnt 模型較其他模型具有明顯的優(yōu)勢，可以有效解決不完整的物種分布情況[9]。MaxEnt 模型已被廣泛應(yīng)用于瀕危物種適生區(qū)預(yù)測、特定物種氣候環(huán)境適生性預(yù)測和物種保護(hù)優(yōu)先級評估[10-13]。

本研究以伊洛河流域為研究區(qū)，以油松為研究對象，建立MaxEnt 模型，模擬伊洛河流域油松潛在空間分布格局劃分油松適生區(qū)等級，確定影響伊洛河流域油松分布的環(huán)境主導(dǎo)因子，為油松在伊洛河流域的引種與栽培提供科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

伊洛河流域是黃河中下游河段右岸重要支流，位于陜西地區(qū)東南部及河南省西北部 地理坐標(biāo)為109°43～ 113°10' E，33°39'～ 34°54' N，見圖1。伊洛河流域地處暖溫帶大陸性季風(fēng)氣候區(qū)，多年平均氣溫為13 ℃，降水集中在6—9 月，年均降水量為710～ 930 mm，年徑流量為13.99 億m3。伊洛河流域位于秦嶺構(gòu)造體系東段與新華夏系的復(fù)合部位，主要分布有花崗巖、玄武巖、第四系黃土、淤泥質(zhì)砂壤土等；地勢由西南向東北降低，分別為山地、丘陵和平原過渡帶，地帶性植被為暖溫帶落葉闊葉林。

圖1 伊洛河流域油松分布樣點位置Fig.1 Sampling sites in distribution area of P.tabuliformis in the Yiluo River Basin

1.2 數(shù)據(jù)來源

1.2.1 油松居群地理數(shù)據(jù) 油松地理分布信息數(shù)據(jù)由河南省林業(yè)局調(diào)查檔案和野外觀測獲取，通過去除無效和重復(fù)數(shù)據(jù)，共收集到油松分布樣本點109個（圖1）。當(dāng)收集區(qū)域比較密集產(chǎn)生空間聚類時，往往導(dǎo)致模型過度擬合，模型性能值下降。本文通過ArcGIS 10.5 以點距離大于5 km 的標(biāo)準(zhǔn)對采樣點進(jìn)行重采樣，刪除具有高度空間自動相關(guān)記錄的原始分布數(shù)據(jù)，并在采樣過程中盡量保留采樣點，最終獲取76 個有效分布點。

1.2.2 環(huán)境變量 從世界氣候數(shù)據(jù)庫（https://www.worldclim.org/）[14]中獲取氣候變量，該生物氣象數(shù)據(jù)集由1970—2000 年世界氣候2.1 版氣候數(shù)據(jù)插值而成；海拔（Altitude）、坡度（Slope）、坡向（Aspect）3 個地形變量通過30 m×30 m 分辨率的DEM 提取；土壤變量從1∶100 萬中華人民共和國土壤圖中獲?。ū?）。

表1 相關(guān)環(huán)境變量名稱及描述Tab.1 Description of environmental variables

為避免模型過度擬合，利用Spearman 和Pearson 相關(guān)檢驗對油松的19 個生物氣象數(shù)據(jù)進(jìn)行相關(guān)性分析。利用油松點隨機(jī)抽取的76 個交叉相關(guān)矩陣，估算出19 個氣候變量子集中的協(xié)方差。在此基礎(chǔ)上，去除相關(guān)系數(shù)大于0.70 的環(huán)境變量，以便于對環(huán)境信息及其與生境的關(guān)系進(jìn)行合理地解釋。另外，利用MaxEnt v3.4.1 的刀切法（Jackknife）直接判斷油松預(yù)測變量的重要性。最終篩選出10 個環(huán)境變量（表2）。

表2 篩選后的環(huán)境變量名稱及描述Tab.2 Description of environmental variables after screening

1.3 MaxEnt 模型構(gòu)建與評價方法

MaxEnt 模型基于生態(tài)位原理，可以只依據(jù)物種“存在點”的環(huán)境變量特征確定物種分布限制條件，通過識別最大熵分布的可能性對物種生境分布進(jìn)行預(yù)測。本文將76 個油松地理分布樣本點和篩選出相關(guān)性低的10 個環(huán)境特征變量導(dǎo)入MaxEnt v3.4.1，隨機(jī)選擇75%的油松分布點進(jìn)行建模，25%的油松分布點進(jìn)行模型驗證。本文基于樣本點隨機(jī)最大選取10 000 個背景點進(jìn)行模擬，在數(shù)據(jù)處理完畢后，其他參數(shù)保持默認(rèn)值。本文采用非閾值依賴的受試者工作特征曲線（receiver operating curve,ROC）下面積AUC 評價模型，AUC 越接近1，說明預(yù)測的結(jié)果越準(zhǔn)確[15-16]。

1.4 油松適生區(qū)等級劃分

將MaxEnt 預(yù)測結(jié)果在ArcGIS 10.5 中轉(zhuǎn)化為柵格數(shù)據(jù)，并利用其空間分析的重分類工具依據(jù)自然分割法（Jenks）和實際情況對適生區(qū)柵格數(shù)據(jù)圖層進(jìn)行合理的劃分，共劃分為以下4 個等級：0～ 0.2、? 0.2～ 0.4、? 0.4～ 0.7、?0.7，其分別代表非適生區(qū)、低適生區(qū)、中適生區(qū)和高適生區(qū)。

2 結(jié)果與分析

2.1 MaxEnt 模型效果

MaxEnt 模型的構(gòu)建采用油松75%的隨機(jī)地理分布樣點作為訓(xùn)練集，25%的分布樣點作為測試集。AUC在隨機(jī)預(yù)測下為0.500，經(jīng)過20 次重復(fù)建模，ROC 曲線下的AUC值為0.957，測試數(shù)據(jù)集的AUC為0.935，預(yù)測值均在0.90 以上，標(biāo)準(zhǔn)差為0.059，說明MaxEnt 模型具有良好的預(yù)測效果，見圖2。

圖2 MaxEnt 模型ROC 曲線Fig.2 ROC curve of MaxEnt model

2.2 伊洛河流域油松潛在空間分布格局

MaxEnt 模型的模擬結(jié)果表明，伊洛河流域內(nèi)模擬油松潛在適生區(qū)的面積總計為9 075 km2，約占伊洛河流域總面積的46.98%，說明在伊洛河流域油松可以大量分布。其中，流域內(nèi)油松的高適生區(qū)面積約為3 064 km2（表3），占流域總面積的15.86%，平均適生生長概率為0.866，中位數(shù)為0.878（圖3），主要位于伊洛河流域的中南部，呈連續(xù)帶狀分布，一小部分呈離散狀分布于伊洛河流域的西北部（圖4）；中適生區(qū)面積約為3 015 km2，占流域總面積的15.61%，平均適生生長概率為0.541，中位數(shù)為0.536，主要分布在高適生區(qū)邊緣，分布格局與高適生區(qū)基本一致，主要分布在流域南部，呈連續(xù)帶狀分布；低適生區(qū)面積約為2 996 km2，占流域總面積的15.51%，平均適生生長概率為0.300，中位數(shù)為0.299，在流域西南部呈離散帶狀分布，與中適生區(qū)交界較大，分布范圍較廣。低、中、高適生區(qū)均以欒川縣、盧氏縣、洛南縣等為主要代表區(qū)域。而非適生區(qū)面積約為10 241 km2，占流域總面積的53.02%，平均適生生長概率為0.047，中位數(shù)為0.022，主要分布在流域中部和東南部，其中有一小部分位于西部，分布區(qū)域多為連續(xù)型，分布于流域的大部分地區(qū)。

圖3 油松適生分布區(qū)存在概率的25%～ 75%分布范圍、最值、中位線和均值Fig.3 Probability for suitable distribution of P.tabuliformis

圖4 油松在伊洛河流域的適生分布區(qū)Fig.4 Suitable distribution area of P.tabuliformis in the Yiluo River Basin

表3 伊洛河流域范圍內(nèi)油松適生區(qū)分布面積Tab.3 Areas of suitable distribution of P.tabuliformis in the Yiluo River Basin

2.3 影響伊洛河流域油松分布的環(huán)境主導(dǎo)因子

Jackknife 檢驗結(jié)果表明，僅使用單變量時，正則化訓(xùn)練增益由高到低前5 個環(huán)境變量依次為降水季節(jié)性變化（bio_15，0.890 6）、平均年降水量（bio_12，0.87 9）、土壤類型（Soil Type，0.875 4）、溫度季節(jié)性變化（bio_4，0.874 3）、最干月降水量（bio_14，0.854 4）（圖5A）；測試增益由高到低前5 個環(huán)境變量依次為平均年降水量（bio_12，1.050 2）、土壤類型（Soil Type，1.025 4）、降水季節(jié)性變化（bio_15，1.004 4）、最干月降水量（bio_14，0.995 7）、溫度季節(jié)性變化（bio_4，0.959 1）（圖5B）。在所有環(huán)境變量中，平均年降水量、土壤類型、海拔、降水季節(jié)性變化、等溫性的貢獻(xiàn)率依次為35%、24.8%、19.9%、9.2%、7.1%，位居貢獻(xiàn)率前五位。而置換重要性前五位依次為降水季節(jié)性變化、海拔、平均年降水量、等溫性、土壤類型，其重要性數(shù)值分別為28.9%、14.9%、14.8%、14.1%、10.8%（表4）。由響應(yīng)曲線結(jié)果表明，油松適宜生長土壤類型為黃褐土和棕壤，降水季節(jié)性變化范圍為70～ 74 mm，平均年降水量為710～ 850 mm，氣溫季節(jié)變化范圍為7.8～8.5 ℃，海拔范圍在1 156 m 以上（圖6）。

表4 相關(guān)環(huán)境變量的貢獻(xiàn)率、累計貢獻(xiàn)率和置換重要性Tab.4 Single and cumulative contribution rate and permutation importance of some environmental variables

圖5 環(huán)境變量訓(xùn)練和測試集的重要程度Fig.5 Importance of training and regularization of environmental variables for P.tabuliformis in Yiluo River Basin by Jackknife test

圖6 Maxent 模型降水季節(jié)性變化、溫度季節(jié)性變化、土壤類型、海拔和平均年降水量的響應(yīng)曲線Fig.6 Response curve of seasonal variation of precipitation,seasonal variation of temperature,soil type,elevation and mean annual precipitation

3 討論

3.1 MaxEnt 模型模擬效果

物種潛在分布模型多種多樣，比如基于統(tǒng)計分析的GLM、GAM、RF、DOMAIN、MaxEnt 模型和GARP、ANN 等模型[11]，其中，MaxEnt 模型能使用分類變量，具有源碼開放、操作簡便、模擬精度高的特點，被廣泛認(rèn)為是生態(tài)建模中性能最好的方法之一[12-20]。此外，油松是大型木本植物，與草本植物相比，生態(tài)位寬度大且保守程度較高。因此，本研究選用MaxEnt 模型來模擬油松的潛在分布。但MaxEnt 難以考慮生物相互作用和偏向于保守生態(tài)位，模型過于復(fù)雜時可能難以用生態(tài)學(xué)原理解釋。因此，本文經(jīng)過采樣偏差處理和變量相關(guān)性分析等方式來降低模型誤差，提高模型的預(yù)測精度。本文基于AUC 參數(shù)檢驗?zāi)Ｐ托阅?，受試者操作曲線（ROC）下的面積（AUC）不受閾值限制[21]。MaxEnt 評價結(jié)果表明，訓(xùn)練和測試AUC 分別為0.957 和0.935，說明模型模擬精度達(dá)到優(yōu)秀水平。但是預(yù)測結(jié)果與實際分布仍然具有一定的擬合誤差，其原因一部分可歸結(jié)為環(huán)境數(shù)據(jù)誤差。生態(tài)位預(yù)測需要充足的環(huán)境數(shù)據(jù)，如果一部分生態(tài)位處于收集的環(huán)境變量范圍以外，則模擬出來的是受空間限制而縮小的生態(tài)位，或模型未包含完整的重要環(huán)境變量時，模擬的生態(tài)位大于現(xiàn)實生態(tài)位，從而導(dǎo)致誤導(dǎo)性預(yù)測的產(chǎn)生。

3.2 MaxEnt 模型優(yōu)化

本文采用空間距離篩選法來處理樣本點，去除63%的采樣點，采樣偏差有多種處理方法，比如聚類分析、空間距離篩選法或空間限制法等[22]。另外，環(huán)境因子之間存在空間線性關(guān)系，尤其降水量、溫度環(huán)境因子相關(guān)的衍生變量之間共線性關(guān)系更為明顯。為了便于更好地解釋油松分布與預(yù)測因子之間的關(guān)系，本文借助Spearman和Pearson 相關(guān)性分析縮減變量，同時保證模型具有足夠的變量來描述油松的生態(tài)位。對于相關(guān)性分析，目前還沒有統(tǒng)一標(biāo)準(zhǔn)，一般閾值取值在0.6～ 0.9[23-24]，當(dāng)預(yù)測變量過多時，閾值可適當(dāng)提高，以選擇生態(tài)學(xué)意義較高的變量，同時保證篩選的變量都能反映物種的有效特征。但當(dāng)樣本數(shù)量小于模擬變量時，篩選標(biāo)準(zhǔn)應(yīng)側(cè)重于其生態(tài)學(xué)意義，而不是相關(guān)性[22]。對氣候因素來說，降水量是影響油松生長的關(guān)鍵因素，溫度次之。與氣候相比，地形、土壤等非氣候環(huán)境因素在物種分布的重要性增加，尤其對于植物，土壤因子在一定的區(qū)域尺度上對植物分布有重要影響，在我們的研究中，海拔、土壤類型都是油松分布的主導(dǎo)因子之一。油松為深根型樹種，根系主要分布在40～ 80 cm 土層，土壤結(jié)構(gòu)對植物生長有重要的影響。對其他土壤因子，土壤含水率、有機(jī)質(zhì)含量等與林分密度有關(guān)[25]，其他土壤養(yǎng)分與油松林齡有關(guān)[26]，且隨氣候變化波動較大，不宜作為預(yù)測油松林生態(tài)位的環(huán)境因子。

3.3 油松分布與環(huán)境變量的相互關(guān)系

對于氣候因素來說，降水量是影響油松生長的關(guān)鍵因素，溫度次之，與梁非凡、何如夢等[25-27]的研究結(jié)果相似。但在關(guān)于降水量、溫度的衍生變量方面，肖敏等[2]認(rèn)為，最暖冷月溫度差、最冷月平均溫度等是其主導(dǎo)因子。造成這種差異的原因，一方面，是數(shù)據(jù)來源、處理方式以及預(yù)測變量等建模數(shù)據(jù)不一致；另一方面，是預(yù)測參數(shù)調(diào)整的差異，MaxEnt 模型參數(shù)較多，優(yōu)化方案多樣，預(yù)報效果不一。

與氣候相比，地形、土壤等非氣候環(huán)境因素在物種分布的重要性增加[28]，尤其對于植物。土壤因子在一定的區(qū)域尺度上對植物分布有重要影響[29]。本研究發(fā)現(xiàn)，土壤類型是油松分布的主導(dǎo)因子之一，土壤結(jié)構(gòu)對植物生長有重要影響[31]。對其他土壤因子，油松林的土壤含水率、有機(jī)質(zhì)含量等與林分密度有關(guān)[32]，其他土壤養(yǎng)分與油松林齡有關(guān)[33]，且隨氣候變化波動較大，不宜作為油松林的預(yù)測因子。

4 結(jié)論

根據(jù)MaxEnt 構(gòu)建的模型模擬伊洛河流域內(nèi)油松潛在適生區(qū)面積總計為9 075 km2，其中，高適生區(qū)面積約為3 064 km2、中適生區(qū)面積約為3 015 km2、低適生區(qū)面積約為2 996 km2，低、中、高適生區(qū)分別以欒川縣、盧氏縣、洛南縣等為主要代表區(qū)域。影響油松適生區(qū)的主要環(huán)境因素包括，年平均氣溫等溫性、溫度季節(jié)性變化、平均年降水量、最干月降水量、降水季節(jié)性變化、土壤類型、坡度、坡向、海拔。在所有環(huán)境變量中，平均年降水量、土壤類型、海拔、降水季節(jié)性變化和等溫性的貢獻(xiàn)率分別為35%、24.8%、19.9%、9.2%和7.1%。而置換重要性前五位依次為降水季節(jié)性變化、海拔、年平均降水量、等溫性、土壤類型。建議在伊洛河流域下游地區(qū)的油松生長適生區(qū)內(nèi)應(yīng)該適當(dāng)增加油松林分的種植密度，建議在流域的西南地區(qū)（盧氏縣、欒川縣等地）適當(dāng)增加人工油松林種植點，尤其在河網(wǎng)密度較大，水土流失嚴(yán)重的地區(qū)，而在流域南部和中部建立油松林保護(hù)基地。