基于MaxEnt和GIS的浙江省杉木良種在湖北省適宜引種生態(tài)區(qū)研究

胡 超，于 靜

(1.湖北省林業(yè)局林木種苗管理總站，武漢 430079；2.嶺南生態(tài)文旅股份有限公司，武漢 430062)

林木良種有其適宜的生態(tài)區(qū)域，如果自然條件不適宜，再好的良種也達不到豐產(chǎn)、穩(wěn)產(chǎn)目標(biāo)。過去我國林木良種不適應(yīng)引種區(qū)自然條件而造成巨大損失的教訓(xùn)深刻，如20世紀(jì)70年代，各地在油茶(Camelliaoleifera)生產(chǎn)發(fā)展過程中盲目引種，導(dǎo)致幼林生長不良、成林產(chǎn)量很低，在人力、物力等方面都造成了不小損失[1]。傳統(tǒng)的林木良種引種適宜生態(tài)區(qū)憑主觀經(jīng)驗判斷較多，如宜林范圍內(nèi)每個氣候帶和垂直氣候帶都分布著特有類型的森林植被。經(jīng)緯度由北向南，由西向東調(diào)運范圍大于相反方向的范圍，海拔高度不超過300～500 m，但是，1958年，湖北引種廣東、福建馬尾松(Pinusmassoniana)種子成功，用事實改變了過去專家認為“馬尾松南種北移的幅度不能超過2°～3°”的定論[2]。1979年李傳志論證認為馬尾松一次北移6°～7°育苗可以成功[2]。所以溫度、降水、土壤等主要環(huán)境因子相似，即為林木良種同一適宜引種生態(tài)區(qū)。

杉木(Cunninghamialanceolata)是包括湖北省在內(nèi)的我國南方主要造林樹種，杉木良種數(shù)量較多，生長周期長，其引種不能像農(nóng)作物良種一樣都進行引種試驗，傳統(tǒng)的良種引種適宜生態(tài)區(qū)也憑主觀經(jīng)驗判斷較多。為了充分發(fā)揮湖北省造林地立地條件的優(yōu)勢，給杉木良種“適地適樹”造林，經(jīng)營管理上使杉木生產(chǎn)達到“速生、豐產(chǎn)、優(yōu)質(zhì)”奠定良好基礎(chǔ)。我們充分科學(xué)利用浙江省優(yōu)良杉木品種資源，使用ARCGIS和MaxEnt生態(tài)學(xué)模型精準(zhǔn)分析了杉木良種同一適宜引種生態(tài)區(qū)[3-5]。

1 數(shù)據(jù)與方法

1.1 基本信息來源

浙江省審定推廣的杉木良種信息概況如表1所示，來源于湖北省林業(yè)局林木種苗管理總站。

表1 浙江省杉木產(chǎn)區(qū)已審定的杉木良種

34個環(huán)境因子數(shù)據(jù)獲取于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)、中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心、國家青藏高原科學(xué)數(shù)據(jù)中心、中國西部環(huán)境與生態(tài)科學(xué)數(shù)據(jù)中心(表2)。中國行政區(qū)劃數(shù)據(jù)、中國海拔高度(DEM)數(shù)據(jù)獲取于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心和湖北省林業(yè)調(diào)查規(guī)劃院。

表2 浙江省杉木良種在湖北省同一適宜引種生態(tài)區(qū)環(huán)境因子

1.2 信息數(shù)據(jù)處理

1.2.1 分布數(shù)據(jù)處理 根據(jù)浙江省杉木產(chǎn)區(qū)現(xiàn)有審定杉木良種，本研究以100 m×100 m=10 000 m2即1 hm2為單元，用34個環(huán)境因子精準(zhǔn)劃分浙江省杉木良種在湖北省同一適宜引種生態(tài)區(qū)。

為避免樣點數(shù)據(jù)在某個地理空間上過度聚集，在浙江省杉木適生范圍內(nèi)，用Arcgis10的Create fishnet 工具生成空間為40行×40列的格網(wǎng)數(shù)據(jù)，以1個格網(wǎng)作為1個采樣單元對杉木良種的分布數(shù)據(jù)進行采樣(圖1)[6]。根據(jù)選育單位確定的杉木良種適宜的自然地理環(huán)境條件范圍，在Excel表中，剔除高程小于0 m、高程大于1 500 m、土壤厚度小于30 cm和異常值的采樣點，全部采樣分布記錄共725條。按照MaxEnt軟件的“Samples”要求整理數(shù)據(jù)，將分布點以“物種+經(jīng)度+緯度(西經(jīng)、南緯的值為負，經(jīng)緯度為十進制小數(shù)格式)”另存為CSV格式文件。

圖1 浙江省杉木產(chǎn)區(qū)采樣點分布圖

1.2.2 環(huán)境因子處理 地形因子(經(jīng)度、緯度、高度)與環(huán)境因子有較好的回歸關(guān)系，利用中國2 160個基本、基準(zhǔn)地面氣象觀測站的觀測數(shù)據(jù)，推算模擬無測站區(qū)域的環(huán)境資源分布情況。建立Bio1～Bio10、Bio13～Bio21、Bio24～Bio27等23個環(huán)境因子的空間分布模型，其表達式為:

Y=f(λ，φ，h)+ε

式中，Y為環(huán)境因子要素，λ為經(jīng)度，φ為緯度，h為海拔高度(m)，函數(shù)f(λ,φ,h)為氣候?qū)W方程，ε為殘差項，可視為小地形因子(坡度、坡向等)及下墊面對環(huán)境的影響。將f(λ,φ,h)展成三維二次趨勢面方程[7]。

f(λ,φ,h)=b0+b1λ+b2φ+b3h+b4λφ+b5φh+b6λh+b7λ2+b8φ2+b9h2

式中b0～b9為待定系數(shù)，利用SAS9.4建立逐步回歸優(yōu)化回歸模型，模擬23個環(huán)境因子的宏觀趨勢項，分別建立23個環(huán)境因子的小網(wǎng)格推算模型(表3)。

表3 環(huán)境因子的小網(wǎng)格推算模型

在中國海拔高度(DEM)數(shù)據(jù)支持下，在ArcGIS10里，用23個環(huán)境因子的小網(wǎng)格推算模型，將環(huán)境因子Bio1～Bio10、Bio13～Bio21、Bio24～Bio27分別插值為100 m×100 m 網(wǎng)格的基礎(chǔ)數(shù)據(jù)[8-9]。用IDW法分別插值其殘差項為100 m×100 m 網(wǎng)格的修正數(shù)據(jù)。用Spatial Analyst工具→數(shù)學(xué)→邏輯→加，將每個環(huán)境因子的基礎(chǔ)數(shù)據(jù)和修正數(shù)據(jù)疊加相加為環(huán)境因子?xùn)鸥駭?shù)據(jù)。23個環(huán)境因子?xùn)鸥駭?shù)據(jù)用投影柵格工具統(tǒng)一為地理坐標(biāo)系D_WGS_1984。以湖北省和浙江省矢量邊界為掩膜，裁剪出這23個環(huán)境因子?xùn)鸥駭?shù)據(jù)圖層。最后，用柵格轉(zhuǎn)ASCII工具將這23個環(huán)境因子?xùn)鸥駭?shù)據(jù)轉(zhuǎn)換保存為MaxEnt所需要的ASCII格式文件。

在ArcGIS10里，將下載的Bio11、Bio12、Bio22、Bio23、Bio28～Bio34等11個環(huán)境因子數(shù)據(jù)通過重采樣工具使其像元大小與Bio1～Bio10、Bio13～Bio21、Bio24～Bio27等23個環(huán)境因子一致[10]。11個環(huán)境因子數(shù)據(jù)統(tǒng)一為地理坐標(biāo)系D_WGS_1984。以湖北省和浙江省矢量邊界為掩膜，裁剪出這11個環(huán)境因子?xùn)鸥駭?shù)據(jù)圖層。最后，用柵格轉(zhuǎn)ASCII工具將這11個環(huán)境因子?xùn)鸥駭?shù)據(jù)轉(zhuǎn)換保存為MaxEnt所要求的ASCII格式文件。

1.3 模型構(gòu)建

1.3.1 MaxEnt軟件建模

(1)物種數(shù)據(jù)

將之前導(dǎo)出的杉木良種分布數(shù)據(jù)(csv格式)的文件，通過Browse加載到MaxEnt軟件“Samples”模塊。

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

將34個ASCII格式的環(huán)境數(shù)據(jù)加載到MaxEnt軟件“Environmental layers”模塊。

(3)參數(shù)設(shè)置

使用auto features選項，根據(jù)自動特征規(guī)則進行計算，所有的要素類型都將用到。結(jié)果以comulative類型和ASCII格式輸出，并定義其輸出位置。設(shè)置界面的選擇 settings里‘Random test percentage、設(shè)置為25，隨機選取75%的樣本點數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)[11]，settings中replicates本試驗選擇3次重復(fù)作為平行試驗，最大迭代次數(shù)設(shè)為500次，收斂閥值設(shè)為0.000 01，取值范圍0～100[12]。選擇‘Do jackknife to measure variable importance’衡量所有變量的重要性，MaxEnt軟件分別對每一個環(huán)境影響因子進行刀切圖繪出。

1.3.2 ROC曲線繪制 繪制響應(yīng)曲線(response curves)評價模型精度。ROC 曲線以真陽性率為縱坐標(biāo)(敏感性，實際存在且被預(yù)測為存在的比率)，以假陽性率(1-特異性，實際不存在但被預(yù)測為存在的比率)為橫坐標(biāo)，AUC值指 ROC曲線與橫坐標(biāo)圍成的面積值，值域為0～1。AUC值越大表示與隨機分布相距越遠，環(huán)境因子變量與預(yù)測的杉木良種同一適宜引種生態(tài)區(qū)之間的相關(guān)性越大，即模型預(yù)測效果越好，反之說明模型預(yù)測效果越差。AUC值在0.5～0.6、0.6～0.7、0.7～0.8、0.8～0.9、0.9～1分別表示模擬效果失敗、較差、一般、好、非常好[4，8]。34個環(huán)境因子預(yù)測模型的訓(xùn)練樣本和測試樣本的AUC值達到0.842和0.806(圖2)，AUC均值在0.8～0.9之間，說明模型預(yù)測效果好。

1.4 模型優(yōu)化

在使用MaxEnt模型進行較大空間范圍的杉木良種同一適宜引種生態(tài)區(qū)預(yù)測時，如果選擇了較多數(shù)量的環(huán)境因子，將導(dǎo)致MaxEnt模型的運行效率降低。因為環(huán)境因子變量過多、變量空間共線性過強，將導(dǎo)致模型的復(fù)雜性增加，隨機誤差增大。所以，過多低貢獻率的環(huán)境因子變量會導(dǎo)致模型運行結(jié)果的準(zhǔn)確性降低。因此，需要對環(huán)境因子進行篩選或降維[13]。

1.4.1 篩選貢獻率高的環(huán)境因子變量 在34個環(huán)境因子中，對于杉木良種同一適宜引種生態(tài)區(qū)分布貢獻較大的環(huán)境因子變量有：Bio2、Bio5、Bio7、Bio9、Bio14、Bio16、Bio18、Bio20、Bio22、Bio25，累計貢獻率為95%。Bio1、Bio3、Bio4、Bio6、Bio8、Bio10～Bio13、Bio15、Bio17、Bio19、Bio21、Bio23、Bio24、Bio26～Bio34等24個環(huán)境變量的貢獻率都小于1%(表4)，對杉木的種植分布影響有限，對這24個環(huán)境因子變量進行剔除[14]。

表4 各環(huán)境因子變量的貢獻率

1.4.2 篩選正規(guī)化訓(xùn)練增益高的環(huán)境因子變量 刀切法(jackknife test)測定各環(huán)境因子變量權(quán)重。刀切法就是每次都忽略一個環(huán)境因子變量，然后基于剩下的環(huán)境因子變量來對杉木良種同一適宜引種生態(tài)區(qū)進行預(yù)測，最后MaxEnt繪出環(huán)境因子變量重要性的柱形圖。紅色條帶表示所有變量的貢獻；深藍色條帶表示該變量的貢獻，條帶越長，說明該變量越重要；淺藍色的條帶長度表示除該變量以外，其他所有變量組合的貢獻。Bio2、Bio5、Bio7、Bio9、Bio14、Bio16、Bio18、Bio20、Bio22、Bio25等10個貢獻較大的環(huán)境因子變量中，Bio7對應(yīng)的深藍色條帶很短(圖3)，說明它本身的增益值小于0.1，表明它對預(yù)測杉木良種同一適宜引種生態(tài)區(qū)并不是重要環(huán)境因子變量，所以，剔除環(huán)境因子變量Bio7。

圖3 刀切法的環(huán)境因子變量重要性分析

1.4.3 篩選多重共線的環(huán)境因子變量 用GIS軟件的值提取至點工具提取有效分布點的環(huán)境因子變量數(shù)值，用SPSS軟件對貢獻較大的Bio2、Bio5、Bio9、Bio14、Bio16、Bio18、Bio20、Bio22、Bio25等9個主導(dǎo)環(huán)境因子進行Spearman相關(guān)分析(表5)，檢驗環(huán)境因子變量之間的多重共線性。Bio16與Bio14的相關(guān)系數(shù)|r|≥0.8，對比初始模型中二者的貢獻率，二者的貢獻率相同，Bio14置換重要性較小，所以，剔除置換重要性較低的變量Bio14，提高模型模擬的精度[14]。

表5 關(guān)鍵環(huán)境因子變量的相關(guān)系數(shù)

2 結(jié)果與分析

2.1 模型的有效性

用剩余的Bio2、Bio5、Bio9、Bio16、Bio18、Bio20、Bio22、Bio25等8個主導(dǎo)環(huán)境因子變量重新建模，重建模型的訓(xùn)練樣本和測試樣本的AUC值達到0.822和0.813(圖4)，AUC均值在0.8～0.9之間，表明重建模型適用性及模擬精度均好，與主導(dǎo)環(huán)境因子變量之間的相關(guān)性大，預(yù)測同一適宜引種生態(tài)區(qū)的結(jié)果好，可以據(jù)此進行引種推廣。

圖4 重建模型的ROC曲線分析及AUC值

2.2 杉木良種同一適宜引種生態(tài)區(qū)預(yù)測及適生等級劃分

MaxEnt進行3次重復(fù)試驗，選取重復(fù)試驗中，AUC值最高的圖層導(dǎo)人ArcGIS軟件進行適宜等級劃分和可視化表達(圖5)。MaxEnt模型輸出的數(shù)據(jù)為ASCⅡ格式，用ArcGIS的ASCII to Raster功能，輸出數(shù)據(jù)類型選FLOAT，使該結(jié)果可在 ArcGIS中顯示[14]。利用“Reclassify”功能, 劃分分布值等級及相應(yīng)分布范圍, 并使用不同顏色表示，劃分標(biāo)準(zhǔn)為：存在概率<0.05 為不適生區(qū); 0.05≤存在概率<0.33為低適生區(qū)；0.33≤存在概率<0.66為中適生區(qū); 存在概率≥0.66為高適生區(qū)[4，8]。整體來看，浙江省杉木產(chǎn)區(qū)的杉木良種在湖北省無高適生區(qū)和中適生區(qū)，只有低適生區(qū)和不適生區(qū)，可以引種的低適生區(qū)面積為986 196 hm2，主要分布在：鄂東的黃梅縣、蘄春縣、英山縣、羅田縣、麻城市、陽新縣、大冶市、咸安區(qū)、赤壁市、通山縣、崇陽縣和通城縣；鄂中的廣水市和大悟縣；鄂西的長陽縣、五峰縣、神龍架和竹溪縣。低適宜區(qū)域在引種杉木良種時，需要選擇適宜的小生境。

圖5 浙江省杉木良種在湖北省同一適宜引種生態(tài)區(qū)分布圖

2.3 杉木良種地理分布與生物氣候變量的關(guān)系

用刀切法(Jackknife Test)檢測8個主導(dǎo)環(huán)境因子變量對分布增益貢獻的結(jié)果(表6)表明：當(dāng)月最長連續(xù)降水日數(shù)(Bio18)對杉木分布的增益最大，當(dāng)月最長連續(xù)降水日數(shù)為12.906～19 d，分布值隨月最長連續(xù)降水日數(shù)的升高而增大(圖6)。年極大風(fēng)速(Bio25)也對杉木分布的影響較大，當(dāng)年極大風(fēng)速為21～28 m·s-1，分布值隨年極大風(fēng)速的升高而增大(圖7)。太陽輻射日均值(Bio2)也對杉木分布的影響較大，當(dāng)太陽輻射日均值為132～155 w·m-2，分布值隨太陽輻射日均值的升高而增大(圖8)。

表6 主導(dǎo)環(huán)境因子變量的貢獻率

圖6 月最長連續(xù)降水日數(shù)(Bio18)反饋曲線

圖7 年極大風(fēng)速(Bio25)反饋曲線

圖8 太陽輻射日均值(Bio2)反饋曲線

3 討論

林木良種引種若沒做到“適地適種”，則可能導(dǎo)致育苗和造林失敗。主要原因是不同地理、氣候和土壤等環(huán)境因子不盡完全相同，對該林木良種所要求的生態(tài)條件滿足程度不同。所以，只有在適生地域內(nèi)和適宜的立地條件下，選擇適合的林木良種，才能發(fā)揮林木良種造林的優(yōu)良特性，實現(xiàn)速生豐產(chǎn)。否則，林木良種表現(xiàn)不好，甚至不如一般的當(dāng)?shù)仄贩N。傳統(tǒng)林木引種是以單個樹種劃出同一適宜生態(tài)區(qū)。然而，隨著自然條件演變和科學(xué)技術(shù)發(fā)展，轉(zhuǎn)抗性基因育種、種間和遠緣雜交育種等林木遺傳改良工作廣泛開展，每年都有新的林木良種通過審定。在相同的立地條件下，同一樹種，不同良種之間的生長好壞存在明顯差別。為了獲得精準(zhǔn)的引種效果，本研究是以單個良種劃出同一適宜生態(tài)區(qū)，精準(zhǔn)預(yù)測浙江省杉木良種在湖北省同一適宜生態(tài)區(qū)。

基于MaxEnt生態(tài)位模型的同一適宜生態(tài)區(qū)研究中，環(huán)境因子數(shù)據(jù)常來自于世界氣候-全球氣候數(shù)據(jù)庫網(wǎng)站，僅有19個環(huán)境因子，空間分辨率僅為5arc-min[5,6,10,14]。為了獲得更精準(zhǔn)的引種效果，本研究選取34個重要環(huán)境因子，用中國2 160個基準(zhǔn)地面氣象觀測站的觀測數(shù)據(jù)，推算模擬無測站區(qū)域的環(huán)境資源分布情況，精準(zhǔn)預(yù)測浙江省杉木良種在湖北省同一適宜引種生態(tài)區(qū)。傳統(tǒng)的杉木良種引種同一適宜生態(tài)區(qū)都是以鄉(xiāng)鎮(zhèn)、縣、市、省等行政單位為單元。然而，影響杉木成活生長的光、熱、水、氣等環(huán)境因子，受太陽輻射、大氣環(huán)流的影響，而使水、熱在地面上呈地帶性的分布。由于山體起伏，垂直森林地帶在實際中并不都是連續(xù)的，而是由斷斷續(xù)續(xù)地呈孤島狀分布的地塊組成。為了獲得精準(zhǔn)的引種效果，本研究是以100 m×100 m=10 000 m2即1 hm2為單元，精準(zhǔn)預(yù)測浙江省杉木良種在湖北省同一適宜生態(tài)區(qū)。

4 結(jié)論

本研究基于ARCGIS數(shù)據(jù)和氣象與土壤等數(shù)據(jù)，運用MaxEnt生態(tài)位模型對浙江省目前推廣的13個杉木良種在湖北省在湖北省適宜引種生態(tài)區(qū)進行了研究分析，結(jié)果為浙江省杉木產(chǎn)區(qū)的13個杉木良種在湖北省無高適生區(qū)和中適生區(qū)，只有低適生區(qū)和不適生區(qū)，可以引種的低適生區(qū)面積為986 196 hm2，主要包括鄂東的12個縣、鄂中的2個縣市、鄂西的4個縣區(qū)，低適宜區(qū)域在引種杉木良種時，需要選擇適宜的小生境。

通過運用MaxEnt生態(tài)位模型對浙江省杉木良種在湖北省同一適宜引種生態(tài)區(qū)進行分析，證明了MaxEnt模型在林木良種引種應(yīng)用方面的可行性以及可信度，同時結(jié)合刀切法探討對杉木良種生長影響最顯著的環(huán)境因子，這對浙江省杉木良種適生性分析提供了更進一步的技術(shù)支撐?；贛axEnt和GIS精準(zhǔn)預(yù)測浙江省杉木良種在湖北省同一適宜引種生態(tài)區(qū)，可進一步避免引種浙江省杉木良種的盲目性，獲得良好的引種效果，值得推廣應(yīng)用。