河北省北部華北落葉松人工林立地指數空間分布預測*

李文博 呂振剛 黃選瑞 張志東

1.河北農業大學林學院 河北省林木種質資源與森林保護重點實驗室 保定 071000； 2.河南省林業調查規劃院 鄭州 451200； 3.華中農業大學資源與環境學院 武漢 430070)

華北落葉松(Larixprincipis-rupprechtii)是河北省重要用材樹種之一，主要分布于河北省北部地區，在維持該區木材生產、水土保持和森林游憩等方面發揮著重大作用(高敏等， 2017)。立地條件是林木生長發育和森林發展的重要基礎，研究發現，氣候、地形和土壤等立地條件會影響華北落葉松人工林生長及其適應性，進而影響立地生產力(王云霓等， 2013； 趙匡記， 2015)。當前，在林分水平關于華北落葉松人工林生產力的研究較多(王冬至等， 2015； 趙匡記， 2015)，但對于河北省北部地區華北落葉松人工林基于較大空間尺度立地生產力的確定及其影響因子的研究較少，不利于該地區華北落葉松人工林可持續經營(Bueisetal.， 2017)。

立地指數(site index, SI)指林分在基準年齡時的優勢木高，其實質是將任意林分的樹高最高值(林分優勢木高)標準化至同一年齡的生長曲線(Littkeetal.， 2016)，因優勢木高生長與材積生長密切相關且受林分結構和森林經營活動影響很小，因此SI被廣泛用作同齡純林的立地生產力評估指標。SI能夠反映出區域土壤、氣候等固有立地特征差異(McKenneyetal.， 2003； Skovsgaardetal.， 2008； Bueisetal.， 2017)，而氣候、土壤和地形等因子的一些不利變化也會導致森林立地生產力下降并對森林生態系統造成較大威脅(Bravo-Oviedoetal.， 2011； Littkeetal.， 2016； Sabatiaetal.， 2014)，因此一些學者(Skovsgaardetal.， 2008； Noordermeeretal.， 2018)基于氣候、地形和土壤等環境因子間接預測森林立地生產力。Nigh等(2004)建立不同氣候變量組與SI的多元回歸模型，得出不同樹種SI對氣候變量的響應程度不同； Bravo-Oviedo等(2011)通過氣候、土壤因子與地中海油松(Pinuspinaster)SI的多元線性回歸關系分析發現，在大的空間尺度，氣候是該樹種生產力的主要驅動力，但在同一氣候區內，SI差異與土壤密切相關； Jiang等(2015)構建樹種SI與氣候和土壤不同組合的隨機森林(random forest，RF)模型，認為基于氣候和土壤的SI模型預測精度稍高于僅基于氣候或土壤因子的SI模型。由此可知，SI預測精度受樹種類型、研究區域和尺度、模型預測變量等因素影響，準確預測華北落葉松人工林SI及其空間變異規律，需綜合分析SI與氣候、地形和土壤等環境因子的關系。

近年來，隨著林地需求的不斷增加，預估森林立地生產力空間變異對森林經營和保護規劃愈發重要(McKenneyetal.， 2003)。3S技術提高了較大區域尺度環境因子數據的精度和時效性，解決了傳統調查成本高、樣本獲取量少的問題(Bravo-Oviedoetal.， 2011)，并能夠繪制出信息明晰、分辨率較高的SI空間分布圖，可為適地適樹和制定精準提高森林質量的營、造林方案提供保證(Falkowskietal.， 2009； Parresoletal.， 2017)，因此很多研究基于3S多源空間環境因子數據構建了多元線性回歸(multiple linear regression，MLR)(Huangetal.， 2017)、隨機森林(RF)(Jiangetal.， 2015)、神經網絡回歸(Hlasnyetal.， 2017)等模型預測區域SI及其與環境因子的關系，并繪制了SI空間分布圖。地統計學方法同樣被用于森林立地生產力預測和繪圖，如Fayad等(2016)基于機載和星載激光雷達提取的樹高數據，構建與環境因子的回歸克里格(regression Kriging, RK)模型預測法屬圭亞那熱帶雨林高度并繪制了其空間分布圖； Benitez等(2016)基于植被指數和土壤數據采用地理加權回歸克里格(geographical weighted regression Kriging, GWRK)模型預測并繪制了亞馬遜熱帶雨林地上生物量空間分布圖。以往研究基于不同環境因子采用多種模型預測森林立地生產力的空間變異并對其結果進行可視化，但因林分類型、研究區域和尺度、模型預測變量等因素不同，模型預測精度存在較大差異。

為準確預測華北落葉松人工林SI及其空間變異規律，本研究基于河北省北部地區多源環境因子數據，篩選主要環境因子，應用MLR、RF、RK和GWRK模型預測研究區華北落葉松人工林SI并繪制其空間分布圖，旨在解決3個問題： 1) 哪些環境因子是影響研究區華北落葉松人工林SI的主要因子； 2) 如何在MLR、RF、RK和GWRK 4種預測模型中選擇最優預測模型，模型精度表現如何； 3) 華北落葉松人工林SI的空間分布格局及其對主要環境因子的響應機制是什么。

1 研究區概況

研究區位于河北省承德市圍場滿族蒙古族自治縣(41°31′—42°37′N，116°41′—118°20′E)(圖1)，北接內蒙古，地處內蒙古高原和冀北山地的過渡帶，海拔663～2 042 m，分壩上和壩下2種自然地貌類型。壩上氣候寒冷，年均溫-1.2～4.9 ℃，年均降水量445 mm； 壩下四季分明，年均溫-0.4～6.9 ℃，年均降水量441 mm，區域降水少且多集中于6—9月。受地形和氣候影響，形成了褐土、棕壤、灰色森林土、風沙土和草甸土等土壤類型。研究區喬木主要有華北落葉松、樟子松(Pinussylvestrisvar.mongolica)、云杉(Piceaasperata)、蒙古櫟(Quercusmongolica)等； 灌木主要有沙棘(Hippophaerhamnoides)、山刺玫(Rosadaverica)、稠李(Prunuspadus)等； 草本植物主要有蒲公英(Herbataraxaci)、苔草(Carextristachya)、老鸛草(Geaniumdaharicumvar.alpinum)等。

圖1 研究區地理位置及SI樣地點

2 數據收集與研究方法

2.1 數據收集

2.1.1 立地指數 采用塞罕壩機械林場和木蘭國有林場管理局森林資源二類調查數據及近6年樣地解析木數據。調查樣木健康且均無損害、無蟲害和人為活動干擾，其中，森林資源二類調查數據包含優勢樹種、林齡、胸徑和樹高等測樹因子，研究主要選取優勢樹種為華北落葉松人工純林的小班樣地； 樣地解析木采取5株優勢木和5株平均木的方法，分別確定能夠代表每塊樣地林分水平的1株優勢木和1株平均木進行解析，按5年一齡階、2 m一分段整化。由于不同立地林地生產力評價結果會因基準年齡誤取而有所偏差(吳恒等， 2015)，故選取研究區內不同立地條件下20塊華北落葉松人工近熟林(≥ 30年且年齡為5的倍數)樣地的優勢木解析資料分析樹高生長過程。由圖2可知，優勢高生長在林齡20年后逐漸趨于穩定，30年時對應的樣本量最大，故確定基準年齡為30年。由于森林資源調查樹高為樣地平均木高，因此篩選研究區內林齡≥ 20年的92塊不同立地條件樣地解析木數據建立華北落葉松人工林樣地優勢木-平均木樹高轉換方程(表1、圖3)，將研究區內基準年齡為30年的華北落葉松人工林樣地平均木高轉換為立地指數，剔除±3倍標準差外的樣地數據，最終對1 179塊樣地數據進行描述性統計分析(表1)。

2.1.2 環境因子 氣候數據來源于世界氣候數據庫，空間分辨率30弧秒(約為1 km空間分辨率)，與研究區基礎地理數據配準、裁剪后產生包括氣溫、降水等19個生物氣候變量數據(Hijmansetal.， 2005； http:∥www.worldclim.org/)； 土壤數據主要來源于北京師范大學中國土壤數據集(Weietal.， 2013)，空間分辨率30弧秒。為探究土壤化學性質與立地指數的相關關系，采用研究區內土壤全氮、全磷等9個土壤因子柵格數據集； 從地理空間數據云(http:∥www.gscloud.cn/)下載研究區30 m分辨率DEM影像，影像預處理后再基于DEM生成坡向、坡度和海拔3個地形因子。

31個環境因子柵格數據均采用Asia_Lambert_Conformal_Conic投影坐標系，空間分辨率統一為30 m。為篩選出對SI影響大且因子間相關性低的變量，首先對31個環境因子和SI進行Pearson相關分析，得到顯著影響(P< 0.05)SI的15個環境因子(表2)； 然后計算環境因子的方差膨脹系數VIF，剔除VIF > 15的環境因子，篩選得到海拔(DEM)、最干月降水(BIO14)、土壤有機質(SOM)、土壤全氮(TN)和土壤全磷(TP)5個環境因子； 最后建立SI與5個環境因子的逐步回歸方程，進一步篩選出與SI顯著相關的DEM、BIO14、TN和TP共4個環境因子，VIF分別為3.01、2.64、1.48和1.29，不存在多重共線性，4個環境因子的空間分布如圖4所示。

圖2 優勢高生長變異系數

圖3 樣地優勢木-平均木樹高轉換方程

表1 森林調查數據描述性統計結果

表2 環境因子描述性統計①

2.2 研究方法

2.2.1 隨機森林(RF) RF是一種非參數決策樹學習算法，每株決策樹依賴于獨立同分布的隨機向量，最初是為分類和回歸問題設計的，在機器學習和數據挖掘方面應用廣泛(Breiman， 2001)。設響應變量(SI)為一維向量Y，環境因子為p維向量X，其中，Y可為連續、二進制、分類或生存變量集合，X可為連續或離散變量集合。本研究SI為連續變量，故實際上是在非參數回歸框架下求X=x(X∈x)時的響應函數h(x)，此估算基于訓練數據λ={(X1，Y1)，…，(Xn，Yn)}，其中，(Xi，Yi)與(X，Y)獨立同分布，h(x)可表示為：

h(x)=E[Y∣X=x]，

(1)

即RF回歸模型輸出為所有決策樹結果的平均值。模型基于R3.5.1調用“randomForest”包實現，決策樹株數(ntree)和樹節點(mtry)是模型的2個重要設置參數，針對RF解決回歸問題，Breiman(2001)建議mtry默認值為全部響應變量的1/3(取整)，ntree也代表重采樣參數，一般選取整體誤差率趨于穩定時對應的值(歐強新等， 2018)。

2.2.2 回歸克里格(RK) RK為通過建立響應變量與輔助變量之間的線性模型并對回歸殘差進行OK(ordinary Kriging)插值從而預測響應變量空間特征的統計模型(Hengletal.， 2007)。設樣地為S，SI為響應變量，則：

SIRK(S)=SIMLR(S)+εRK(S)；

(2)

(3)

式中： SIMLR(S)和εRK(S)分別為模型趨勢項和殘差項，趨勢項主要通過建立MLR模型獲得，將MLR模型殘差基于ArcGIS進行OK插值獲得殘差項；k為環境因子總數，k=4(圖4)；bi為第i個環境因子Fi的模型系數；IMLR為截距。

2.2.3 地理加權回歸克里格(GWRK) GWRK與RK模型的最大區別在于模型趨勢項不同，GWRK首先采用局部回歸方法構建模型趨勢項，然后將局部回歸殘差進行OK插值(楊順華等， 2015)，其表達式為：

SIGWRK(S)=SIGWR(S)+εGWRK(S)。

(4)

式中： SIGWR(S)為響應變量在樣地S處的地理加權回歸(geographical weighted regression，GWR)預測值。

構建GWR的目的是探究響應變量與輔助變量之間是否存在空間非平穩性關系，可由下式表示：

(5)

式中：Fi(S)為第i個環境因子在樣地S處的實測值；βi(S)為第i個環境因子在樣地S處的擬合系數；IGWR(S)為截距。

圖4 環境因子空間分布

在ArcGIS中建立GWR模型，對其局部回歸殘差進行OK插值，得到εGWRK(S)，二者相加即為GWRK模型最終預測結果。

RK和GWRK模型殘差空間相關性通過半變異函數模型進行分析，采用最優半變異函數模型的插值參數進行OK插值。

2.2.4 模型精度評價 模型驗證基于10折交叉驗證法，即每次隨機選取1 179塊樣地的90%作為模型訓練集，10%作為測試集。采用模型決定系數(R2)、均方根誤差(root mean square error，RMSE)、絕對平均誤差(mean absolute error，MAE)和赤池信息準則(Akaike information criterion，AIC)進行評價和選擇，計算公式如下：

(6)

(7)

(8)

AIC=-2lgL+2p。

(9)

2.2.5 統計分析 因子篩選、MLR和RF模型建立、RK和GWRK預測模型殘差半變異函數分析以及模型精度指標計算均基于R3.5.1 進行，地統計學的RK和GWRK預測模型實現及SI空間預測主要基于ArcGIS 10.2進行，SI與主要環境因子的相關和偏相關分析主要在R3.5.1軟件中完成。

3 結果與分析

3.1 模型精度評價

由MLR、RF、RK和GWRK模型預測值與實測值殘差平均分布(表3)可知，整體上，4種模型殘差均為正態分布(P>0.05)； MLR和RF模型的空間自相關性明顯(P<0.05)，RK和GWRK模型殘差空間分布隨機(P>0.05)。由表4可知，MLR、RF、RK和GWRK模型預測精度交叉驗證結果相差較大： MLR模型R2最低，RMSE、MAE和AIC最大，模型精度最低； GWRK模型精度最高； RF模型精度一般； RK與GWRK模型精度相近。相比MLR模型，RF、RK和GWRK模型的平均預測精度均有較大提升，R2分別提高0.45、1.00和1.07倍，RMSE、MAE和AIC依次呈下降趨勢； 地統計學模型顯著優于RF模型，基于空間局部回歸的GWRK模型對SI空間變異的預測能力優于基于空間全局回歸的RK模型(R2提高0.4%，RMSE、MAE和AIC分別降低0.5%、1.0%和1.6%)。分別選取10折交叉驗證中精度最高的MLR、RF、RK和GWRK模型預測研究區SI空間分布，雖然4種模型預測SI分布趨勢大致相同，但MLR與RF模型預測結果“鋸齒化”明顯，而RK和GWRK模型預測結果高低值邊緣曲線更加細膩，GWRK模型不同SI的小區域明顯增多(圖5)。綜上所述，選擇GWRK作為研究區SI的最終預測模型。

表3 模型殘差描述性統計和正態性檢驗

3.2 SI與主要環境因子的空間分布格局

由表5可知，Moran’I的分布范圍為0.28～0.84，表明模型系數具有顯著空間自相關性(Z>2.58，P<0.01)，環境因子對SI的影響隨地理位置改變而變化。由GWRK模型預測的SI空間分布(圖5D)可知，立地指數整體上西北高、東南低，與DEM和BIO14的空間分布格局一致(圖4A、B)，表明DEM和BIO14共同主導SI的東西差異，但DEM空間局部“聚類”較BIO14更顯著(Moran’I: 0.65 > 0.38)，DEM對SI東西變異的主導性更強； TP和TN的空間分布格局較復雜，除西部邊緣及中南部較高外，均呈南低北高的空間分布格局(圖4C、D)，二者共同主導SI的南北差異。

表4 預測模型交叉驗證統計量

將SI自動劃分為Ⅰ～Ⅴ級，分別代表高(14.66～20.43 m)、較高(12.53～14.66 m)、中(10.39～12.53 m)、較低(8.40～10.39 m)和低(< 8.40 m)立地生產力(圖5D)。第Ⅰ等級占研究區總面積的10.7%，主要分布于西北部壩上地區； 第Ⅱ等級占研究區總面積的15.3%，絕大部分分布于西部，北部分布適中，極少數(約占第Ⅱ等級的1.1%)分布于東部； 第Ⅲ等級占研究區總面積的25.7%，主要分布于第Ⅰ、Ⅱ等級和Ⅳ、Ⅴ等級邊緣區域及東部； 第Ⅳ、Ⅴ等級占研究區總面積比例較大(48.4%)，集中于東部和南部的壩下地區。研究區華北落葉松人工林立地生產力區域兩極分化較為明顯，西北壩上和東南壩下地區SI差距較大，為7.99～15.43 m。

圖5 不同預測模型SI空間分布

表5 GWR模型描述性統計

3.3 SI與主要環境因子的關系

研究區華北落葉松人工林SI預測值與主要環境因子的相關分析(表6)表明，在大的空間尺度上，SI與DEM、BIO14、TN和TP呈顯著正相關關系(P< 0.01)； 偏相關分析表明，控制其他環境因子不變，各立地生產力等級內SI均對BIO14的變化響應較弱； 第Ⅰ、Ⅴ等級SI主要受限于TN，第Ⅱ、Ⅲ等級SI與DEM的變化均呈顯著正相關關系(P< 0.01)，而第Ⅱ、Ⅲ等級SI分別隨TP、TN增加而增大(P< 0.01)。 結合圖5、6可知，華北落葉松人工林喜海拔較高(約> 1 450 m)、土壤TP含量(0.07～0.08 g·kg-1)以及TN含量(0.2～0.3 g·kg-1)適中的區域，土壤TN、TP含量過高/低均會限制華北落葉松人工林生長。

表6 各立地生產力等級SI與環境因子的偏相關關系①

圖6 SI與其敏感環境因子的關系

4 討論

本研究探討研究區SI隨環境因子的變化規律，進一步驗證了重要環境因子作為輔助變量預測立地生產力的可行性(Bravo-Oviedoetal.， 2011； Jiangetal.， 2015)； 但不同立地生產力因研究區域以及林分類型等差異，篩選出的重要環境因子不盡相同(Huangetal.， 2017； 雷相東等， 2018)。綜合考慮SI樣本數量和環境因子數量，本研究最終選取與立地指數相關性高(P< 0.01)、彼此之間相關性低(VIF < 4)的海拔、最干月降水、土壤全氮和全磷作為影響研究區立地指數的重要環境變量，最干月降水在數值上存在較大變異(CV=43.75%)，與立地指數呈顯著正相關關系(P< 0.01)，在空間上也具有較強空間相關性(Moran’I=0.54，P< 0.01)。為進一步探討其對立地指數的影響，本研究保留該解釋變量，剔除多重共線性問題嚴重(VIF > 15)的生長季溫度和降水等氣候因子。

立地生產力模型能夠應用于森林經營實踐至少需滿足2個條件： 解釋立地指數變異至少50%； 輔助變量易于獲取(Blythetal.， 1981)。本研究除MLR模型外，其他3種預測模型均滿足該條件，能夠解釋立地指數變異的54.4%～78.0%，其中地統計學模型(RK和GWRK)的解釋能力更強(表4)。RK和GWRK模型預測殘差均采用屬于線性克里格的普通克里格進行插值，其主要差異是趨勢項建模方法不同，RK建立的是環境因子與立地指數的全局回歸模型(MLR)，而GWRK建立的是局部回歸的地理加權回歸模型(GWR)。由于MLR模型僅為整個研究區數據的平均化呈現，丟棄了大量關系空間變異和模型性能的潛在有價值信息(Fotheringhametal., 1998)，而GWR模型基于連續空間權重函數估測除觀測樣地外其他位置參數，克服了全局回歸模型中不同空間位置環境因子權重相同的問題(Georganosetal.， 2017)，因此GWRK模型精度更高(表4)； 但也有研究得出相反結論(Jinetal.， 2016)，可能是受數據完整性和空間尺度大小等因素的影響。本研究為保證所有模型輸入變量一致，只選取4個主要環境因子進行建模，雖然RF模型適用于數據維度較高的情況，且不需要預先選擇解釋變量(Farrellyetal.， 2011)，但當解釋變量為與立地指數顯著相關(P< 0.01)的15個環境因子時，其預測精度(R2=0.671, RMSE=1.817, MAE=1.346, AIC=281.824)低于地統計學模型，并不是最優模型。相比RF和地統計學模型，MLR模型只能解釋預測變量的部分方差變異，與以往研究結論(Zhangetal.， 2017)一致。因此本研究認為，當預測變量具有顯著空間變異規律或預測變量與解釋變量之間并不是簡單的全局線性關系時，應避免使用MLR模型。

立地生產力整體分布趨勢與海拔和最干月降水分布規律一致，均由西北向東南遞減(圖5D)。研究區地處華北針葉林帶，森林覆蓋面積大，而東西海拔變化再分配影響林分生長的光照和地表土壤水熱條件，進而對土壤養分及酶活性、土壤微生物、凋落物質量和根系產生影響(Fiereretal.， 2011； Wangetal.， 2009)，導致華北落葉松人工林生長差異。受地形影響，研究區最干季(一般為冬季)均溫-17.88～-7.07 ℃，西北低、東南高。溫帶落葉喬木光合作用最低溫度為-3～-1 ℃，而冬季木本植物的越冬器官-25 ℃仍未停止呼吸(潘瑞熾， 2008)，因此，在最干月華北落葉松呼吸作用會消耗較多的儲存養分以維持基本生命活動，且在低溫脅迫條件下可能會造成細胞缺水(Bray， 1997)。水是植物細胞中最重要的分子，降水會影響整個土壤-植物-大氣連續體(soil-plant-atmosphere continuum，SPAC)的物質循環(Pengetal.， 2007； Philip， 1996)，而最干月降水量大的地區能夠增強華北落葉松人工林對低溫干旱條件的耐受能力(Bray， 1997)，因此立地指數高，但主要影響宏觀區域尺度(表6)。土壤氮磷是植物生長發育的兩大限制因子和必需的礦質元素(盧同平等， 2017)，土壤全氮和全磷含量增加在一定程度上可促進植物吸收氮磷營養元素，而植物氮磷含量反映植物養分的地理分布格局和生物地球化學循環的區域特征(樊江文等， 2014)，是對區域立地因子、生物群落以及人類活動等生物和非生物因子長期適應的結果(盧同平等， 2017)，如西北邊緣地區的華北落葉松人工林生長受海拔和土壤全磷含量的共同限制(圖5、表5)。溫帶華北落葉松人工林土壤全氮、全磷含量過高或過低均會導致生態系統中氮磷的不平衡輸入，可能誘發植物從氮限制向磷限制轉變(Dengetal.， 2016)，不利于華北落葉松生長發育，本研究進一步驗證了該結果(表6、圖6)。總體來看，華北落葉松人工林立地生產力等級高的地區海拔約高于1 400 m(圖6)，且在土壤全氮和全磷含量適中的地區競爭優勢更為突出。在森林經營實踐中，管理者應充分認識到立地指數是林木生長特性與多個環境因子綜合作用的結果(Littkeetal.， 2016)，從而采取有效的經營措施提高立地質量，促進樹木生長。

5 結論

空間明晰的森林立地生產力圖的繪制以及立地生產力與重要環境因子之間關系的確定一直是森林經營實踐中的熱點。本研究建立的立地指數與主要環境因子之間的地理加權回歸克里格模型，可用于繪制立地指數空間分布圖，探討研究區華北落葉松人工林立地生產力對主要環境因子的響應機制。在景觀尺度上，華北落葉松人工林立地指數對主要環境因子響應較強，海拔、最干月降水較高以及土壤氮磷含量適中的研究區西北部，是華北落葉松人工林適生區，而海拔、最干月降水較低以及土壤氮磷含量不平衡的東、南部，華北落葉松人工林的適應性較差。因此，建議在東、南部地區的森林經營實踐中，可通過控制林分密度(約2 170株·hm-2)(任麗娜等， 2012)、營造混交林來緩解氣候變化和海拔對樹木生長的限制，或向土壤中添加適量的氮磷肥以平衡土壤氮磷比，從而促進華北落葉松人工林生長以及生產力提高。