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

1.河北農(nóng)業(yè)大學(xué)林學(xué)院 河北省林木種質(zhì)資源與森林保護(hù)重點(diǎn)實(shí)驗(yàn)室 保定 071000； 2.河南省林業(yè)調(diào)查規(guī)劃院 鄭州 451200； 3.華中農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院 武漢 430070)

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

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

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

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

1 研究區(qū)概況

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

圖1 研究區(qū)地理位置及SI樣地點(diǎn)

2 數(shù)據(jù)收集與研究方法

2.1 數(shù)據(jù)收集

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

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

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

圖2 優(yōu)勢(shì)高生長(zhǎng)變異系數(shù)

圖3 樣地優(yōu)勢(shì)木-平均木樹高轉(zhuǎn)換方程

表1 森林調(diào)查數(shù)據(jù)描述性統(tǒng)計(jì)結(jié)果

表2 環(huán)境因子描述性統(tǒng)計(jì)①

2.2 研究方法

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

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

(1)

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

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

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

(2)

(3)

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

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

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

(4)

式中： SIGWR(S)為響應(yīng)變量在樣地S處的地理加權(quán)回歸(geographical weighted regression，GWR)預(yù)測(cè)值。

構(gòu)建GWR的目的是探究響應(yīng)變量與輔助變量之間是否存在空間非平穩(wěn)性關(guān)系，可由下式表示：

(5)

式中：Fi(S)為第i個(gè)環(huán)境因子在樣地S處的實(shí)測(cè)值；βi(S)為第i個(gè)環(huán)境因子在樣地S處的擬合系數(shù)；IGWR(S)為截距。

圖4 環(huán)境因子空間分布

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

RK和GWRK模型殘差空間相關(guān)性通過半變異函數(shù)模型進(jìn)行分析，采用最優(yōu)半變異函數(shù)模型的插值參數(shù)進(jìn)行OK插值。

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

(6)

(7)

(8)

AIC=-2lgL+2p。

(9)

2.2.5 統(tǒng)計(jì)分析 因子篩選、MLR和RF模型建立、RK和GWRK預(yù)測(cè)模型殘差半變異函數(shù)分析以及模型精度指標(biāo)計(jì)算均基于R3.5.1 進(jìn)行，地統(tǒng)計(jì)學(xué)的RK和GWRK預(yù)測(cè)模型實(shí)現(xiàn)及SI空間預(yù)測(cè)主要基于ArcGIS 10.2進(jìn)行，SI與主要環(huán)境因子的相關(guān)和偏相關(guān)分析主要在R3.5.1軟件中完成。

3 結(jié)果與分析

3.1 模型精度評(píng)價(jià)

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

表3 模型殘差描述性統(tǒng)計(jì)和正態(tài)性檢驗(yàn)

3.2 SI與主要環(huán)境因子的空間分布格局

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

表4 預(yù)測(cè)模型交叉驗(yàn)證統(tǒng)計(jì)量

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

圖5 不同預(yù)測(cè)模型SI空間分布

表5 GWR模型描述性統(tǒng)計(jì)

3.3 SI與主要環(huán)境因子的關(guān)系

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

表6 各立地生產(chǎn)力等級(jí)SI與環(huán)境因子的偏相關(guān)關(guān)系①

圖6 SI與其敏感環(huán)境因子的關(guān)系

4 討論

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

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

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

5 結(jié)論

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