王冬至，劉紅艷，張冬燕，b，宋曉彤，徐安陽，李 靜

（河北農(nóng)業(yè)大學(xué) a.林學(xué)院；b.經(jīng)濟(jì)管理學(xué)院，河北 保定 071000）

樹高和胸徑是森林資源調(diào)查工作中最重要的兩個(gè)觀測(cè)因子，不僅是構(gòu)建森林生長與產(chǎn)量模型的基礎(chǔ)變量[1-2]，而且在預(yù)測(cè)林分結(jié)構(gòu)[3]、材積分布以及蓄積量[4-5]等研究中具有重要意義。Vargas等[6]研究表明，相比樹高而言胸徑的觀測(cè)不僅快速、準(zhǔn)確，并且觀測(cè)成本較低，而樹高的觀測(cè)則耗時(shí)耗力且費(fèi)用較高[7]。因此，在森林資源調(diào)查工作中，通常對(duì)樣地內(nèi)所有立木胸徑進(jìn)行實(shí)測(cè)，通過樹高與胸徑關(guān)系模型來預(yù)測(cè)樣地內(nèi)樹高，從而在滿足預(yù)測(cè)精度前提下可降低調(diào)查成本[8-9]。

目前，國內(nèi)外許多學(xué)者[10-15]大都基于不同形式線性和非線性模型來研究不同林分類型樹高與胸徑關(guān)系預(yù)測(cè)模型，其中非線性模型能夠更好地描述樹高與胸徑的關(guān)系[13，15]。最初樹高與胸徑關(guān)系模型中僅包含胸徑一個(gè)自變量[16-17]，隨著林齡的增加，樹高與胸徑關(guān)系通常會(huì)受到多種林分因子及立地因子的影響[18]。為了提高模型預(yù)測(cè)精度及適用性，Sánchez 等學(xué)者[19-21]在研究中加入了立地因子、林分因子及環(huán)境因子，構(gòu)建了包含不同變量及變量組合的樹高與胸徑關(guān)系模型。然而這些模型大都基于最小二乘法來估計(jì)模型參數(shù)[10,22]，在研究中由于調(diào)查數(shù)據(jù)經(jīng)常會(huì)存在觀測(cè)誤差、時(shí)間相關(guān)性和空間異質(zhì)性[3,23-24]，導(dǎo)致模型在應(yīng)用過程中產(chǎn)生較大預(yù)測(cè)誤差，而非線性混合效應(yīng)不僅能有效解決這些問題，而且在應(yīng)用過程中具有更好的靈活性[22-23]。Sharma 等學(xué)者[15,25-26]研究杉木純林和油松純林樹高與胸徑關(guān)系時(shí)，表明固定效應(yīng)模型預(yù)測(cè)及擬合精度均低于混合效應(yīng)模型。當(dāng)前，大多數(shù)學(xué)者均基于非線性混合效應(yīng)模型來構(gòu)建純林樹高與胸徑關(guān)系模型，然而針對(duì)多樹種、結(jié)構(gòu)復(fù)雜的針闊混交林，如何基于主要影響因子來構(gòu)建不同樹種樹高與胸徑關(guān)系的模型研究較少。因此，基于非線性混合效應(yīng)模型構(gòu)建針闊混交林不同樹種樹高與胸徑關(guān)系模型具有重要意義。

在研究多樹種樹高與胸徑關(guān)系模型中，通常采用啞變量[20]來解決樹種效應(yīng)對(duì)模型預(yù)測(cè)精度影響，Eerik?inen 等[27]研究表明包含啞變量預(yù)測(cè)模型比分樹種建立模型具有更好的預(yù)測(cè)效果。因此，本研究以河北省塞罕壩華北落葉松Larix principisrupprechtii與白樺Betula platyphylla針闊混交林為研究對(duì)象，首先基于候選樹高與胸徑關(guān)系模型，確定能夠描述樹高與胸徑關(guān)系的最優(yōu)基礎(chǔ)模型；其次通過相關(guān)性分析確定影響樹高生長主要因子；最后基于最優(yōu)候選模型和主要因子，構(gòu)建包含啞變量的多樹種非線性樹高與胸徑關(guān)系混合效應(yīng)模型，為多樹種混交林立地質(zhì)量評(píng)價(jià)及森林可持續(xù)經(jīng)營提供科學(xué)依據(jù)。

1 研究區(qū)概況

河北省塞罕壩機(jī)械林場(chǎng)（41°22′～42°58′N，116°53′～118°31′E）位于河北省最北部，地勢(shì)北高南低，處于陰山山脈與大興安嶺余脈的交接地帶，屬華北暖溫帶立地類型區(qū)，林場(chǎng)總面積約為9.2×104hm2，總蓄積約為8.1×106m3。研究區(qū)海拔在1 010～1 940 m 范圍內(nèi)，年均氣溫-1.2 ℃，年均最高氣溫33.4 ℃，年均最低氣溫-43.3 ℃；年均降水量約452.2 mm，年均蒸發(fā)量1 339.2 mm；土壤類型主要以褐色森林土、棕色森林土、砂壤土、壤土、風(fēng)沙土為主。研究區(qū)主要喬木樹種包括華北落葉松、白樺、蒙古櫟Quercus mongolica、樟子松Pinus sylvestris、云杉Picea asperata、山楊Populus davidiana、棘皮樺Betual davurica等，主要灌木樹種有胡枝子Lespedeza bicolor、繡線菊Spiraea salicifolia、稠李Prunus padus、忍冬Lonicera japonica等，草本植物主要有曼陀羅Datura stramonium、地榆Sanguisorba officinalis、唐松草Thalictrum aquilegifolium、歪頭菜Vicia unijuga、鐵線蓮Clematis florida、紫斑風(fēng)鈴草Campanula punctata、風(fēng)毛菊Saussurea japonica等。

2 研究方法

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

研究數(shù)據(jù)源于2015—2018年華北落葉松與白樺針闊混交林（華北落葉松為人工栽植、白樺為天然更新）標(biāo)準(zhǔn)地調(diào)查數(shù)據(jù)，在河北省塞罕壩機(jī)械林場(chǎng)下屬的陰河、第三鄉(xiāng)、北曼店和大喚起四個(gè)分林場(chǎng)，共設(shè)置華北落葉松與白樺針闊混交林標(biāo)準(zhǔn)地83 塊（30 m×30 m）。在樣地內(nèi)對(duì)海拔、土壤質(zhì)地、土層厚度、坡位、坡度、坡向等立地因子進(jìn)行觀測(cè)，并對(duì)樣地內(nèi)胸徑≥5 cm 立木進(jìn)行每木檢尺，共調(diào)查立木3 586 株（白樺1 698 株，華北落葉松1 888 株）。研究過程中，分樹種將觀測(cè)數(shù)據(jù)分別按3∶1 比例分為建模數(shù)據(jù)（62 塊標(biāo)準(zhǔn)地）和檢驗(yàn)數(shù)據(jù)（21 塊標(biāo)準(zhǔn)地），數(shù)據(jù)基本信息如表1、表2所示。

2.2 基礎(chǔ)模型選擇

當(dāng)前有近百種適用于不同林分類型的樹高與胸徑關(guān)系預(yù)測(cè)模型表達(dá)式[28]，本研究根據(jù)以往學(xué)者[16，28-31]的研究結(jié)果，選擇了5 個(gè)應(yīng)用范圍較廣且具有生物學(xué)意義的候選樹高-胸徑模型（表3），作為研究華北落葉松與白樺針闊混交林樹高與胸徑關(guān)系模型的基礎(chǔ)模型。

表2 模型檢驗(yàn)數(shù)據(jù)Table 2 Data of test model

表3 樹高與胸徑基礎(chǔ)模型Table 3 Basic model of height-diameter function

2.3 非線性混合效應(yīng)模型

樹木在生長過程中，樹高與胸徑關(guān)系通常會(huì)受到多種因子影響[10,32]，研究中采用非線性混合效應(yīng)建模技術(shù)[30-35]來構(gòu)建混交林中不同樹種樹高與胸徑關(guān)系模型，混合效應(yīng)模型一般形式可表示為：

式（1）中：yijk為第i個(gè)樣地中第j個(gè)樹種的第k株立木樹高值，M是樣地水平分組數(shù)量，Mi是樣地水平內(nèi)樹種分組數(shù)量，nij為第i個(gè)樣地第j個(gè)樹種株數(shù)，f為包含參數(shù)向量ψijk和變量xijk的非線性胸徑-樹高曲線模型，eijk為誤差項(xiàng)，β為固定效應(yīng)參數(shù)向量，Aijk、Bi,jk、Bijk分別為β、bi、bij的設(shè)計(jì)矩陣，μi、μij分別為第i個(gè)樣地中第j個(gè)樹種隨機(jī)效應(yīng)，ψ1、ψ2分別為μi、μij的方差-協(xié)方差矩陣且為正態(tài)分布，Rij為第i個(gè)樣地中第j個(gè)樹種的方差-協(xié)方差矩陣。

非線性混合效應(yīng)模型中內(nèi)方差-協(xié)方差殘差效應(yīng)由包括相關(guān)因子和加權(quán)因子的矩陣Rij來平衡誤差，其誤差項(xiàng)方差-協(xié)方差矩陣結(jié)構(gòu)如公式（2）所示。

式（2）中：Gij為描述異構(gòu)方差的nij×nij對(duì)角陣；Gij為誤差自相關(guān)矩陣nij×nij；s2為誤差擴(kuò)散的比例因子即模型剩余方差值[35]。

本研究中構(gòu)建樹高與胸徑非線性混合效應(yīng)模型，采用廣義正定矩陣表示隨機(jī)效應(yīng)的方差-協(xié)方差矩陣（D），樣地及混角比方差-協(xié)方差隨機(jī)效應(yīng)矩陣結(jié)構(gòu)如公式（3）所示。

式（3）中：為樣地間隨機(jī)效應(yīng)參數(shù)方差；為不同混交比隨機(jī)效應(yīng)參數(shù)方差；s1s2和s2s1為隨機(jī)效應(yīng)參數(shù)協(xié)方差（s1s2=s2s1）。

2.4 模型檢驗(yàn)與評(píng)價(jià)

研究中基于模型確定系數(shù)（R2）、絕對(duì)誤差（Bias）、均方根誤差（RMSE）、赤池信息準(zhǔn)則（AIC）及貝葉斯信息準(zhǔn)則（BIC）對(duì)模型擬合結(jié)果進(jìn)行檢驗(yàn)與評(píng)價(jià)。R2越大，絕對(duì)誤差、均方根誤差、赤池信息準(zhǔn)則及貝葉斯信息準(zhǔn)則越小表明模型擬合精度越高。

式（4）～（8）中：n為樣地中第r個(gè)樹種的第i株樹；n為樣地?cái)?shù)；hi為樹高觀測(cè)值（m）；為樹高預(yù)測(cè)值（m）；為林分平均樹高（m）。

3 結(jié)果與分析

3.1 基礎(chǔ)模型擬合與評(píng)價(jià)

在華北落葉松與白樺針闊混交林中，華北落葉松是優(yōu)勢(shì)樹種且為經(jīng)營目的樹種，因此研究中根據(jù)華北落葉松擬合結(jié)果來確定構(gòu)建華北落葉松與白樺混交林基礎(chǔ)模型。候選模型擬合結(jié)果及評(píng)價(jià)指標(biāo)如表4所示，其中Richard (M5)模型擬合效果較好，其確定系數(shù)（R2=0.918 6）最高，均方根誤差（RMSE=2.405 8）以及絕對(duì)誤差（Bias=0.194 5）最小，因此選擇Richard 模型作為本研究基礎(chǔ)模型。

表4 基礎(chǔ)模型擬合與評(píng)價(jià)Table 4 Basic model fitting and evaluation

3.2 相關(guān)性分析

主要林分因子（密度、優(yōu)勢(shì)種斷面積、林分總斷面積）與立地因子（海拔、坡度、土層厚度）和樹高Pearson 相關(guān)性分析如表5所示，其中海拔和林分總斷面積與樹高生長在P＜0.05 水平上存在顯著相關(guān)性，因此將海拔和林分總斷面積作為主要影響因子來構(gòu)建華北落葉松與白樺樹高曲線，從而提高模型預(yù)測(cè)精度及適用性。

表5 相關(guān)性分析?Table 5 Correlation analysis of different factors

3.3 非線性混合效應(yīng)模型構(gòu)建與評(píng)價(jià)

利用擬合精度較高的樹高與胸徑關(guān)系模型（Richard），基于逐步回歸分析和方差膨脹因子檢驗(yàn)，最終確定林分總斷面積及海拔是影響樹高主要因子，因此將這兩個(gè)因子加入基礎(chǔ)模型參數(shù)a中，生長速率參數(shù)b與形狀參數(shù)c受林分總斷面積及海拔影響較小，那么基于啞變量構(gòu)建的包含林分總斷面積及海拔的針闊混交林樹高與胸徑關(guān)系模型可以表示為：

式（9）～（11）中：a0、a1、a2為模型參數(shù)；bi和ci為啞變量，當(dāng)b0=0，b1=1，c0=0，c1=1 為白樺，當(dāng)b0=1，b1=0，c0=1，c1=0 為華北落葉松；BAL為林分?jǐn)嗝娣e（m2/hm-2）；ASL 為海拔（m）；μ1、μ2是隨機(jī)效應(yīng)參數(shù)。

基于主要因子構(gòu)建包含啞變量的樹高與胸徑關(guān)系非線性混合效應(yīng)模型擬合結(jié)果及評(píng)價(jià)指標(biāo)如表6～7 所示。隨機(jī)效應(yīng)樹高與胸徑關(guān)系模型擬合精度高于固定效應(yīng)模型，隨機(jī)效應(yīng)作用于林分?jǐn)嗝娣e擬合精度高于隨機(jī)效應(yīng)作用于海拔，其模型確定系數(shù)、均方根誤差和絕對(duì)誤差分別為：0.948 6、2.165 7、0.062 6。基于最優(yōu)模型對(duì)不同樹種預(yù)測(cè)殘差進(jìn)行了分析（圖1），各樹種殘差均不存在異方差現(xiàn)象，表明模型具有較好預(yù)測(cè)效果。研究中，當(dāng)隨機(jī)效應(yīng)同時(shí)作用在海拔和林分?jǐn)嗝娣e上時(shí)，模型不收斂。

表6 非線性混合效應(yīng)模型參數(shù)估計(jì)Table 6 Parameters estimation of nonlinear mixed effects model

表7 非線性混合效應(yīng)模型評(píng)價(jià)Table 7 Test of nonlinear mixed effects model

圖1 不同樹種的殘差分布Fig.1 The residual distribution of different species

4 結(jié)論與討論

本研究以方法研究為主，確定Richards 為構(gòu)建華北落葉松與白樺針闊混交林基礎(chǔ)模型，隨著林齡的增加，樹高與胸徑關(guān)系受多種林分因子影響，通過Pearson 相關(guān)性分析確定了海拔、林分總斷面積是影響華北落葉松與白樺混交林樹高與胸徑關(guān)系的主要因子；基于參數(shù)主要影響因子及最優(yōu)基礎(chǔ)模型，構(gòu)建了包含啞變量的樹高與胸徑關(guān)系混合效應(yīng)預(yù)測(cè)模型，當(dāng)隨機(jī)效應(yīng)作用在林分?jǐn)嗝娣e上時(shí)，包含啞變量的樹高與胸徑非線性混合效應(yīng)模型具有更高的預(yù)測(cè)精度。

為了研究針闊混交林中多樹種樹高與胸徑關(guān)系預(yù)測(cè)模型，本研究基于5 種具有生物學(xué)意義的基礎(chǔ)模型進(jìn)行擬合與評(píng)價(jià)，結(jié)果表明Richards 模型擬合精度最高，并與部分學(xué)者研究結(jié)果相同[29,34-36]。然而，有學(xué)者[10,37]認(rèn)為樹木生長會(huì)受到林分特征、樹木競爭關(guān)系、立地條件及經(jīng)營管理水平等[29,31]多種因素制約，Marziliano[17]研究表明林齡和密度是影響樹高與胸徑關(guān)系模型主要因子，并構(gòu)建了包含林齡和密度的樹高與胸徑關(guān)系預(yù)測(cè)模型；而Crecente-Campo 等學(xué)者[38-40]則用含有優(yōu)勢(shì)木胸徑、優(yōu)勢(shì)高以及林分密度等變量的樹高與胸徑關(guān)系模型來描述樹木樹高與胸徑的關(guān)系。本研究通過相關(guān)性分析確定林分?jǐn)嗝娣e和海拔是影響針闊混交林樹高與胸徑關(guān)系模型的主要因子，構(gòu)建了包含林分?jǐn)嗝娣e和海拔因子的樹高與胸徑關(guān)系模型，其中海拔會(huì)影響溫度從而進(jìn)一步對(duì)樹高和胸徑生長產(chǎn)生影響[1,7,36]，而林分?jǐn)嗝娣e在一定程度上反映了林分結(jié)構(gòu)對(duì)樹木樹高和胸徑生長產(chǎn)生影響[28,30]。

本研究將林分?jǐn)嗝娣e和海拔作為組合變量，構(gòu)建了包含啞變量的多樹種樹高與胸徑關(guān)系混合效應(yīng)模型?；趩∽兞繕?gòu)建的預(yù)測(cè)模型其擬合精度高于分樹種構(gòu)建的預(yù)測(cè)模型[20]，包含隨機(jī)效應(yīng)的預(yù)測(cè)模型精度高于未包含隨機(jī)效應(yīng)預(yù)測(cè)模型，其他部分學(xué)者[30-42]通過研究不同林分類型樹高與胸徑關(guān)系模型，均表明混合效應(yīng)模型在擬合精度和預(yù)測(cè)精度上均高于基礎(chǔ)模型，Huang 等[40]認(rèn)為隨機(jī)效應(yīng)能夠有效解決樣地及單木水平間差異性。樹木在生長過程中，不同林分類型其樹高與胸徑關(guān)系會(huì)受到多種因子影響，有學(xué)者[27,35-38]基于林分密度、優(yōu)勢(shì)高、林分?jǐn)嗝娣e及空間位置等影響變量，構(gòu)建了不同變量或變量組合的非線性混合效應(yīng)模型，均在不同程度上提高了模型預(yù)測(cè)精度。當(dāng)隨機(jī)效應(yīng)參數(shù)作用于林分?jǐn)嗝娣e時(shí)的預(yù)測(cè)精度高于隨機(jī)效應(yīng)參數(shù)作用于海拔，當(dāng)構(gòu)建林分?jǐn)嗝娣e和海拔兩水平非線性混合效應(yīng)模型時(shí)，不同模型采用了同一參數(shù)初始值和收斂準(zhǔn)則，導(dǎo)致模型參數(shù)不收斂。在今后構(gòu)建多樹種的樹高與胸徑關(guān)系模型研究中，多水平非線性混合效應(yīng)建模技術(shù)與方法還需進(jìn)一步深入研究和探討。