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

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

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

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

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

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

1 研究區(qū)概況

研究區(qū)位于河北省承德市圍場滿族蒙古族自治縣(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樣地點

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

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

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

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

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

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

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

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

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

2.2 研究方法

2.2.1 隨機森林(RF) RF是一種非參數(shù)決策樹學(xué)習(xí)算法，每株決策樹依賴于獨立同分布的隨機向量，最初是為分類和回歸問題設(shè)計的，在機器學(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ù)回歸框架下求X=x(X∈x)時的響應(yīng)函數(shù)h(x)，此估算基于訓(xùn)練數(shù)據(jù)λ={(X1，Y1)，…，(Xn，Yn)}，其中，(Xi，Yi)與(X，Y)獨立同分布，h(x)可表示為：

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

(1)

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

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

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

(2)

(3)

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

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

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

(4)

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

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

(5)

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

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

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

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

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

(6)

(7)

(8)

AIC=-2lgL+2p。

(9)

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

3 結(jié)果與分析

3.1 模型精度評價

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

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

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

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

表4 預(yù)測模型交叉驗證統(tǒng)計量

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

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

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

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

研究區(qū)華北落葉松人工林SI預(yù)測值與主要環(huán)境因子的相關(guān)分析(表6)表明，在大的空間尺度上，SI與DEM、BIO14、TN和TP呈顯著正相關(guān)關(guān)系(P< 0.01)； 偏相關(guān)分析表明，控制其他環(huán)境因子不變，各立地生產(chǎn)力等級內(nèi)SI均對BIO14的變化響應(yīng)較弱； 第Ⅰ、Ⅴ等級SI主要受限于TN，第Ⅱ、Ⅲ等級SI與DEM的變化均呈顯著正相關(guān)關(guān)系(P< 0.01)，而第Ⅱ、Ⅲ等級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含量過高/低均會限制華北落葉松人工林生長。

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

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

4 討論

本研究探討研究區(qū)SI隨環(huán)境因子的變化規(guī)律，進(jìn)一步驗證了重要環(huán)境因子作為輔助變量預(yù)測立地生產(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)，在空間上也具有較強空間相關(guān)性(Moran’I=0.54，P< 0.01)。為進(jìn)一步探討其對立地指數(shù)的影響，本研究保留該解釋變量，剔除多重共線性問題嚴(yán)重(VIF > 15)的生長季溫度和降水等氣候因子。

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

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

5 結(jié)論

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