黃金金，劉曉彤，張逸如，李?？?/p>

(1. 中國林業(yè)科學研究院資源信息研究所，北京 100091；2. 國家林業(yè)和草原局森林經(jīng)營與生長模擬重點實驗室，北京 100091)

森林蓄積量和生物量是反映森林資源水平的基本指標，精準估測森林蓄積量和生物量，可以為林業(yè)資源管理決策提供科學依據(jù)[1]。國內(nèi)外關(guān)于蓄積量和生物量模型的研建主要集中于單木水平，而林分水平的較少，尚未形成我國林業(yè)行業(yè)標準[2]。林分水平的生長模型能直接預測林分蓄積量和生物量，減少單木的誤差積累[3-4]。目前估計區(qū)域尺度蓄積量和生物量的建模方法有多種，應(yīng)用較多的一是以異速生長方程為基礎(chǔ)，根據(jù)林分總胸高斷面積、林分優(yōu)勢高、林分平均高等估計林分以及區(qū)域的地上生物量[5-6]；二是生物量轉(zhuǎn)換因子法，張會儒等[7]將全國森林類型按照林分優(yōu)勢樹種歸并為9 類，構(gòu)建了全國主要樹種組林木與材積相容的生物量回歸模型；三是通過單木模型推算，傅煜等[8]通過Monte Carlo 法反復模擬由單木生物量模型推算區(qū)域尺度地上生物量，估計了江西省杉木地上總生物量，秦立厚等[9]分析量化單木生物量模型中各種不確定性來源，估計區(qū)域地上生物量。

大量研究表明，林木胸徑與林木年齡之間存在緊密的生長關(guān)系，為建立有關(guān)的胸徑生長模型奠定了基礎(chǔ)[10-14]。曹磊[15-16]和龍時勝等[17]采用Richards生長方程擬合樹木直徑與年齡的生長關(guān)系，模型的擬合優(yōu)度較高，其中龍時勝等[17]利用模型估計林分平均年齡時，擬合年齡與真實年齡的絕對誤差最大值為2 年，相對誤差最大值為8.2%，兩者間的差異性較小，說明利用林木多期直徑估計異齡林林分平均年齡具有較高的可行性。大量研究[18-20]認為，林分蓄積量、生物量生長主要受樹木自身大小、競爭因子和立地條件的影響。

本研究以廣東省3 個針葉樹種155 塊樣地為對象，利用5 期連清數(shù)據(jù)，基于林分平均胸徑與年齡的生長關(guān)系估計樣地林齡，進而建立林分水平的蓄積量和生物量生長模型，以體現(xiàn)生長潛力的參數(shù)a分級和與生長速度有關(guān)的參數(shù)b分級、考慮林分競爭與否、分步建模（一元非線性回歸）和聯(lián)合建模（非線性誤差變量聯(lián)立方程組）共8 種組合，構(gòu)建模型系，評價模型的擬合優(yōu)度和預測效果，篩選最優(yōu)模型，以期實現(xiàn)對廣東省針葉林蓄積量和生物量的準確估計。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

廣東省位于20°13′～25°31′ N 和109°39′～117°19′ E，處于東亞季風區(qū)，包括了熱帶、中亞和南亞熱帶3 種氣候，年均氣溫22.3℃，月降水量約313 mm，年日照時數(shù)約1 746 h，太陽輻射總量約4 200～5 400 MJ·m?2，生長期水熱充足，利于植被生長。地貌類型復雜多樣，以低山、丘陵為主，地勢北高南低。

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

建模數(shù)據(jù)：從廣東省第5—9 次連清數(shù)據(jù)中，選擇馬尾松（Pinus massonianaLamb.）、濕地松（Pinus elliottiiEngelm.）、杉木（Cunninghamia lanceolata(Lamb.) Hook.）3 個針葉樹種，篩選5 次清查中皆存在且保留木大于5 株的樣地，共155 塊（馬尾松林91 塊、濕地松林16 塊、杉木林48 塊），樣地分布見圖1。固定樣地為0.066 7 hm2，樣木調(diào)查因子包括樣木號、胸徑，單株材積按廣東省地方一元材積表[21]計算獲得，單株生物量按行業(yè)標準生物量模型[22-24]獲得。為避免株數(shù)差異造成樣地蓄積量和生物量的差異，本研究中樣地蓄積量和生物量為樣地內(nèi)5 期保留木的蓄積量和生物量算數(shù)平均值，統(tǒng)計量如表1 所示。

表1 建模數(shù)據(jù)統(tǒng)計量Table 1 Modeling data statistics

圖1 建模數(shù)據(jù)固定樣地分布示意圖Fig. 1 Distribution map of fixed sample plots of modeling data

樣地胸徑采用林分平均胸徑：

式中，Dg為林分平均胸徑，n為樣木株數(shù)，di為樣木胸徑。

驗證數(shù)據(jù)：選取2002 年達到起測胸徑、存在于4 次連清中，且保留木大于5 株的樣地，共126塊，其中馬尾松59 塊、濕地松22 塊、杉木45 塊。各統(tǒng)計量如表2 所示。

表2 驗證數(shù)據(jù)統(tǒng)計量Table 2 Verification data statistics

2 方法

2.1 模型構(gòu)建

采用Mitscherlich 理論生長方程，擬合5 期林分平均胸徑和相應(yīng)時間間隔之間的關(guān)系；以Richards理論生長方程，擬合樣地平均蓄積量和生物量與相應(yīng)時間間隔之間的關(guān)系。

式中：Dgij、Vij、BIOij分別為第i個樣地第j期的樣地平均胸徑、蓄積量和生物量，agei是第i個樣地的初期林齡，adj為第j期與初期（1997 年）的時間間隔，a1、a2、a3、b1、b2、b3、c2、c3為參數(shù)。

基于上述胸徑、蓄積量和生物量3 個生長方程，分別以體現(xiàn)生長潛力的參數(shù)a分級和與生長速度有關(guān)的參數(shù)b分級、考慮競爭與否、分步或聯(lián)合建模等8 種組合，構(gòu)建模型系（表3）。

表3 8 種模型系Table 3 Eight model systems

2.1.1 含林木競爭的胸徑生長方程 對于林分水平模型來說, 與樹冠有關(guān)的競爭指數(shù)能較好地反映樹木間的競爭關(guān)系，因此選用基于交角且對象木與競爭木距離小于10 m 的競爭指數(shù)ju10，詳見惠剛盈等[25]。由于林木競爭對單木胸徑的生長具有較大的影響，但對林分平均胸徑的影響并不算太大，為解釋林木競爭對林分生長的影響，故將競爭指數(shù)引入了胸徑生長方程：

式中：b0為競爭指數(shù)參數(shù)，ju10ij為第i個樣地第j期的競爭指數(shù)（樣地競爭指數(shù)為樣地內(nèi)所有保留木競爭指數(shù)的算數(shù)平均值），其余同上。

2.1.2 參數(shù)分級 為了兼顧區(qū)域尺度上立地質(zhì)量的多樣性以提高模型在大尺度范圍內(nèi)的適用性，對參數(shù)a或b分級代替立地質(zhì)量的差異，以參數(shù)a分級為例：

將式（2）作為立地質(zhì)量劃分等級的分級方程，已知林分平均胸徑Dgij和時間間隔adj，以樣地號和樹種號構(gòu)成的唯一編號作為參數(shù)a和age的啞變量，樹種作為參數(shù)b的啞變量，擬合方程（6）。將參數(shù)a1i從低到高排序，并按照有序樣本聚類分為6 類。

式中：a1i(i=1,2,...,155)是第i個樣地的參數(shù)，b1k(k=1,2,3時分別代表馬尾松、濕地松和杉木)是第k個樹種的參數(shù)，其余同上。

2.1.3 分步建模 以Model 5 為例，基于式（2），以參數(shù)a區(qū)別樹種，以參數(shù)b區(qū)別立地質(zhì)量差異，以唯一編號作為初期林齡age（未知）的啞變量，根據(jù)5 期林分平均胸徑Dgij和時間間隔adj之間的關(guān)系來估計樣地初期林齡agei，見式（7）。將agei作為式（3）和（4）的自變量，以參數(shù)c區(qū)別樹種，其余參數(shù)處理同上，分別擬合樣地平均蓄積量Vij和生物量BIOij與林齡agei的關(guān)系，即式（8）和（9）。

式（7）～（9）中：a1k、a2k、a3k(k=1,2,3)分別為胸徑、蓄積量和生物量生長方程中的第k個樹種的參數(shù)，b1m、b2m、b3m(m=1,2,…,6)分別為胸徑、蓄積量和生物量生長方程中的第m個立地質(zhì)量等級的參數(shù)，c2k、c3k(k=1,2,3)分別為蓄積量和生物量生長方程中的第k個樹種的參數(shù)，其余同上。

2.1.4 聯(lián)合建模 為了比較建模方法對生長模型擬合優(yōu)度和估測精度的差異，除分步建模法外，本研究還采用了聯(lián)合建模法，以Model 7 為例。由于林齡未知，且3 個方程中年齡參數(shù)agei是一致的，基于式（2）～（4），參數(shù)處理同分步建模法，根據(jù)林分平均胸徑、蓄積量和生物量與時間間隔之間的關(guān)系，采用非線性誤差變量聯(lián)立方程組法同時擬合3 個方程，即式（10）。

2.2 模型評價

考慮各模型參數(shù)不同，為了檢驗8 種模型系的擬合效果，采用以下4 項模型評價指標：確定系數(shù)R2、估計值標準誤差SEE、平均預估誤差MPE、總體相對誤差TRE，公式詳見文獻[26]。

2.3 模型驗證

（1）確定立地質(zhì)量等級和林齡（以參數(shù)b 分級為例）

根據(jù)胸徑生長方程反推林齡計算公式，如下：

式中，pageij為第i個樣地第j期的樣地林齡，當不引入競爭指數(shù)時b0=0，其余同上。

由于驗證樣本所在樣地的立地質(zhì)量等級m和初期林齡皆未知，需先確定立地質(zhì)量等級，再由式（11）求出林齡，才可估計蓄積量和生物量。當立地質(zhì)量分級時，至少需要兩期胸徑才能確定樣地的立地質(zhì)量等級，薛春泉等[27]發(fā)現(xiàn)長周期確定的等級比短周期穩(wěn)定，因此本研究采用期初（2002年）和期末（2017 年）2 期胸徑來確定樣地的立地質(zhì)量等級，確定方法如下：

式中，pagei1和pagei4為第i個樣地期初和期末的樣地林齡，該方法的直觀解釋就是使推算的2002 年和2017 年林齡間隔與實際林齡間隔（15 年）相差最小的立地質(zhì)量等級，即為該樣地的立地質(zhì)量等級。

根據(jù)式（12）確定樣地的立地質(zhì)量等級后可知參數(shù)bm，再由式（11）計算初期（2002 年）林齡age，則2007 年、2012 年、2017 年時樣地的林齡分別由初期林齡相應(yīng)加上5、10、15 年。

（2）模型估測效果評價指標

由蓄積量和生物量生長方程估計樣地平均蓄積量和生物量，樣地總蓄積量和總生物量可由平均蓄積量和生物量乘以株數(shù)獲得，按式（13）～（14）計算各樹種4 個時期下區(qū)域尺度蓄積量和生物量估計的總相對誤差EV和EBIO。

式中：n為樣地數(shù)，Ni為第i個樣地的保留木株數(shù)，PVi、Vi分別為第i個樣地平均蓄積量的估計值和真實值，PBIOi、BIOi分別為第i個樣地平均生物量的估計值和真實值。

2.4 生長量函數(shù)

連年生長量函數(shù)Z(t)為總生長過程曲線y(t)的一階導數(shù)，為樹木在某年的實際生長速度，即連年生長量隨年齡t變化的函數(shù)。Richards 生長方程的連年生長量方程為：

平均生長量函數(shù)θ(t)用來表示樹木在某一時刻的平均生長速度，總生長過程曲線y(t)除以年齡t則可得平均生長量依年齡變化的函數(shù)。Richards生長方程的平均生長量方程為：

3 結(jié)果

3.1 模型擬合

表4 列出了8 種模型系的擬合評價指標。整體上看，8 種模型系的3 個生長方程擬合確定系數(shù)R2皆達0.980 以上，總相對誤差TRE在±0.5%以內(nèi)，平均預估誤差MPE在±1%以內(nèi)，擬合效果良好，但各模型間的4 項指標差異不大。當控制建模方法和立地質(zhì)量分級的參數(shù)相同、僅比較競爭影響時，引入競爭指數(shù)對模型擬合優(yōu)度的提升沒有明顯的作用，如Model 5 和Model 6 擬合指標間差異甚微、Model 2 較Model 1 的蓄積量和生物量生長方程的擬合效果甚至出現(xiàn)了一定程度的下降；類似地，可以發(fā)現(xiàn)參數(shù)b分級模型的擬合優(yōu)度較參數(shù)a分級模型稍好；當控制競爭影響和分級參數(shù)相同時，相對于分步建模，采用含度量誤差的聯(lián)立方程組時，蓄積量和生物量生長方程的擬合效果有所提高，但是胸徑生長方程的擬合效果稍差。

表4 模型擬合評價指標Table 4 Model fitting evaluation index

3.2 模型驗證

8 種蓄積量和生物量生長方程的模型驗證指標分別如圖2 和圖3 所示?？傮w而言，各模型的預測效果都較好，從圖2 可知，區(qū)域尺度蓄積量估計誤差最大的僅有17.62%，說明蓄積量估計精度可達82.38%以上；從圖3 可知，區(qū)域尺度生物量估計誤差最大的僅有12.51%，則生物量估計精度可達87.49%以上。8 個模型中以基于參數(shù)b分級的分步模型（即Model 5 和Model 6）預測效果最優(yōu)，其對3 個樹種4 期蓄積量和生物量的估計誤差較其余模型稍小，且在不同樹種和不同時期的表現(xiàn)都較為穩(wěn)定，說明模型適用性較好。以Model 5 為例，3 個樹種中對杉木的估計效果稍好，估計誤差最大為?5.25%，而估計馬尾松和濕地松時誤差最大為10.00%和10.36%。8 個模型對3 個針葉樹種的4 期區(qū)域尺度蓄積量和生物量估計的誤差，大都表現(xiàn)出中期誤差高于首末2 期，初期估計精度較高，但明顯出現(xiàn)低估的現(xiàn)象，伴隨著樣地林齡的增大，中期出現(xiàn)高估且精度稍低，但到末期（2017 年）時估計誤差降低，精度有所提高。

圖2 8 種蓄積量生長模型的區(qū)域尺度估計誤差Fig. 2 Regional-scale estimation errors of eight volume growth models

圖3 8 種生物量生長模型的區(qū)域尺度估計誤差Fig. 3 Regional-scale estimation errors of eight biomass growth models

3.3 最優(yōu)模型

綜合8 個模型的擬合效果（表4）和驗證效果（圖2、圖3），Model 5 最優(yōu)，Model 6 次之。表5 和表6 列出了基于參數(shù)b分級的分步模型的模型參數(shù)估計值。由表5 可知，馬尾松胸徑、蓄積量和生物量的生長潛力均為最大，杉木最小，其中馬尾松蓄積量和生物量的生長極限值可達濕地松和杉木的2 倍以上。由表6 可知，立地質(zhì)量等級越高，樹木的生長速度越快，但是各等級下的林分蓄積量生長速度均為最快，其次為生物量，胸徑生長速度最慢。

表5 樹種啞變量參數(shù)和競爭指數(shù)參數(shù)Table 5 Dummy variable parameters of tree species and competition index parameters

表6 基于參數(shù)b分級的立地質(zhì)量等級啞變量參數(shù)Table 6 Site quality grade dummy variable parameters based on parameterbclassification

3.4 蓄積量和生物量生長量函數(shù)

基于參數(shù)b分級不考慮競爭的分步建模（即Model 5）最優(yōu)，為節(jié)約篇幅，選擇樣地立地質(zhì)量等級占比較多的等級2、3 和4，繪制Model 5 下蓄積量、生物量的連年生長量和平均生長量的關(guān)系曲線圖（圖4）。蓄積量和生物量的連年生長量、平均生長量均顯示出立地質(zhì)量越好，生長量能達到的極限值越大，且所需年限更短，各相鄰等級間的差異隨著等級的增加而逐漸加大。對于不同的樹種來說，同一立地質(zhì)量等級下的生長量均為馬尾松>濕地松>杉木，生長量達到最大時的林齡為杉木>濕地松>馬尾松，各樹種間生長量差異隨著等級的增加而逐漸加大。立地質(zhì)量等級越高，平均生長量和連年生長量曲線相交時的樣地年齡越小。

4 討論

處理多樹種問題是建立林分模型的一個難點，通常是將組成樹種進行分組或合并[28]。曾偉生等[2]建立了東北林區(qū)林分水平的三儲量模型，篩選出10 種主要森林類型并引入啞變量體現(xiàn)差異后，蓄積量模型、生物量模型和碳儲量模型的確定系數(shù)分別從0.945、0.805、0.839 提高到0.959、0.949、0.951。本研究對存在多個組成樹種的樣地按優(yōu)勢樹種分屬不同的樣地，并構(gòu)建啞變量區(qū)別廣東省的3 個主要針葉樹種，從表5 和圖4 可知，各樹種間的生長潛力和生長量有較大的差異，表現(xiàn)為馬尾松>濕地松>杉木，其中馬尾松遠超2 倍左右，這是因為建模數(shù)據(jù)（表1）中馬尾松保留木大多為天然林，而杉木的起源大多為人工林，Zeng 等[29]研究結(jié)果表明同一直徑的單木地上生物量估計值在不同起源上存在一定差異，天然林的生物量大于人工林。此外，在同一立地質(zhì)量下，不同樹種間生長量的差異主要是由于樹種特性的差異引起的，不同樹種對同一區(qū)域的生長適應(yīng)程度并不一致，馬尾松可能更適宜廣東的立地條件。各樹種間區(qū)域尺度蓄積量和生物量的預測精度也有較大的差異，對杉木的估計效果較馬尾松和濕地松稍好，可能原因是不同樹種的建模樣本量和檢驗樣本量不一致，另一個可能的原因是不同的樹種有不同最優(yōu)的基礎(chǔ)生長模型，如吳宏煒等[30]研究發(fā)現(xiàn)，Schumacher 模型是濕地松林分蓄積量模型的最優(yōu)基礎(chǔ)模型。對于估計誤差，大都表現(xiàn)出中期高于首末兩期的現(xiàn)象，一方面可能是由于采用首末兩期胸徑數(shù)據(jù)來確定驗證樣本所在樣地的立地質(zhì)量等級，相較于中期，推算的各樣地立地質(zhì)量等級更貼合期初和期末時的實際情況，因此用于林齡確定、估計蓄積量和生物量時，3 個生長方程中反映立地質(zhì)量差異的參數(shù)更為穩(wěn)定，模型的預測效果會更好。另一方面隨著林齡增大，蓄積量和生物量逐漸積累，公式（13）和（14）分母變大，這可能使后期的估計誤差降低，估計效果更好。

非線性獨立回歸估計法的估計效果優(yōu)于含度量誤差變量聯(lián)立方程組法，因為分步建模時，3 個方程都是獨立擬合給出參數(shù)的最優(yōu)估計，而采用聯(lián)立方程組法聯(lián)立求解時，由于方程中林齡一致的約束，為保證方程之間的聯(lián)系，得到更優(yōu)化的參數(shù)估計，各方程較獨立擬合時都有了不同的變化。從擬合評價指標（表4）可以看出，對比分步建模法，采用聯(lián)立方程組法可以提高蓄積量和生物量生長方程的擬合效果，但是卻降低了胸徑生長方程的擬合優(yōu)度，導致林分年齡的估計精度有所下降，使用驗證數(shù)據(jù)集檢驗?zāi)Ｐ偷倪m用性時，由于林分年齡未知，確定首期林齡成為關(guān)鍵的一步，所以胸徑生長方程的擬合效果直接影響了驗證時的估計精度。由于林分大部分為混交林，具有年齡不一、樹種結(jié)構(gòu)關(guān)系復雜等特點，因此探究林木間競爭作用對林分結(jié)構(gòu)和生長變化的影響至關(guān)重要。一般而言，模型中的參數(shù)個數(shù)越多，擬合效果會更優(yōu)，但是引入競爭指數(shù)后對于模型的擬合效果（表4）和估計精度（圖2、3）沒有明顯的提高。雖然從理論上講包含樹木空間位置的競爭指數(shù)對預測林木生長更有效, 但Biging 等[31]研究發(fā)現(xiàn)與距離無關(guān)的競爭指數(shù)比與距離有關(guān)的競爭指數(shù)預測單木生長更好，臧灝等[32]探究人工杉木林胸徑生長對競爭的響應(yīng)時也發(fā)現(xiàn)對比于與距離有關(guān)的單木競爭指標，采用與距離無關(guān)的單木競爭指標的生長模型估計精度更高。競爭指數(shù)表現(xiàn)的預測能力依賴于樹種、立地條件等因素，這些因素的變化會導致不同類型的競爭指數(shù)表現(xiàn)不一致，且不具有決定性[33-34]。本研究引入競爭指數(shù)對林分模型的適用性沒有明顯影響，可能與選取的競爭指數(shù)有關(guān)，也有可能和立地的限制有關(guān)，因為在較大范圍內(nèi)密度對林分蓄積量、生物量的影響不大[35]。

模型的應(yīng)用范圍和模型的預測精度相關(guān)聯(lián)，當模型應(yīng)用于大尺度的森林資源測算時，若忽略立地質(zhì)量的差異必然導致模型估計結(jié)果的不準確。立地質(zhì)量是指在某一立地上既定森林或其他植被類型的生產(chǎn)潛力[35]，既定立地上的不同樹種和相同的樹種在不同的立地上的立地質(zhì)量都是有差異的。本研究為體現(xiàn)廣東省樣地間的立地質(zhì)量多樣性，分別對與生長潛力有關(guān)的參數(shù)a和與生長速度有關(guān)的參數(shù)b進行分級，結(jié)果表明以參數(shù)b分級體現(xiàn)立地質(zhì)量的差異優(yōu)于參數(shù)a分級。薛春泉等[26]研究發(fā)現(xiàn)生長速度分級的胸徑生長方程與地上生物量的異速生長方程構(gòu)成的模型系，不僅擬合效果達到了0.95以上，區(qū)域尺度生物量的估計誤差也下降到了8%以內(nèi)，這與本研究參數(shù)b分級模型預測效果較好具有一致性。另外，林分蓄積量和生物量的生長量隨著立地質(zhì)量的提高而增加，說明立地質(zhì)量等級越高，森林蓄積量和生物量密度越大。立地質(zhì)量與森林生產(chǎn)力相關(guān)，不同的立地質(zhì)量代表著不同的生產(chǎn)力，立地質(zhì)量越高，林分的生長速率越快，達到最大生長量所需的年限也更短。趙菡等[36]用優(yōu)勢木樹高等級代替立地等級估計江西省不同立地質(zhì)量等級的馬尾松林生物量，研究發(fā)現(xiàn)不同立地等級的地上生物量均值估計結(jié)果隨著立地質(zhì)量等級的升高而增大，這與本研究結(jié)果一致。

5 結(jié)論

本研究以廣東省3 種針葉優(yōu)勢樹種為研究對象，利用1997—2017 年的5 期森林資源連續(xù)清查數(shù)據(jù)，基于林分平均胸徑、蓄積量和生物量與相應(yīng)調(diào)查時間間隔的關(guān)系建立林分生長模型，實現(xiàn)了對廣東省蓄積量和生物量的準確估計。本研究結(jié)論如下：（1）基于與生長速度有關(guān)的參數(shù)b分級不考慮競爭的分步模型最優(yōu)，引入競爭指數(shù)與否差異不大。（2）最優(yōu)模型的區(qū)域尺度蓄積量和生物量估計誤差最大分別為10.36%和10.22%，模型的適用性較好，4 期的估計誤差表現(xiàn)為中期高于首末兩期。（3）不同樹種中以馬尾松的生長潛力最大，杉木最小，但是模型對杉木的估計效果較馬尾松和濕地松稍好。（4）立地質(zhì)量等級越高，生長量極值越大，達到最高峰所需年限也更短；同一立地質(zhì)量下馬尾松的生長量最大，其次為濕地松，杉木最小。