朱光玉 ，康 立 ，何海梅 ，呂 勇 ，銀勇平 ，吳 毅

（1.中南林業(yè)科技大學(xué)，湖南 長沙 410004；2.中國林業(yè)科學(xué)研究院 a.森林生態(tài)環(huán)境與保護(hù)研究所；b.資源信息研究所，北京 100091）

基于樹高-年齡分級(jí)的杉木人工林多形立地指數(shù)曲線模型研究

朱光玉1，2a，2b，康 立1，何海梅1，呂 勇1，銀勇平1，吳 毅1

（1.中南林業(yè)科技大學(xué)，湖南 長沙 410004；2.中國林業(yè)科學(xué)研究院 a.森林生態(tài)環(huán)境與保護(hù)研究所；b.資源信息研究所，北京 100091）

考慮相同或者相近年齡的杉木純林，其優(yōu)勢(shì)木高生長受立地條件的影響較大，對(duì)年齡、優(yōu)勢(shì)高進(jìn)行分級(jí)，基于樹高-年齡分級(jí)構(gòu)造啞變量，采用基于啞變量的非線性回歸分析，建立區(qū)域性杉木多形立地指數(shù)曲線模型，為區(qū)域性立地質(zhì)量評(píng)價(jià)和生產(chǎn)力預(yù)估提供研究思路和方法；同時(shí)提出一種合理的構(gòu)建樹高-年齡分級(jí)啞變量的思路和方法，并進(jìn)一步研究分級(jí)啞變量如何加在模型參數(shù)上較為合理。以湖南不同地區(qū)606組杉木純林優(yōu)勢(shì)木平均高-年齡數(shù)據(jù)為研究對(duì)象，采用理查德等4種常用非線性方程進(jìn)行模擬與分析，找出最優(yōu)基礎(chǔ)模型。采用循環(huán)迭代法結(jié)合最優(yōu)基礎(chǔ)模型，對(duì)樹高和年齡進(jìn)行分級(jí)，找出最佳的分級(jí)結(jié)果并構(gòu)建啞變量：(1)先給定初始齡階間距值3年，對(duì)各齡階樹高進(jìn)行20種分級(jí)方法，得到不同樹高分級(jí)方法所對(duì)應(yīng)的啞變量，采用含啞變量的非線性回歸分析方法進(jìn)行樹高-年齡相關(guān)關(guān)系模擬，并分析、對(duì)比建模精度與效果，找出最優(yōu)的樹高分級(jí)方法及結(jié)果；(2)基于最優(yōu)樹高分級(jí)結(jié)果，考慮17種不同的齡階劃分方法，得到17種不同齡階所對(duì)應(yīng)的啞變量，采用含啞變量的非線性回歸方程模擬樹高-年齡相關(guān)關(guān)系，并分析、對(duì)比建模精度與效果，找出最優(yōu)的齡階劃分方法及結(jié)果；(3)基于(2)中最優(yōu)齡階，重復(fù)1)的研究工作，得到最優(yōu)的樹高分級(jí)結(jié)果及所對(duì)應(yīng)的啞變量，即最優(yōu)啞變量?；谧顑?yōu)啞變量，考慮到啞變量加在不同的參數(shù)上模型模擬預(yù)測(cè)效果可能不一致，采用7種不同加啞變量的方法，進(jìn)行樹高-年齡生長曲線模擬，并對(duì)比、分析其模擬效果，模型檢驗(yàn)采用確定系數(shù)(R2)、絕對(duì)平均誤差(MAE)、均方根誤差(RMSE)和殘差4個(gè)指標(biāo)，找出最優(yōu)模型，并繪制多形立地指數(shù)曲線圖。通過對(duì)4種非線性模型的建模精度分析與比較得到理查德方程最優(yōu)，模型表達(dá)式為：H=a×[1-exp(-b×Age)]c，其確定系數(shù)(R2)為0.799 644；基于最優(yōu)模型，采用含啞變量的非線性回歸分析，利用循環(huán)迭代法分析不同樹高、年齡分級(jí)的模型模擬精度與結(jié)果，指出齡階間距為5、樹高分9級(jí)所得啞變量最優(yōu)；基于最優(yōu)啞變量，考慮7種啞變量模型模擬方法，得到啞變量加在參數(shù)a、b上最優(yōu)，其確定系數(shù)(R2)、絕對(duì)平均誤差(MAE)、均方根誤差(RMSE)分別為：0.953 029、0.836 93、1.119 938，殘差值均勻分布在橫軸(預(yù)測(cè)值)兩側(cè)，表明模型模擬效果較好。本研究采用含啞變量的回歸分析與循環(huán)迭代法，系統(tǒng)的構(gòu)建了湖南杉木多形立地指數(shù)曲線模型，不僅明顯提高了模型的精度：確定系數(shù)從0.799 644增加到0.953 029，而且為區(qū)域性立地質(zhì)量評(píng)價(jià)和生產(chǎn)力預(yù)估模型的構(gòu)建提供了可行思路和方法。

杉木人工林；樹高-年齡分級(jí)；多形立地指數(shù)模型；非線性回歸；啞變量

立地指數(shù)的概念最早產(chǎn)生于18世紀(jì)，并被作為評(píng)價(jià)立地生產(chǎn)潛力的指標(biāo)，即利用林分中平均優(yōu)勢(shì)高和年齡的關(guān)系，以基準(zhǔn)年齡時(shí)的優(yōu)勢(shì)高作為衡量立地生產(chǎn)潛力的指標(biāo)[1-3]。立地指數(shù)模型的建立研究是評(píng)價(jià)森林立地質(zhì)量的核心和關(guān)鍵。目前，國內(nèi)外學(xué)者廣泛采用立地指數(shù)法評(píng)價(jià)純林立地質(zhì)量[4-6]。立地指數(shù)曲線模型可分為同形（單形）和多形兩種，利用同形立地指數(shù)模型進(jìn)行小范圍內(nèi)、單一樹種間的立地質(zhì)量評(píng)價(jià)時(shí)可得到較好的模擬結(jié)果，但應(yīng)用到大區(qū)域甚至全國范圍內(nèi)時(shí)效果并不理想[7]。基于此科學(xué)問題，多形立地指數(shù)曲線模型的研究成為一種趨勢(shì)[8]。Bull利用多形曲線評(píng)價(jià)了紅松林立地質(zhì)量[9]；Carmean等在加拿大安大略省西北部建立了黑云杉和美國山楊多形立地指數(shù)模型[10]；Cieszewski等利用147塊樣地970組黑松數(shù)據(jù)建立了黑松多形立地指數(shù)模型，大幅提高了模型的精度和適用范圍[11]；Jerez-Rico M等基于柚木臨時(shí)樣地和固定樣地?cái)?shù)據(jù)，研究了基于非線性混合效應(yīng)的柚木立地指數(shù)模型，在委內(nèi)瑞拉西部平原取得較好的效果[12]。國內(nèi)，李希菲等將啞變量的回歸分析方法用于多形立地指數(shù)曲線構(gòu)建[13]；段愛國等以多期固定樣地?cái)?shù)據(jù)為基礎(chǔ)，采用差分方程法構(gòu)建了杉木人工林優(yōu)勢(shì)高生長模型，并探討了多形立地指數(shù)模型的構(gòu)建[14]。目前，關(guān)于立地指數(shù)曲線參數(shù)的求解主要有2種思路：一是利用差分方程法求解，其對(duì)數(shù)據(jù)要求比較嚴(yán)格，要求對(duì)固定樣地進(jìn)行多期觀測(cè)，對(duì)于區(qū)域性立地質(zhì)量評(píng)價(jià)，需要大量的財(cái)力、物力；另一種是采用經(jīng)驗(yàn)或者理論方程，利用回歸分析方法直接進(jìn)行參數(shù)求解，此方法與差分方程法相比，不僅數(shù)據(jù)要求沒那么嚴(yán)格，而且有大量成熟的統(tǒng)計(jì)軟件可以借助，簡單易行。

含啞變量的回歸分析與考慮混合效應(yīng)模型的回歸分析，被證明對(duì)于多形立地指數(shù)曲線模型的構(gòu)建是切實(shí)可行的，均是有效的[15]。但是對(duì)于區(qū)域性立地質(zhì)量評(píng)價(jià)與生產(chǎn)力預(yù)估，關(guān)于啞變量或混合效應(yīng)變量的構(gòu)建與劃分方法、以及啞變量或混合效應(yīng)變量在模型中如何表達(dá)仍需進(jìn)一步探索。

李?？仍跇?gòu)建全國主要樹種的樹高-胸徑曲線模型時(shí)，先對(duì)樹高和胸徑的分級(jí)得到啞變量，然后建立了含啞變量的樹高-胸徑曲線模型，取得了良好的效果，但是對(duì)于適宜啞變量的分級(jí)方法以及啞變量在模型中如何表達(dá)比較合理，并沒有進(jìn)行深入研究[16]。

第八次全國森林資源清查結(jié)果顯示，我國森林每公頃蓄積量為89.79 m3，遠(yuǎn)遠(yuǎn)低于歐美發(fā)達(dá)國家水平。我國森林質(zhì)量不高已成為影響林業(yè)發(fā)展的一個(gè)突出問題，未能做到“適地適樹”和依據(jù)立地條件開展森林經(jīng)營與培育是一個(gè)重要原因。因此，開展區(qū)域性森林立地質(zhì)量評(píng)價(jià)和生產(chǎn)力預(yù)估研究，為區(qū)域性森林經(jīng)營與培育提供基礎(chǔ)理論依據(jù)和支持，迫在眉睫。

基于上述現(xiàn)實(shí)與科學(xué)問題，本研究以湖南26個(gè)地區(qū)的杉木純林優(yōu)勢(shì)木平均高和年齡配對(duì)數(shù)據(jù)，開展含啞變量的多形立地指數(shù)曲線模型研究，為區(qū)域性森林立地質(zhì)量評(píng)價(jià)提供基礎(chǔ)理論研究。

1 材料與方法

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

以湖南杉木純林為對(duì)象，研究數(shù)據(jù)包括兩種類型：解析木和標(biāo)準(zhǔn)地（固定樣地和臨時(shí)樣地）數(shù)據(jù)，覆蓋湘東、南、西、北和湘中地區(qū)。其中：解析木數(shù)據(jù)采集于湖南衡陽縣岣嶁峰林場(chǎng)、湖南衡東縣紫金山和四方山林場(chǎng)；標(biāo)準(zhǔn)地?cái)?shù)據(jù)取自于湖南江華、雙排、汝城、通道、隆回、湘鄉(xiāng)、懷化、武岡、會(huì)同、龍山、新邵、株洲、黔陽、衡東、攸縣、臨湘、江永、城步、靖州、資興、瀘溪縣、綏寧、桃源、新邵等26縣。建模樣本數(shù)據(jù)共計(jì)606組。建模樣本的平均優(yōu)勢(shì)木年齡（Age）、和樹高（H）分布情況見表1。

表1 建模樣本統(tǒng)計(jì)Table 1 Summary statistics for modeling

1.2 研究方法

1.2.1 基礎(chǔ)模型選擇

采用4種常用的立地指數(shù)曲線模型[2-5]，對(duì)杉木優(yōu)勢(shì)木高和年齡的相關(guān)關(guān)系進(jìn)行模擬，模型表達(dá)式列于表2，其中：H為優(yōu)勢(shì)木高、Age為年齡；a、b、c為待求模型參數(shù)。對(duì)模型模擬結(jié)果進(jìn)行分析對(duì)比，選出最優(yōu)模型。模型擬合效果評(píng)價(jià)采用確定系數(shù)（R2）、平均絕對(duì)誤差（MAE）和均方根誤差（RMSE）3個(gè)指標(biāo)，選擇確定系數(shù)最大、平均絕對(duì)誤差和均方根誤差最小的模型作為基礎(chǔ)模型。

表2 4種樹高-年齡曲線模型Table 2 Four kinds of height-age models

1.2.2 啞變量分級(jí)及最優(yōu)啞變量的劃分

由于樹高受立地條件影響較大，從而導(dǎo)致相同或者相近年齡的樹木在不同的立地條件下樹高差別很大，基于此，對(duì)年齡進(jìn)行等級(jí)劃分（齡階），同時(shí)對(duì)不同齡階的樹高進(jìn)行分級(jí)，構(gòu)造反映不同立地條件對(duì)優(yōu)勢(shì)高影響的分級(jí)變量，并在構(gòu)建優(yōu)勢(shì)高及其年齡生長曲線時(shí)，將分級(jí)變量作為啞變量進(jìn)行回歸分析，求解立地指數(shù)曲線模型。對(duì)樹高和年齡分級(jí)時(shí)采用循環(huán)迭代法，并利用含啞變量的回歸分析方法判斷、比較年齡與樹高分級(jí)的結(jié)果，探索、找出最佳的年齡和樹高分級(jí)方法。

1.2.2.1 年齡、優(yōu)勢(shì)高分級(jí)啞變量構(gòu)造及選優(yōu)

對(duì)于建模樣本，令其樣本量為n（n=606），對(duì)于第i個(gè)樣本，其優(yōu)勢(shì)高為hi，其年齡為Agei。如果分別將樣本中的年齡分m（m≥1）級(jí)，優(yōu)勢(shì)高分為p級(jí)；則對(duì)于第i個(gè)樣本，其齡階標(biāo)志為其優(yōu)勢(shì)高分級(jí)標(biāo)志為均為整數(shù)。

由于不同樹高分級(jí)數(shù)可能對(duì)模型擬合結(jié)果有不同的影響，且樹高分級(jí)數(shù)過大不利于實(shí)際應(yīng)用，本研究對(duì)杉木最優(yōu)樹高分級(jí)數(shù)進(jìn)行了探索研究。在本研究中給定的樹高分級(jí)數(shù)m的定義域：m為整數(shù)，最小值為2，最大值為平均分級(jí)間距接近0.5時(shí)所對(duì)應(yīng)的分級(jí)數(shù)。

樹高級(jí)確定后，不同立地條件下，不同的齡階劃分方法有著不同的擬合結(jié)果。在本研究中，樣本年齡的最大值為49 a，所以定義齡階間距最小值為2 a，齡階間距最大值為50 a。由于各齡階樣本數(shù)量分布不均勻，可能會(huì)出現(xiàn)某個(gè)齡階只有一個(gè)樣本的情況，此時(shí)該齡階最大樹高值與最小樹高值相等，分級(jí)間距為0，從而導(dǎo)致無法正確判斷該樣本的樹高等級(jí)。針對(duì)這種特殊情況，本研究對(duì)于齡階中只有一個(gè)樣本的數(shù)據(jù)不參與建模，而是在建模后將該樣本代入各樹高等級(jí)曲線中，將該樣本劃分到預(yù)測(cè)值與實(shí)測(cè)值最接近的樹高等級(jí)曲線。

具體算法如下：

（1）對(duì)所有Agei，按l年一個(gè)齡階等分為m個(gè)齡階，得到如：齡階為5 a，則1—5 a為第一個(gè)齡階，6—10 a為第二個(gè)齡階，以此類推；

（3）對(duì)每個(gè)hi，計(jì)算其hoi：

（1）給定初始齡階間距為3 a（我國杉木解析木齡階間距通常為3 a），將平均優(yōu)勢(shì)木年齡按3 a一個(gè)齡階整化為Q個(gè)齡階；將相同齡階的樹高等分為m級(jí)，即將相同齡階的樹高最大值與樹高最小值的差值等分為m級(jí)。根據(jù)上述啞變量分級(jí)方法，可得每個(gè)樣本觀測(cè)點(diǎn)的初次樹高分級(jí)，1級(jí)是最低等級(jí)樹高曲線。利用統(tǒng)計(jì)之林軟件中非線性回歸模塊，按不同樹高級(jí)劃分得到的不同啞變量，進(jìn)行含啞變量的非線性回歸分析?；A(chǔ)模型以理查德模型為例：

將不同樹高級(jí)劃分后的啞變量DJm加到基礎(chǔ)模型參數(shù)a上，模型變?yōu)椋?/p>

式中：a、b、c為模型參數(shù)，H為平均優(yōu)勢(shì)木高，Age為平均優(yōu)勢(shì)木年齡，DJm為樹高級(jí)為m時(shí)對(duì)應(yīng)的各樣本樹高等級(jí)。選取模擬效果最好的樹高級(jí)S0級(jí)。

（2）利用啞變量分級(jí)方法，按樹高級(jí)S0級(jí)，不同齡階N年再次進(jìn)行啞變量劃分，可以得到不同的啞變量分級(jí)結(jié)果。

利用統(tǒng)計(jì)之林軟件中非線性回歸模塊，按不同齡階劃分得到的不同啞變量結(jié)果，進(jìn)行回歸分析，將不同齡階劃分得到的啞變量DJN加到基礎(chǔ)模型式參數(shù)a中，模型變?yōu)椋?/p>

式中：a、b、c為模型參數(shù)，H為優(yōu)勢(shì)木平均高，Age為平均優(yōu)勢(shì)木年齡，DJN為齡階N年時(shí)對(duì)應(yīng)的各樣本樹高等級(jí)。選取模擬效果最好的齡階間距K年。

（3）確定齡階間距為K0年后，重復(fù)（1）的工作，選出最優(yōu)樹高級(jí)值S1。

（4）以樹高分級(jí)級(jí)數(shù)S1為基礎(chǔ)，繼續(xù)對(duì)齡階間距進(jìn)行劃分，選出最優(yōu)的齡階間距為K1。

（5）重復(fù)(3)和(4)，直至齡階間距K0=K1，為常數(shù)，樹高分級(jí)級(jí)數(shù)，S1=S0，為常數(shù)，則終止循環(huán)，得到最優(yōu)的啞變量劃分結(jié)果，其齡階間距為K1，其樹高分級(jí)級(jí)數(shù)為S1，所對(duì)應(yīng)樹高-年齡分級(jí)啞變量為最優(yōu)啞變量，最優(yōu)啞變量的級(jí)數(shù)與S1一致。

1.2.2.2 最優(yōu)啞變量模型選擇

為了進(jìn)一步研究啞變量如何加在參數(shù)上較為合理，按照齡階K年、樹高級(jí)值S0級(jí)的啞變量劃分結(jié)果，將啞變量DJ分別加在基礎(chǔ)模型的不同參數(shù)上，進(jìn)行7種情況的建模分析與比較，通過對(duì)模型進(jìn)行建模精度檢驗(yàn)，選出最優(yōu)啞變量模型。7種模型的表達(dá)式如下：

1.2.3 模型檢驗(yàn)

模型檢驗(yàn)采用確定系數(shù)（R2）、平均絕對(duì)誤差（MAE）、均方根誤差（RMSE）和殘差4個(gè)評(píng)價(jià)指標(biāo)，判別標(biāo)準(zhǔn)為：R2越大，精度越高；MAE、RMSE越小，精度越高；殘差值分布越均勻，模型擬合效果越好。計(jì)算公式如下：

其中：hi為樹高實(shí)測(cè)值，為樹高預(yù)測(cè)值，為樹高實(shí)測(cè)值的平均值，N為所有樣本總數(shù)。

2 結(jié)果與分析

2.1 基礎(chǔ)模型選取

利用統(tǒng)計(jì)之林軟件中一元非線性回歸模塊對(duì)4種候選模型進(jìn)行擬合分析，將擬合結(jié)果列于表3。結(jié)果指出理查德模型確定系數(shù)（R2）在4種候選模型中最大，為0.799 644；平均絕對(duì)誤差（MAE）、均方根誤差在4者中最小，分別為1.480 003、2.313 032，故選用理查德模型作為構(gòu)建多形立地指數(shù)模型的基礎(chǔ)模型。

表3 4種候選模型的擬合結(jié)果Table 3 Fitting results of four candidate models

2.2 基于循環(huán)迭代法的最優(yōu)啞變量的劃分

（1）令初始齡階為3 a，將平均優(yōu)勢(shì)木年齡按3 a一個(gè)齡階整化為Q個(gè)齡階；將相同齡階的樹高等分為m級(jí)，即將相同齡階的樹高最大值與樹高最小值的差值等分為m級(jí)；如表4所示，樹高分級(jí)數(shù)為21時(shí)，其平均分級(jí)間距為0.491，與本研究要求的精度（0.5 m 1級(jí)）是一致的，所以最低樹高級(jí)值為2，最高樹高級(jí)值為21，一共有20種不同的啞變量構(gòu)造結(jié)果。

表4 樹高級(jí)間距值Table 4 Interval value of height classification

利用統(tǒng)計(jì)之林軟件的非線性回歸模塊，將20種不同的啞變量構(gòu)造結(jié)果DJ，加在基礎(chǔ)模型的參數(shù)a上分別進(jìn)行擬合，建模精度檢驗(yàn)采用MAE、RMSE和確定系數(shù)（R2）進(jìn)行比較，20次擬合的精度結(jié)果如表5：

表5 不同樹高分級(jí)數(shù)擬合精度Table 5 Fitting accuracy of different height classification

由各樹高級(jí)值擬合結(jié)果可知：樹高分9級(jí)時(shí)，其確定系數(shù)最大、MAE和RMSE的值最小，表明分9級(jí)時(shí)，模型模擬最優(yōu)，所以選用9級(jí)為初始樹高級(jí)。

（2）利用啞變量分級(jí)方法，按樹高級(jí)9級(jí)，不同齡階N年再次進(jìn)行啞變量劃分，可以得到17種不同的啞變量分級(jí)結(jié)果。

利用統(tǒng)計(jì)之林軟件中非線性回歸模塊，加在基礎(chǔ)模型的參數(shù)a上分別對(duì)17種模型進(jìn)行擬合分析，擬合結(jié)果見表6。

由擬合結(jié)果可知，齡階間距為5 a時(shí)，其確定系數(shù)最大，MAE和RMSE的值最小，模型擬合效果最好，所以初步確定齡階間距為5 a最優(yōu)。

（3）確定齡階間距為5年后，重復(fù)（1）的研究，20次模型擬合的結(jié)果見表7。

擬合結(jié)果表明：樹高級(jí)為9級(jí)時(shí)，其確定系數(shù)最大，并且其MAE和RMSE的值最小，表明其模型模擬效果最好，與（1）結(jié)果一致，所以，再次確定樹高分9級(jí)最優(yōu)，再基于最優(yōu)樹高級(jí)數(shù)劃分齡階，將得到（2）的齡階劃分結(jié)果。至此，已選出最優(yōu)齡階間距為5 a，最優(yōu)樹高分級(jí)級(jí)數(shù)為9級(jí)。

表6 不同齡階擬合精度Table 6 Fitting accuracy of different age class

表7 不同樹高分級(jí)數(shù)擬合精度Table 7 Fitting accuracy of different height classification

通過循環(huán)迭代的年齡樹高分級(jí)和回歸模擬比較分析可知：齡階間距為5 a、樹高分9級(jí)所構(gòu)造的啞變量DJ最優(yōu)。

2.3 最優(yōu)含啞變量模型的選取與多形立地指數(shù)曲線模型構(gòu)建

考慮到分級(jí)啞變量加在模型的不同參數(shù)上，對(duì)模型模擬效果可能有影響，將啞變量DJ分別加在基礎(chǔ)模型的不同參數(shù)上，利用統(tǒng)計(jì)之林含啞變量的非線性回歸分析模塊，進(jìn)行了7種情況的建模分析，探索最優(yōu)啞變量模型（見圖1）。

（1）啞變量DJ加在參數(shù)a上

其模型表達(dá)式為：H=(a×DJ)(1-e-b×Age)c。

通過模擬，模型的確定系數(shù)(R2)為0.945 278，平均絕對(duì)誤差（MAE）為0.861 106，均方根誤差（RMSE）為1.178 589，擬合結(jié)果見表8。

圖1 啞變量DJ加在a上模型的樹高實(shí)測(cè)值-預(yù)測(cè)值與預(yù)測(cè)值-殘差Fig.1 Observed value-predicted value and predicted value-residuals of add dummy variable to parameter a

表8 啞變量加在參數(shù)a上的模型模擬參數(shù)Table 8 Model Simulation parameters with adding dummy variable to parameter a

（2）啞變量DJ加在參數(shù)b上

其模型表達(dá)式為：H=(a×DJ)(1-e-b×DJ×Age)c。

通過擬合可得，模型參數(shù)a的擬合檢驗(yàn)值P值等于1，表明參數(shù)a的擬合沒有達(dá)到預(yù)期，結(jié)果不可信（見表9）。

（3）啞變量DJ加在參數(shù)c上

其模型表達(dá)式為：H=(a×DJ)(1-e-b×Age)c×DJ。

通過擬合得，模型的確定系數(shù)（R2）為0.924 894，平均絕對(duì)誤差（MAE）為1.041 839，均方根誤差（RMSE）為1.416 174，擬合結(jié)果見表10和圖2。

（4）啞變量DJ加在參數(shù)a、b上：

其模型表達(dá)式為：H=(a×DJ)(1-e-b×DJ×Age)c。

通過模型擬合，模型的確定系數(shù)（R2）為0.953 029，平均絕對(duì)誤差（MAE）為0.836 93，均方根誤差（RMSE）為1.119 938，擬合結(jié)果見表11和圖3。

表9 啞變量加在參數(shù)b上的模型模擬參數(shù)Table 9 Model simulation parameters with adding dummy variable to parameter b

表10 啞變量加在參數(shù)c上的模型模擬參數(shù)Table 10 Model simulation parameters with adding dummy variable to parameter c

圖2 啞變量DJ加在參數(shù)c上模型的樹高實(shí)測(cè)值-預(yù)測(cè)值與預(yù)測(cè)值-殘差Fig.2 Observed value--predicted value and predicted value-residuals of add dummy variable to parameter c

（5）啞變量DJ加在參數(shù)b、c上

表11 啞變量加在參數(shù)a、b上的模型模擬參數(shù)Table 11 Model simulation parameters with adding dummy variable to parameters a and b

圖3 啞變量DJ加在參數(shù)a、b上模型的樹高實(shí)測(cè)值-預(yù)測(cè)值與預(yù)測(cè)值-殘差Fig.3 Observed value-predicted value and predicted value-residuals of add dummy variable to parameters a, b

其模型表達(dá)式為：H=(a×DJ)(1-e-b×DJ×Age)c×DJ。

通過擬合可得，模型參數(shù)a的擬合檢驗(yàn)值P值等于1，表明參數(shù)a的擬合沒有達(dá)到預(yù)期，結(jié)果不可信（見表12）。

表12 啞變量加在參數(shù)b、c上的模型模擬參數(shù)Table 12 Model simulation parameters with adding dummy variable to parameters b and c

（6）啞變量DJ加在參數(shù)a、c上

其模型表達(dá)式為：H=(a×DJ)(1-e-b×Age)c×DJ。

通過模型擬合，模型的確定系數(shù)（R2）為0.949 347，平均絕對(duì)誤差MAE為0.839 631，均方根誤差RMSE為1.139 34。參數(shù)估計(jì)和預(yù)測(cè)值分布情況見表13。

表13 啞變量加在參數(shù)a、c的模型模擬參數(shù)Table 13 Model simulation parameters with adding dummy variable to parameters a and c

（7）啞變量DJ加在參數(shù)a、b、c上

其模型表達(dá)式為：H=(a×DJ)(1-e-b×DJ×Age)c。

通過擬合可得，模型參數(shù)a的第9個(gè)等級(jí)的擬合檢驗(yàn)值P值等于1，表明參數(shù)a的擬合沒有達(dá)到預(yù)期，結(jié)果不可信（見表14）。

（8）最優(yōu)啞變量模型選取

分別在不同參數(shù)上添加啞變量，利用統(tǒng)計(jì)之林中的非線性回歸模塊進(jìn)行7次不同的擬合分析。從確定系數(shù)、平均絕對(duì)誤差、均方根誤差、樹高實(shí)測(cè)值-預(yù)測(cè)值和樹高預(yù)測(cè)值來看，在參數(shù)a、b上添加啞變量結(jié)果最優(yōu)，所以最優(yōu)啞變量模型為：H=(a×DJ)(1-e-b×DJ×Age)c。

（9）多形立地指數(shù)曲線模擬與分析

圖4 啞變量DJ加在參數(shù)a、c上模型的樹高實(shí)測(cè)值-預(yù)測(cè)值與預(yù)測(cè)值-殘差Fig.4 Observed value-predicted value and predicted value - residuals of add dummy variable to parameters a, c

表14 啞變量加在參數(shù)a、b、c上的模型模擬參數(shù)Table 14 Model simulation parameters with adding dummy variable to parameters a, b and c

利用最優(yōu)模型 H=(a×DJ)(1-e-b×DJ×Age)c，模擬研究區(qū)杉木平均優(yōu)勢(shì)木樹高和年齡的相關(guān)關(guān)系，參數(shù)預(yù)估值見表12， 模型的確定系數(shù)（R2）為0.953 029，平均絕對(duì)誤差（MAE）為0.836 93，均方根誤差（RMSE）為1.119 938?；诖藢?dǎo)向曲線繪制其立地質(zhì)量預(yù)估效果圖，即多形立地指數(shù)曲線圖，圖5中的散點(diǎn)表示觀測(cè)值，即建模用數(shù)據(jù)。

圖5 多形立地指數(shù)曲線Fig.5 Polymorphic site index curve

各立地質(zhì)量等級(jí)對(duì)應(yīng)的樹高-年齡散點(diǎn)圖如圖6。

通過各立地等級(jí)對(duì)應(yīng)的樹高-年齡散點(diǎn)圖，可知各個(gè)立地等級(jí)的立地指數(shù)曲線模擬效果均較好，體現(xiàn)了模型精度較高的特點(diǎn)，進(jìn)一步證實(shí)了模型的可行性和適用性，同時(shí)也驗(yàn)證了啞變量劃分與選擇的合理性與可行性；啞變量級(jí)數(shù)值越高，優(yōu)勢(shì)木生長越好，立地條件越好，立地質(zhì)量越高。

3 結(jié) 論

本研究以湖南東、南、西、北、中26個(gè)縣的杉木人工純林為對(duì)象，收集包括固定樣地、臨時(shí)樣地和解析木等606組樣品，采用循環(huán)迭代的思維和含啞變量的非線性回歸分析方法，構(gòu)建了湖南杉木人工林多形立地指數(shù)曲線模型。

采用理查德（Richards）、坎派茲（Gompertz）、考爾夫（Korf）和單分子（Mitscherlich）4種非線性模型對(duì)優(yōu)勢(shì)木高和年齡的相關(guān)關(guān)系進(jìn)行模擬，選出最優(yōu)模型為理查德（Richards）函數(shù)：H=a×[1-exp(-b×Age)]c。

圖6 9種樹高-年齡分級(jí)曲線Fig.6 Nine kinds of curve with classi fi cation of height and age

考慮到不同條件下杉木生產(chǎn)潛力的差異性，首次采用循環(huán)迭代的思維對(duì)平均優(yōu)勢(shì)木的年齡和樹高進(jìn)行等級(jí)劃分，借助含啞變量的非線性回歸進(jìn)行分級(jí)結(jié)果檢驗(yàn)，指出年齡分級(jí)時(shí)，齡階間距為5年，樹高分9級(jí)，其所得的分級(jí)啞變量最佳，分級(jí)最佳啞變量即為杉木立地質(zhì)量高低定性劃分的依據(jù)，分級(jí)啞變量的等級(jí)值越大，優(yōu)勢(shì)木生長越好，立地質(zhì)量越好。

考慮樹高-年齡分級(jí)啞變量對(duì)理查德（Richards）生長方程3個(gè)參數(shù)的影響可能不一致，采用了啞變量分別加在不同參數(shù)上的7種情況進(jìn)行了7種含啞變量的非線性回歸模型模擬，并利用模型的確定系數(shù)（R2）、平均絕對(duì)誤差（MAE）、均方根誤差（RMSE）等，來比較分析不同模型的精度，指出最優(yōu)模型為：H=(a×DJ)(1-e-b×DJ×Age)c。并繪制了該模型的多形立地指數(shù)曲線圖，及分9級(jí)的各級(jí)曲線圖，圖示表明模型模擬效果較好。與基礎(chǔ)模型相比，將模型的確定系數(shù)從0.799 644增至0.953 029，大大提高了模型的精度。同時(shí)也驗(yàn)證了樹高-年齡分級(jí)的方法是正確、可行的。

研究結(jié)果為區(qū)域性立地指數(shù)曲線的構(gòu)建提供一種新思路，有助于區(qū)域性立地質(zhì)量評(píng)價(jià)。

[1]Nanos N, Calama R, Ca?adas N, et al. Spatial stochastic modelling of cone production from stone pine (Pinus pinea L.)stands in the Spanish Northern Plateau[C]// Modelling forest systems. Workshop on the interface between reality, modelling and the parameter estimation processes, Sesimbra, Portugal, 2-5 June 2002.

[2]Payandeh B, Wang Y. Relative accuracy of a new base-age invariant site index model[J]. Forest Science, 1994, 40(2):341-348.

[3]Harry Eriksson, Ulf Johansson, Andres Kiviste. A site-index model for pure and mixed stands of Betula pendula and Betula pubescens in Sweden[J]. Scandinavian Journal of Forest Research, 1997, 12(2):149-156.

[4]Yue C, M?kinen H, Kl?dtke J, et al. An approach to assessing site index changes of Norway spruce based on spatially and temporally disjunct measurement series[J]. Forest Ecology &Management, 2014, 323(7):10-19.

[5]Batho A, García O. A Site Index Model for Lodgepole Pine in British Columbia[J]. Forest Science, 2014(2): 60.

[6]Kimberley M O, Ledgard N J. Site index curves for Pinus nigra grown in the South Island high country, New Zealand[J]. New Zealand Journal of Forestry Science, 1998(1): 389-399.

[7]Devan J S, Burkhart H E. Polymorphic site index equations for loblolly pine based on a segmented polynomial differential model[J]. Forest Science, 1982, 28(3):544-555.

[8]陳紹玲. 馬尾松人工林多形地位指數(shù)曲線模型的建模方法[J].中南林業(yè)科技大學(xué)學(xué)報(bào), 2008, 28(2): 125-128.

[9]Bull H. The use of polymorphic curves in determining site quality in young red pine plantations[M]. US Government Printing Of fi ce, 1931.

[10]Carmean W H, Hazenberg G, Deschamps K C. Polymorphic site index curves for black spruce and trembling aspen in northwest Ontario[J]. Forestry Chronicle, 2006, 82(2): 231-242.

[11]Cieszewski C J, Bella I E. Polymorphic height and site index curves for lodgepole pine in Alberta[J]. Canadian Journal of Forest Research, 2011, 19(9): 1151-1160.

[12]Jerez-Rico M, Moret-Barillas A Y, Carrero-Gámez O E, et al. Site index curves based on mixed models for teak (Tectona grandis L. F.) plantations in the Venezuelan plains[J]. Agrociencia, 2011,45(1):135-145.

[13]李希菲,洪玲霞.用啞變量法求算立地指數(shù)曲線族的研究[J].林業(yè)科學(xué)研究, 1997,10(2):108-112.

[14]段愛國,張建國.杉木人工林優(yōu)勢(shì)高生長模擬及多形地位指數(shù)方程[J]. 林業(yè)科學(xué), 2004,40(6):13-19.

[15]Palahí M, Tomé M, Pukkala T, et al. Site index model for Pinus sylvestris in north-east Spain[J]. Forest Ecology & Management,2004, 187(1):35-47.

[16 李?？?法 蕾.基于分級(jí)的全國主要樹種樹高-胸徑曲線模型[J].林業(yè)科學(xué), 2011,47(10):83-90.

Study on polymorphic site index curve model based on height-age classi fi cation for Cunninghamia lanceolata plantation

ZHU Guangyu1,2a,2b, KANG Li1, HE Haimei1, LV Yong1, YIN Yongping1, WU Yi1
(1.College of Forestry, Central South University of Forestry and Technology, Changsha 410004, Hunan, China; 2a. Institute of Forest Ecology, Environment and Protection; 2b. Research Institute of Forestry Resource Information Techniques, Chinese Academy of Forestry, Beijing 100091, China)

In this study, the tree age and its dominant height of Cunninghamia lanceolata plantation are fi rstly classi fi ed in consideration of the pure stands of Cunninghamia lanceolata in the same or similar age whose growth of dominant tree is greatly affected by site conditions. Secondly, the dummy variable will be constructed on the basis of classification of tree height and age, Finally, the polymorphic site index curve model for regional Cunninghamia lanceolata will be established, employing the non-linear regression analysis based on the dummy variable. The above research will provide new ideas and methods for evaluating forest site quality and estimating productivity in regional forest. In the same time, a new reasonable method of constructing dummy variables of tree height and age classi fi cation is also proposed and how to reasonably add the dummy variable to the model will be further studied in the research. In this study,the data of average height and age of dominant tree is taken as the study subject, which are from 606 groups of Cunninghamia lanceolata pure forest in different regions of Hunan province. In order to fi nd the optimal basic model, the commonly used nonlinear equations as Richards, Gompertz, Korf and Mitscherlich are employed to simulate and analyze. In the purpose of fi nding out the ideal results of classi fi cation to construct the dummy variables, the cyclic iterative method integrated with the optimal basic model is taken to classify the tree height and age. The speci fi c implementation steps are as follows: Firstly, setting the three-year as the initial tree age graduation, classifying the tree height of various age graduation by nine methods, then the corresponding dummy variables to different trees’ classi fi cation methods will be obtained, employing non-linear regression analysis method with the dummy variables to simulate and analyze the relationship between tree height and tree age, comparing the accuracy and validity of models to fi nd out the optimal classifying method and the ideal result of tree height. Secondly, 17 kinds of different tree age classi fi cation methods and their corresponding dummy variables are obtained from the above ideal result, employing non-linear regression analysis method with the dummy variables to simulate and analyze the relationship between tree height and tree age, comparing the accuracy and validity of models to fi nd out the optimal classifying method and the ideal result of tree-age graduation. Thirdly, the optimal corresponding dummy variable and tree-height classi fi cation are obtained from the above optimal tree-age graduation by repeating the fi rst step. Based on the optimal dummy variable, the growth curve of the tree height and age is simulated by taking seven different methods with the dummy variable in consideration of that adding the dummy variable to different parameters’ model to simulate would obtain different results.And then comparing and analyzing these simulation effects to fi nd out the optimal model to draw the polymorphic site index curve. The model is tested by the following four indexes as the determine coef fi cient(R2), mean absolute error(MAE), root mean square error(RMSE)and the residual error. The fi rst one is that the Richard equation is the optimal one by comparing and analyzing the modeling accuracy of the four kinds of non-linear models, and its expression is ,whose determinant parameter R2 is 0.79964,And the second one is that the optimal dummy variable is obtained when the tree-age graduation set as fi ve-year and the tree-height classi fi ed into nine classi fi cations by adopting the non-linear regression with dummy variables and employing the cyclic iterative methods to analyze the model simulation accuracy and validity of different tree height classi fi cation and tree age graduation on the basis of the optimal model. And the third one is that the dummy variable added to the parameters of a、b is the optimal one by simulating seven kinds of models with dummy variables on the basis of the optimal dummy variable, whose determinant coef fi cient(R2), mean absolute error(MAE), root mean square error(RMSE)are 0.953 029、0.836 93、1.119 938 respectively, and the residual values evenly distribute on both sides of the horizontal axis, which shows that the validity of the model is good. In this study, the polymorphic site index curve model for Cunninghamia lanceolata of Hunan province is systematically established by employing the methods of regression analysis with dummy variables and cyclic iteration method,which not only increased the accuracy of determinant coef fi cients from 0.799 644 to 0.953 029, but also provide new ideas and methods for constructing the model of evaluating forest site quality and estimating productivity prediction in regional places.

Cunninghamia lanceolata plantation; height-age classi fi cation; polymorphic site index model; nonlinear regression; dummy variable

S757.2

A

1673-923X(2017)07-0018-12

10.14067/j.cnki.1673-923x.2017.07.003

2016-07-03

國家自然科學(xué)基金（31570631，31100476）；國家林業(yè)局項(xiàng)目（1692016-06，SFA2130218）；中國博士后科學(xué)基金面上項(xiàng)目（2014M550103）；湖南省教育廳項(xiàng)目（17C1664）

朱光玉，副教授，博士；E-mail：zgy1111999@163.com

朱光玉，康 立，何海梅，等.基于樹高-年齡分級(jí)的杉木人工林多形立地指數(shù)曲線模型研究[J].中南林業(yè)科技大學(xué)學(xué)報(bào)，2017, 37(7): 18-29.

[本文編校：吳 毅]