年齡無關(guān)的生長模型研究—以落葉松平均高為例

國 紅，雷淵才，郎璞玫

(中國林業(yè)科學(xué)研究院資源信息研究所，北京 100091)

了解和掌握森林資源的動態(tài)變化，對制定森林撫育經(jīng)營方案和措施，促進森林質(zhì)量精準(zhǔn)提升具有重要的意義。森林生長模型是了解和掌握森林資源動態(tài)變化的一項重要工具。這類方程的形式很多，可以分為理論生長函數(shù)、經(jīng)驗生長函數(shù)和Box-Lucas 生長函數(shù)，而著名的理論生長函數(shù)有Richards模型，Schumacher 模型，Schnute 模型等[1-3]。這些模型一直被廣泛應(yīng)用，如模擬樹高生長，優(yōu)勢高生長，斷面積生長等。森林生長模型在構(gòu)建和應(yīng)用研究中，通常使用年齡作為重要因子，但林木年齡因子往往存在兩方面的問題，一是測定困難，一般需要用生長錐測定具有清晰年輪樹種的從髓心開始的整個年輪的數(shù)據(jù)，對于沒有清晰年輪的樹種，測定往往非常困難，而絕大多數(shù)闊葉樹不具有清晰的年輪；二是年齡不一定有效，在樹木年齡差異較大的林分中，年齡不能用于預(yù)測林木或林分生長和產(chǎn)量，特別是對于天然林的生長和動態(tài)監(jiān)測更加依賴于年齡無關(guān)的生長模型。因此，年齡無關(guān)的森林生長模型的研究變得越來越必要。

國內(nèi)外許多學(xué)者開展了一些與年齡無關(guān)的森林生長模型的相關(guān)研究。Martin 等[4]，Monserud 等[5]，Trasobares 等[6]通過使用與年齡無關(guān)的經(jīng)驗線性模型模擬了單木的生長，但是這類模型可能不能確保單木或林分生長有適合的生長曲線特性，即生長曲線達到最大生長量后生長逐漸下降。而應(yīng)用理論生長函數(shù)這個特性進行外推預(yù)測林木生長時，可以保證單木或林分的生物學(xué)特性。

Vanclay[7-8]考慮單木或林分生長曲線的生物學(xué)特性，構(gòu)建了具有這種性質(zhì)的線性經(jīng)驗函數(shù)，并已經(jīng)被許多作者成功地使用。這類生長函數(shù)的表達式如下：

其中Y 為林木變量生長量，dY 和dt 為林木變量生長率，k 可以取為1 或2，常數(shù)β0可以表達為環(huán)境變量和(或) 競爭變量的函數(shù)。然而，這些參數(shù)沒有特定生物學(xué)意義的含義。

Huang 和Titus[9]為了模擬阿爾伯達北部混交林中白云杉(Picea glauca (Moench) Voss) 的直徑增量生長變化，比較了前面生長函數(shù)的非線性形式與能描述單木或林分生長曲線形狀的非線性Box-Lucas 函數(shù)，結(jié)論是Box-Lucas 函數(shù)模擬直徑生長效果比非線性的經(jīng)驗生長函數(shù)好。Box-Lucas 函數(shù)表達式為：

其中：θ1和θ2是估計增量曲線形狀的參數(shù)，dY 和dt 為林木變量生長率。這些函數(shù)和其他類似的函數(shù)可以用來對單木或林分變量的增量進行建模。這些參數(shù)也不具有特定生物學(xué)意義的含義。

理論生長函數(shù)通常建立在公認的生物學(xué)關(guān)系的基礎(chǔ)上，因此它們的模型參數(shù)可以從生物學(xué)的角度進行解釋。不管樹種和林分類型，樹木的生長都遵循某種S 型或凸型生長曲線，它顯示出樹木壽命的3 個階段。模擬樹木和林分生長的一種非常成功的方法是使用用差分方程表示的生長函數(shù)，如Bailey 等[10]和Borders 等[11]提出的，已被許多學(xué)者使用[12-14]。差分方程表示一組生長函數(shù)，除了一個“自由” 參數(shù)外，所有參數(shù)都是公共的。但多數(shù)這類差分方法并沒有隱去年齡。Tomé等[15]利用隱去年齡的差分方法將Lundqvist 模型和Richards 模型分別應(yīng)用于葡萄牙桉樹優(yōu)勢高生長和軟木櫟胸徑生長的研究；Gea-lzquierdo 等[16]基于Hossfeld IV 和korf模型的差分形式模擬樹木直徑的生長。葛宏立[17]以Johnso -Schumacher 模型為基礎(chǔ)，導(dǎo)出年齡隱含的模型應(yīng)用于森林資源連續(xù)清查數(shù)據(jù)的生長研究，鄒奕巧等[18]用麗水市的馬尾松樣木胸徑數(shù)據(jù)，導(dǎo)出了Richards 模型隱去年齡的形式，并采用單點估計法和雙點估計法2 種方法進行研究。但是，此類研究目前為止相對較少，而且這些研究往往沒有或者較少考慮到環(huán)境因子對于生長的影響。本研究以理查德模型為基礎(chǔ)研究隱含年齡的差分形式，并建立差分方程的系數(shù)參數(shù)與環(huán)境因子和林分因子的相關(guān)函數(shù)，以便更好的解決年齡不可及或者不是有效變量的問題和更全面解釋單木和林分生長的生物學(xué)特性，為混交林質(zhì)量精準(zhǔn)提升提供一種有效的經(jīng)營工具。

1 材料和方法

1.1 材料

研究采用2004 和2009 年吉林省兩期177 塊落葉松固定樣地數(shù)據(jù)，測量間隔為5 a。固定樣地為20 m×30 m，主要調(diào)查因子有胸徑、林分年齡、平均樹高、郁閉度、地位級、坡位、坡向、坡度、海拔高度、土層厚度和土壤類型等（見表1）。落葉松樣地的選取原則是樹種組成中，落葉松的胸徑斷面積超過70%，平均高通過每塊樣地3 株平均木的樹高取平均得到，平均年齡通過生長錐鉆取3 株平均木的樹芯計算的年齡平均值得到。研究區(qū)氣候數(shù)據(jù)從WorldClim 數(shù)據(jù)庫下載得到（網(wǎng)址：www.worldclim. org）。每個樣地可以下載19 個用于構(gòu)建氣候敏感的生物量模型的候選氣候因子（表2），分辨率為1 km×1 km。這些因子可以分為兩大類，即溫度有關(guān)的變量10 個和降水量有關(guān)的變量9 個。每個氣候因子是1990—2000 年的平均值。

表1 林分和氣候因子 Table 1 Stand and climate factors

利用雙重篩選逐步回歸方法，建立兩期平均高和19 個氣候因子的線性關(guān)系，入選的因子為最濕季度的平均溫度BIO8、年降水量BIO12 和最濕季度的降水量BIO13，這3 個氣候因子在建模中最為敏感，本研究予以選用。具體含義見表2。

表2 入選氣候因子含義Table 2 Implications for selected climate factors

1.2 方法

理查德方程認為一個器官的生長可以表達成有機物合成和降解之差。假設(shè)合成可以被表達成質(zhì)量的異速生長關(guān)系，而降解和質(zhì)量是成比率的，因此生長可以被表達為：

式中：β1、β2和m是參數(shù)。生長方程是不基于年齡的，因此對于異齡林樣地或者年齡未知的異齡林樣地數(shù)據(jù)結(jié)構(gòu)模擬單木或林分生長是非常適用的。實際應(yīng)用中，上式瞬時生長量用單木或林分收獲量表達即形式可以由式（3）積分為：

式中：參數(shù)可以分別由下式表達：

A是當(dāng)t趨于無窮時的漸近線，k是和生長率相關(guān)的參數(shù)，m是形狀參數(shù)。c= e-kt0，t0是Y=0 時的年齡。大多數(shù)應(yīng)用中，如果t0=0 時C=1。這里取C=1.

與年齡無關(guān)的差分形式可以通過類似于通常用于將生長函數(shù)公式化為差分方程的方法(如Amaro等[12]的過程來獲得）?？紤]到吉林省177 塊天然落葉松林的連續(xù)清查數(shù)據(jù)缺乏年齡，為了建立單木或林分的樹高生長模型，本研究基于理查德方程（4）推導(dǎo)不含年齡的差分生長方程。根據(jù)式（4），第i年和第i+a年的平均高為：

在這里，a為預(yù)測期，或為連續(xù)清查時間間隔，Yi+a是預(yù)測期（或連續(xù)清查時間間隔）的單木或平均高，Yi是連續(xù)清查時間初期的單木或平均高，A是當(dāng)t趨于無窮時的漸近線，k是和生長率相關(guān)的參數(shù)，m是形狀參數(shù)。式（12）表達比較復(fù)雜，但是當(dāng)擬合單木樹高或平均樹高生長時參數(shù)很容易收斂。

為了模擬平均樹高生長，本研究分別嘗試采用k為固定參數(shù)和自由參數(shù)兩種形式。也考慮了m為自由參數(shù)的情況，但獲得的結(jié)果是m受林分因子、氣候因子的影響不明顯。當(dāng)k為自由參數(shù)時，將k表達為林分因子和氣候因子的線性函數(shù)關(guān)系。

其表達形式分別有如下4 個形式，并進行比較。

在這里，b0～b5是估計參數(shù)，BA是樣地斷面積，S是土壤厚度，BIO8 為最濕季度的平均溫度，BIO12 是年降水量，BIO13 是最濕月份的降水量。

1.3 模型檢驗

選擇平均誤差(eˉ)，均方根誤差(RMSE)，總相對誤差(TRE)，確定系數(shù)(R2)4 種統(tǒng)計量來檢驗?zāi)Ｐ偷臄M合優(yōu)度和誤差，并觀察殘差分布是否有異質(zhì)性。利用修正確定系數(shù)AjustedR2對模型進行檢驗和評價（式22）。

1.4 模型比較

對模型間的差異采用F 統(tǒng)計檢驗，F(xiàn)統(tǒng)計指標(biāo)按照下式計算：

其中：sses和dfs分別為Richards 模型1 的回歸模型殘差平方和和自由度，ssec和dfc分別為Richard模型2 的回歸模型殘差平方和和自由度。

利用上式計算的F值與自由度分別為dfsdfc和dfs的F臨界值進行比較，以確定兩個模型之間是否存在顯著差異。

2 結(jié)果

表3 是利用最小二乘法擬合的模型式（13）～（15）、式（12）和式（4）的參數(shù)值。從表中可以看到，所有模型的漸近線A值在17.53～22.72 之間，年齡無關(guān)的平均高生長模型式（13）～（15）和式（12）的漸近線在21.85～22.72 之間，更符合落葉松最大高的實際情況。而年齡相關(guān)的模型式（4）的漸近線只有17.53，明顯偏低。

在表4 中，利用刀切法對模型進行驗證，每次從建模樣本中隨機抽取50% 的樣本，使用同一模型擬合，重復(fù)10 次，進行模型參數(shù)的靈敏度分析來評價模型；擬合結(jié)果使用驗證樣本年齡無關(guān)的模型式（13）～（16）和式（12）的修正確定系數(shù)在0.92～0.94 之間，而年齡相關(guān)的模型修正確定系數(shù)為0.85，時間無關(guān)的模型修正確定系數(shù)提高了8%～10%；模型式（13）～（16）和式（12）的均方根誤差RMSE在1.30～1.59 之間，遠遠小于年齡相關(guān)的模型式（4）的20.48。如圖1a 和圖1b 所示，時間和平均高的關(guān)系很難用一條曲線來描述，而用年齡無關(guān)的模擬擬合的平均高和實測平均高相關(guān)系數(shù)卻高達0.97。

表3 模型參數(shù)Table 3 Model parameters

表4 線性函數(shù)參數(shù)Table 4 Parameters of linear function

年齡無關(guān)的幾個模型進行對比，從表5 的F值來看，F(xiàn)0.05在2.42～3.89 之間，F(xiàn)值在0.15～81.05之間。從表中可以看出，模型式（15）和式（16）的F值是0.15，F(xiàn)0.05是3.89，說明這兩個模型差異不顯著；同樣模型式（15）和式（13）的F值是2.98，F(xiàn)0.05是3.89，說明這兩個模型差異也不顯著。其它模型兩兩對比，都差異顯著。

模型式（13）和式（14）對比，式（14）引入了氣候因子，氣候因子的選用參考了José Tomé等[12] 的文章，從F值來看，差異顯著，同時式（14）的修正決定系數(shù)0.941 9，也高于式（13）的0.935 2（表6），而式（14）的均方根誤差RMSE和TRE較式（13）也有明顯的降低，說明氣候因子的引入提高了模型的模擬效果。圖2 是6 個模型的殘差分布圖，從圖上看出年齡無關(guān)的模型（式12～16）的殘差分布均優(yōu)于年齡相關(guān)的模型(式4)。

圖1 a 時間和平均高的關(guān)系圖 ，b 模型式（16）實測和預(yù)測數(shù)據(jù)關(guān)系圖Fig. 1 a The relationship between mean height and Age b The relationship between fitting data and measured data of formula 16

本研究在模型式（13）～（16）中采用自由參數(shù)k的形式，將k表達成立地因子和氣候因子的線性形式，而模型式（12）采用固定參數(shù)的形式，對比兩類模型發(fā)現(xiàn)，模型式（13）～（16）的修正可絕系數(shù)均大于式（12），而誤差均小于式（12）。

3 討論

本研究基于Richards 函數(shù)建立了年齡無關(guān)的差分方程，從而可用于年齡未知或者年齡不是有效因子的情況下的林分。選擇年齡已知的固定樣地，對比了年齡無關(guān)模型和年齡有關(guān)模型的模擬精度。年齡無關(guān)的生長模型比年齡相關(guān)的模型精度高而誤差小。葛宏立等(17)采用年齡無關(guān)模型方法研究馬尾松生長模型時預(yù)測精度高達99.80%，高東啟等[19]用年齡無關(guān)的模型擬合蒙古櫟生長，模型的R都大于0.98，預(yù)估精度在99.40%～99.72% 之間，邱思玉等[20]分別用不同的年齡無關(guān)樹高模型建立優(yōu)勢高模型也得到了較高的預(yù)估精度，研究結(jié)果與本研究基本一致。分析原因，可能是因為年齡無關(guān)的生長模型直接以單木自身的觀測值為自變量來預(yù)測后期的生長情況，這在一定程度上淡化了林相、立地、競爭等因子的影響。

表5 模型F 值Table 5 F parameter of models

表6 模型統(tǒng)計值Table 6 Model statistics

另外，從擬合得到的漸近線來看，年齡無關(guān)模型擬合的漸近線更符合落葉松高的實際情況，分析原因，作者認為一方面年齡無關(guān)的模型具有較好的預(yù)測能力。另一方面外業(yè)調(diào)查獲得的林分年齡存在一定的偏差，這也和年齡因子難于獲得有關(guān)。

圖2 模型殘差分布圖(a. 式12, b. 式13, c. 式14, d. 式15, e. 式16, f. 式4)Fig. 2 Distribution of residuals errors (a. Formula 12, b. Formula 13, c. Formula 14,d. Formula 15, e. Formula 16, f. Formula 4)

同時，本研究采用自由參數(shù)的方式將與年齡無關(guān)的立地因子和氣候因子引入模型的擬合中，對比固定參數(shù)和自由參數(shù)的兩類模型，可以發(fā)現(xiàn)自由參數(shù)的模型表現(xiàn)更好。因此認為參數(shù)k 與立地因子和氣候因子相關(guān)性強，這和Tomé等[15]的結(jié)論一致。引入氣候因子的模型比僅引入林分因子的模型表現(xiàn)更好，說明考慮氣候?qū)ιL的影響更有利于優(yōu)勢高生長的模擬。

但是，以5 a 為預(yù)測間隔期是比較短的，樹木生長短期內(nèi)不會有太大的生長變化，使得模型的預(yù)測誤差較小，從而保證了模型有較高的預(yù)測精度。模型仍然需要三期及以上的數(shù)據(jù)進行進一步的檢驗，從而分析在年齡間隔較大時能否保證較高的預(yù)測精度。