董云飛 孫玉軍 許 昊 梅光義

(北京林業(yè)大學(xué)，北京，100083)

責(zé)任編輯:王廣建。

樹(shù)高是森林經(jīng)營(yíng)管理以及生長(zhǎng)收獲研究中必需的因子。森林調(diào)查中，通常是先測(cè)量部分樣木的胸徑與樹(shù)高來(lái)建立樹(shù)高曲線，然后對(duì)缺失的樹(shù)高進(jìn)行預(yù)測(cè)。普通的樹(shù)高曲線，僅以胸徑為自變量，不適用于各種不同類型的林分，需要為每個(gè)林分建立不同的模型，因此，普通的樹(shù)高曲線的應(yīng)用范圍非常局限，而標(biāo)準(zhǔn)樹(shù)高曲線，以胸徑和樹(shù)木及林分因子為自變量，可以用于更廣的區(qū)域［1-3］。然而樹(shù)高曲線隨林分的變化而變化，并且同一林分的樹(shù)高與胸徑的關(guān)系隨年齡也發(fā)生變化［4］。另一方面用于建立樹(shù)高曲線的數(shù)據(jù)通常來(lái)源于不同林分內(nèi)的多個(gè)不同的樣地，或是固定樣地內(nèi)不同期間的多次測(cè)量。這樣的隨機(jī)嵌套結(jié)構(gòu)導(dǎo)致數(shù)據(jù)缺乏獨(dú)立性。目前利用非參數(shù)估計(jì)方法及廣義加性模型(GAM)方法構(gòu)造的樹(shù)高曲線方程也只能反映平均水平，而不能體現(xiàn)樣地或林分之間的差異性?；旌夏Ｐ偷陌l(fā)展為詳細(xì)描述這種嵌套的隨機(jī)結(jié)構(gòu)提供了可行的統(tǒng)計(jì)方法，并且構(gòu)造的隨機(jī)效應(yīng)參數(shù)可以校正或縮小樣地之間及林分之間的差異［5］。混合模型在預(yù)估因變量時(shí)，可以通過(guò)測(cè)量某一個(gè)樣地少量樣本來(lái)預(yù)測(cè)這個(gè)樣地的隨機(jī)參數(shù)向量，讓大區(qū)域的標(biāo)準(zhǔn)樹(shù)高曲線運(yùn)用到某一個(gè)具體的樣地，因此，可以運(yùn)用到森林資源清查和經(jīng)營(yíng)決策中［6］。

國(guó)外開(kāi)展了大量的樹(shù)高曲線的非線性混合模型研究［7-9］，國(guó)內(nèi)已經(jīng)有許多關(guān)于杉木、思茅松、馬尾松、楊樹(shù)、云杉、冷杉和紅松等標(biāo)準(zhǔn)樹(shù)高曲線的研究［10-13］，但基于非線性混合模型的標(biāo)準(zhǔn)樹(shù)高曲線的研究還很少，李春明等［5］基于非線性混合模型研究了栓皮櫟的樹(shù)高和胸徑的關(guān)系。

本文選取了在將樂(lè)縣國(guó)有林場(chǎng)設(shè)置的29 塊杉木純林樣地(2 577 株)，并將這29 塊樣地分為兩部分:一部分(20 塊)作為建模數(shù)據(jù)，另一部分(9 塊)作為檢驗(yàn)數(shù)據(jù);先利用建模數(shù)據(jù)和非線性最小二乘法選出基礎(chǔ)模型，然后結(jié)合建模數(shù)據(jù)構(gòu)建非線性混合效應(yīng)模型，通過(guò)選擇混合參數(shù)，選擇模型收斂且擬合優(yōu)度最好的混合模型作為最優(yōu)模型，并利用檢驗(yàn)數(shù)據(jù)，及選擇相關(guān)指標(biāo)，與傳統(tǒng)的非線性最小二乘法進(jìn)行精度比較。

1 研究區(qū)概況

實(shí)驗(yàn)區(qū)位于福建省將樂(lè)縣國(guó)有林場(chǎng)，地處武夷山脈東南部，地理坐標(biāo)為26°26'～27°04'N，117°05'～117°40'E。以中、低山為主，海拔高度180～500 m。屬亞熱帶季風(fēng)氣候，具有海洋性和大陸性氣候特點(diǎn)，年均氣溫18.7 ℃，年降水量1 672.3 mm，平均相對(duì)濕度為81%。境內(nèi)氣溫較高，夏季時(shí)間長(zhǎng)，冬天較暖和，霜凍較少，生長(zhǎng)期長(zhǎng)［14］。

2 研究方法

本文選取了在將樂(lè)縣國(guó)有林場(chǎng)設(shè)置的29 塊杉木純林樣地(2 577 株)，方形樣地面積為(20 m ×20 m 或20 m ×30 m)0.04 hm2或0.06 hm2，實(shí)測(cè)樣地樹(shù)木的胸徑和樹(shù)高，起測(cè)胸徑為2.0 cm，樹(shù)高精度為0.1 m(見(jiàn)表1)。

表1 建模和檢驗(yàn)數(shù)據(jù)概況

根據(jù)樣地所有杉木的實(shí)測(cè)數(shù)據(jù)，得到樹(shù)高和胸徑關(guān)系的散點(diǎn)(見(jiàn)圖1)。

圖1 樹(shù)高與胸徑的關(guān)系

2.1 基礎(chǔ)模型

本文并沒(méi)有按照傳統(tǒng)建模方法構(gòu)建模型結(jié)構(gòu)，而是根據(jù)國(guó)內(nèi)外關(guān)于標(biāo)準(zhǔn)樹(shù)高曲線的相關(guān)研究，選用了9 個(gè)常用方程［15-17］對(duì)杉木的樹(shù)高生長(zhǎng)進(jìn)行模擬(見(jiàn)表2)，這些方程都是前人建立起來(lái)的，且效果較好，可以直接進(jìn)行運(yùn)用。

2.2 非線性混合模型構(gòu)建

由于數(shù)據(jù)來(lái)自一個(gè)大范圍內(nèi)的不同樣地，有必要考慮樣地間的差異性，而包含隨機(jī)參數(shù)的混合模型可以模擬這種差異性。SAS 中NLMIXED［18］模塊可以提供相關(guān)的分析。

非線性混合模型的一般表達(dá)式為:

式中:yi和xi分別為i 樣地的因變量和自變量，且都為ni×1 維向量，Φi為r ×1 維的參數(shù)向量，具體到每一個(gè)樣地，r 為模型中參數(shù)的個(gè)數(shù)，f 是非線性方程，ei是ni×1 維的殘差向量。

表2 標(biāo)準(zhǔn)樹(shù)高曲線方程

混合模型的參數(shù)向量是由固定部分和隨機(jī)部分組成，其最主要的特點(diǎn)就是允許參數(shù)向量通過(guò)隨機(jī)效應(yīng)具體到每一個(gè)樣地。

式中:λ 為p ×1 維的固定效應(yīng)向量(p 為模型中固定參數(shù)的個(gè)數(shù))，bi為與i 樣地相關(guān)的q ×1 維隨機(jī)效應(yīng)向量(q 為模型中隨機(jī)參數(shù)的個(gè)數(shù))，Ai和Bi分別為r×p 維的固定效應(yīng)和r ×q 維的隨機(jī)效應(yīng)的設(shè)計(jì)矩陣，且具體到每一個(gè)樣地，其元素通常為0、1，或與固定效應(yīng)和隨機(jī)效應(yīng)相關(guān)的協(xié)方差值。

本文依據(jù)Fang and Bailey［19］提出的混合模型構(gòu)建方法:

(1)混合參數(shù)的確定。確定模型中的固定參數(shù)、隨機(jī)參數(shù)是構(gòu)建混合模型的基礎(chǔ)。根據(jù)Sharma and Parton 的研究，首先將基礎(chǔ)方程中所有的參數(shù)當(dāng)作隨機(jī)參數(shù)，如果無(wú)法收斂，相應(yīng)的減少隨機(jī)參數(shù)的個(gè)數(shù)來(lái)達(dá)到收斂，最后選擇收斂的模型為模擬結(jié)果，比較其負(fù)二倍對(duì)數(shù)似然值(-2Log Likelihood)、AIC(Akaike Information Criterion 赤池信息量準(zhǔn)則)和BIC(Bayesian Information Criterion 貝葉斯信息準(zhǔn)則)值，這3 個(gè)值越小，模擬效果越好，并進(jìn)行似然比檢驗(yàn)。

(2)確定樣地內(nèi)方差協(xié)方差結(jié)構(gòu)。一般數(shù)據(jù)存在殘差方差值的異方差性以及不同測(cè)量之間存在關(guān)聯(lián)性，為了避免這些問(wèn)題的出現(xiàn)，引入了包含關(guān)聯(lián)效應(yīng)和權(quán)重因子的樣地內(nèi)方差協(xié)方差矩陣Ri(λ，bi，p)，其一般形式為:

式中:σ2為誤差方差值，由模型的殘差方差值給定;Gi為ni×ni維的對(duì)角矩陣來(lái)解釋方差異質(zhì)性，其對(duì)角元素為相應(yīng)誤差項(xiàng)的標(biāo)準(zhǔn)差;Ini為ni×ni維的矩陣，來(lái)描述樣地i 不同觀測(cè)之間的關(guān)系。由于本文同一樣地并沒(méi)有多次測(cè)量，數(shù)據(jù)間不存在關(guān)聯(lián)性，無(wú)須考慮自相關(guān)性，因此，Ini為ni×ni維的單位矩陣。

本文試圖通過(guò)對(duì)殘差方差增加權(quán)重消除異方差。從自變量為胸徑的指數(shù)函數(shù)和冪函數(shù)中，根據(jù)負(fù)二倍對(duì)數(shù)似然值、AIC 和BIC 值，以及似然比檢驗(yàn)，確定一個(gè)效果最好的殘差方差模型。

式中:β 為待估參數(shù)，D 為胸徑。

(3)確定樣地間方差協(xié)方差結(jié)構(gòu)。隨機(jī)效應(yīng)的方差協(xié)方差矩陣(D)對(duì)于每一個(gè)樣地來(lái)說(shuō)都是一樣的，用來(lái)解釋樣地間的可變性，由于本研究中包含了兩個(gè)隨機(jī)參數(shù)，因此，矩陣D 為2 ×2 維的方差協(xié)方差矩陣，且表達(dá)式為:

式中:σ2u和σ2v分別為隨機(jī)參數(shù)u 和v 的方差，σuv是兩個(gè)隨機(jī)效應(yīng)的協(xié)方差。

2.3 隨機(jī)參數(shù)的估計(jì)

基礎(chǔ)模型在添加混合參數(shù)之后，在預(yù)估過(guò)程中要考慮以下兩種情況:

(1)如果某一個(gè)樣地只有林分因子變量和胸徑數(shù)據(jù)，那么在估測(cè)樹(shù)高時(shí)，就用只含有固定參數(shù)的模型，它代表了樹(shù)高變異的一般趨勢(shì)。

(2)如果某一個(gè)樣地除了林分因子變量和胸徑數(shù)據(jù)已知外，子樣本k 棵樹(shù)高數(shù)據(jù)已知，那么就可以預(yù)估樣地水平的隨機(jī)參數(shù)向量bi，且bi預(yù)測(cè)有以下表達(dá)式［20］:

式中:^D 為q×q 維的樣地之間的方差協(xié)方差矩陣，q為模型中隨機(jī)參數(shù)的個(gè)數(shù)為k ×k 維的樣地內(nèi)的方差協(xié)方差矩陣，為k ×1 維的殘差向量，為子樣本的樹(shù)高實(shí)測(cè)值與固定效應(yīng)模型的預(yù)測(cè)值之間的差值，為k×q 維矩陣且可以用來(lái)表示:

式中:Φi為r×1 維的參數(shù)向量;λ1，…，λq是預(yù)估固定效應(yīng)向量中混合系數(shù)中固定的部分，dik為i 樣地k 棵樹(shù)的胸徑值。

用這種方法預(yù)估出來(lái)的隨機(jī)參數(shù)和固定參數(shù)一起構(gòu)成了某樣地的參數(shù)值，因此，那些只測(cè)量了胸徑的樹(shù)木的樹(shù)高可以根據(jù)這一樣地的參數(shù)、胸徑以及林分因子變量求得。

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

本文選取決定系數(shù)R2和均方根差RMSE以及平均絕對(duì)殘差||和AIC 值比較其擬合效果。R2越接近1，RMSE和||越接近于0，其擬合效果越好。在模型參數(shù)相同時(shí)，AIC 值越小模型越優(yōu)。為避免模型參數(shù)過(guò)多，本文利用似然比檢驗(yàn):

式中:L1和L2分別為模型1 和模型2 的最大似然函數(shù)值，LRT服從自由度為k1-k2的χ2分布。給定可靠性α=0.05，當(dāng)LRT≥χ20.05(k1-k2)拒絕原假設(shè)，說(shuō)明這2 個(gè)模型差異顯著;反之，2 個(gè)模型差異不顯著，應(yīng)選擇含隨機(jī)效應(yīng)參數(shù)較少的模型。

式中:hij為i 樣地第j 棵樹(shù)樹(shù)高實(shí)際觀測(cè)值為i 樣地第j 棵樹(shù)樹(shù)高預(yù)測(cè)值，為所有樣地樹(shù)高平均值，m為樣地個(gè)數(shù)，ni為i 樣地內(nèi)樹(shù)木株數(shù)，n 為樣本總數(shù)。

3 結(jié)果與分析

3.1 基礎(chǔ)模型的選擇

用SAS 中PROC NLMIXED 模塊運(yùn)用表1中建模數(shù)據(jù)，對(duì)表2中9 個(gè)標(biāo)準(zhǔn)樹(shù)高曲線方程進(jìn)行擬合，擬合結(jié)果和參數(shù)的估計(jì)值(見(jiàn)表3)。由表3可以看出:除了方程M4和方程M9模擬精度明顯較低外，其它7 個(gè)方程擬合精度都較高。由于各方程擬合結(jié)果之間差異并不是很明顯，為了進(jìn)一步分析各個(gè)方程的擬合效果，根據(jù)參數(shù)個(gè)數(shù)和AIC 值，選出方程M1、M6和M8，分別代表了三個(gè)參數(shù)、四個(gè)參數(shù)、五個(gè)參數(shù)的方程，并按胸徑分為不同的徑階(2 cm 為一個(gè)徑階)，根據(jù)每一個(gè)徑階求算出樹(shù)高估計(jì)值的平均殘差以及相應(yīng)的標(biāo)準(zhǔn)差(見(jiàn)表4)。

從表中可以看出:方程M6對(duì)每個(gè)徑階擬合效果最好，和表3的結(jié)果基本一致。通過(guò)查表可得都明顯小于LRT，不同參數(shù)方程間差異顯著。因此選用方程M6作為混合模型的基礎(chǔ)模型。為了檢驗(yàn)方程M6是否有意義，利用F 檢驗(yàn)方法對(duì)模型進(jìn)行檢驗(yàn)。結(jié)果表明:模型F =30 713.2 ＞F0.95(5.159 1)=2.21，表明方程M6對(duì)于描述杉木樹(shù)高生長(zhǎng)具有顯著意義。

表3 各標(biāo)準(zhǔn)樹(shù)高曲線方程擬合結(jié)果

表4 各方程分別徑階的平均誤差和標(biāo)準(zhǔn)差

3.2 非線性混合模型構(gòu)建

3.2.1 混合參數(shù)的確定

在方程M6的基礎(chǔ)上，考慮樣地效應(yīng)，利用SAS中NLMIXED 模塊對(duì)20 塊杉木建模數(shù)據(jù)進(jìn)行模擬計(jì)算。并根據(jù)構(gòu)建非線性混合模型步驟，對(duì)不同情況進(jìn)行模擬，結(jié)果見(jiàn)表5。

從表中可以看出當(dāng)a1和a3同時(shí)作為混合參數(shù)時(shí)，3 個(gè)模型擬合指標(biāo)最小，通過(guò)比較LRT值與χ2值，差異顯著，拒絕原假設(shè)，不同參數(shù)方程之間差異顯著，故選擇a1和a3作為混合參數(shù)。這與Sharma and Parton 用同一方程模擬加拿大北方樹(shù)種標(biāo)準(zhǔn)樹(shù)高曲線時(shí)，所建立的混合效應(yīng)方程的結(jié)果相似。

表5 模型收斂的模擬結(jié)果

3.2.2 誤差項(xiàng)方差協(xié)方差(R)結(jié)構(gòu)

加權(quán)殘差方差模型和不加權(quán)(獨(dú)立等方差)時(shí)模型的評(píng)價(jià)指標(biāo)(見(jiàn)表6)?？紤]殘差加權(quán)的兩種模型與不加權(quán)模型差異顯著，表明兩種結(jié)構(gòu)都可以消除異方差，且冪函數(shù)的三個(gè)指標(biāo)最小，因此，選它作為模型的殘差方差模型。

表6 各殘差方差模型評(píng)價(jià)指標(biāo)

3.2.3 模型參數(shù)估計(jì)

運(yùn)用建模數(shù)據(jù)建立相應(yīng)方程，各參數(shù)估計(jì)值(見(jiàn)表7)，所以基于非線性混合模型的杉木標(biāo)準(zhǔn)樹(shù)高曲線方程表達(dá)式為:

表7 混合模型和固定模型的參數(shù)估計(jì)值

3.3 模型驗(yàn)證

根據(jù)建模結(jié)果，利用檢驗(yàn)數(shù)據(jù)，對(duì)比分析方程M6與方程11 的預(yù)測(cè)效果。固定模型的驗(yàn)證:直接將檢驗(yàn)數(shù)據(jù)帶入用建模數(shù)據(jù)建立的模型中。混合模型的驗(yàn)證:根據(jù)前人的研究，本文在每一塊樣地隨機(jī)抽取兩棵樹(shù)并實(shí)測(cè)其樹(shù)高，利用SAS 中PROC IML模塊，按照公式5 和公式6 計(jì)算其隨機(jī)參數(shù)值，再帶入到方程11 中，預(yù)估樹(shù)高值。兩種模型的預(yù)測(cè)效果(見(jiàn)表8)。

表8 兩種模型預(yù)測(cè)效果

從表8可以看出，混合模型的決定系數(shù)(R2=0.941 5)大于固定模型的決定系數(shù)(R2=0.918 3);均方根誤差(RMSE=1.361 8)比固定模型減少了15.4%;平均絕對(duì)殘差(|ˉE| =0.989 8)比固定模型減小了13.4%。這3 個(gè)評(píng)價(jià)指標(biāo)說(shuō)明混合模型的預(yù)測(cè)精度比固定模型的預(yù)測(cè)精度高。兩個(gè)模型的樹(shù)高預(yù)測(cè)值的殘差分布(見(jiàn)圖2)，由圖2可知，混合模型的殘差分布略優(yōu)于固定模型。為了進(jìn)一步分析兩個(gè)模型對(duì)不同胸徑大小的樹(shù)高預(yù)測(cè)效果，將檢驗(yàn)數(shù)據(jù)按胸徑分為不同的徑階(2 cm 為一個(gè)徑階)，并計(jì)算每個(gè)徑階內(nèi)的平均絕對(duì)殘差(見(jiàn)圖3a)，可以看出混合模型的平均絕對(duì)殘差不管胸徑大小基本都小于固定模型。隨機(jī)抽取了一塊樣地?cái)?shù)據(jù)，對(duì)比兩個(gè)模型的預(yù)測(cè)情況對(duì)比(見(jiàn)圖3b)，混合模型覆蓋的點(diǎn)明顯多于固定模型。這些都說(shuō)明基于混合模型的標(biāo)準(zhǔn)樹(shù)高曲線優(yōu)于普通的標(biāo)準(zhǔn)樹(shù)高曲線。

圖2 兩種模型樹(shù)高預(yù)測(cè)值的殘差分布

4 結(jié)論與討論

本研究中選用的9 個(gè)標(biāo)準(zhǔn)樹(shù)高曲線方程中，其中7 個(gè)模擬效果都較好，這是由于方程中的那些林分因子變量都很好的反映了林分的實(shí)際生長(zhǎng)狀況。表現(xiàn)突出的方程M6考慮了優(yōu)勢(shì)樹(shù)高和平均胸徑，分別代表了林分的立地條件和密度狀況，在實(shí)際運(yùn)用中也只需測(cè)量胸徑和優(yōu)勢(shì)樹(shù)樹(shù)高，具有較強(qiáng)的實(shí)用性。

通過(guò)模型的驗(yàn)證表明，利用基礎(chǔ)方程建立的非線性混合模型與固定模型相比，考慮樣地隨機(jī)影響時(shí)，混合模型的預(yù)測(cè)效果明顯優(yōu)于固定模型，能夠更好地描述當(dāng)?shù)厣寄緲?shù)高生長(zhǎng)趨勢(shì)，為樹(shù)高預(yù)測(cè)提供了依據(jù)。且利用負(fù)二倍對(duì)數(shù)似然值和AIC、BIC 評(píng)價(jià)非線性混合模型的效果表明，在考慮樣地效應(yīng)時(shí)，a1、a3同時(shí)作為混合參數(shù)時(shí)模擬效果最好，這兩個(gè)參數(shù)分別決定其漸近線和生長(zhǎng)率系數(shù)。

圖3 兩種模型預(yù)測(cè)殘差及效果圖

本研究是基于單水平研究非線性混合模型對(duì)樹(shù)高生長(zhǎng)的模擬，并沒(méi)有研究嵌套多水平對(duì)模型精度的影響，但是，為后續(xù)的研究以及混合模型在林業(yè)上的廣泛應(yīng)用提供了依據(jù)和參考。

［1］ 丁貴杰.貴州杉木人工林標(biāo)準(zhǔn)樹(shù)高曲線模型［J］.貴州農(nóng)學(xué)院學(xué)報(bào)，1996，15(4):16 -21.

［2］ 王明亮，唐守正.標(biāo)準(zhǔn)樹(shù)高曲線的研制［J］.林業(yè)科學(xué)研究，1997，10(3):259 -264.

［3］ Adame P，Del Rio M，Canellas I.A mixed nonlinear height-diameter model for pyrenean oak (Quercus pyrenaica Willd.)［J］.Forest ecology and management，2008，256(1/2):88 -98.

［4］ Paulo J A，Tomé J，Tomé M.Nonlinear fixed and random generalized height-diameter models for portuguese cork oak stands［J］.Annals of Forest Science，2011，68(2):295 -309.

［5］ 李春明，李利學(xué).基于非線性混合模型的栓皮櫟樹(shù)高與胸徑關(guān)系研究［J］.北京林業(yè)大學(xué)學(xué)報(bào)，2009，31(4):7 -12.

［6］ Sharma M，Parton J.Height-diameter equations for boreal tree species in ontario using a mixed-effects modeling approach［J］.Forest Ecology and Management，2007，249(3):187 -198.

［7］ Lappi J.A longitudinal analysis of height/diameter curves［J］.Forest science，1997，43(4):555 -570.

［8］ Calama R，Montero G.Interregional nonlinear height-diameter model with random coefficients for stone pine in spain［J］.Canadian Journal of Forest Research，2004，34(1):150 -163.

［9］ Vargas-Larreta B，Castedo-Dorado F，Alvarez-Gonzalez J G，et al.A generalized height-diameter model with random coefficients for uneven-aged stands in El Salto，Durango (Mexico)［J］.Forestry，2009，82(4):445 -462.

［10］ 胥輝，全宏波，王斌.思茅松標(biāo)準(zhǔn)樹(shù)高曲線的研究［J］.西南林學(xué)院學(xué)報(bào)，2000，20(2):74 -77.

［11］ 丁貴杰.馬尾松人工林標(biāo)準(zhǔn)樹(shù)高曲線模型的研究［J］.浙江林學(xué)院學(xué)報(bào)，1997，14(3):15 -20.

［12］ 孫圓，王萬(wàn)江.江蘇省楊樹(shù)樹(shù)高曲線模型的研制［J］.林業(yè)科技開(kāi)發(fā)，2005，19(5):31 -34.

［13］ 趙俊卉，亢新剛，張慧東，等.長(zhǎng)白山3 個(gè)主要針葉樹(shù)種的標(biāo)準(zhǔn)樹(shù)高曲線［J］.林業(yè)科學(xué)，2010，46(10):191 -194.

［14］ 魏曉慧，孫玉軍.基于Richards 方程的杉木樹(shù)高生長(zhǎng)模型［J］.浙江農(nóng)林大學(xué)學(xué)報(bào)，2012，29(5):661 -666.

［15］ Schr?der J，González J G á.Comparing the performance of generalized diameter-height equations for maritime pine in northwestern spain［J］.Forstwissenschaftliches Centralblatt vereinigt mit Tharandter forstliches Jahrbuch，2001，120(1/6):18 -23.

［16］ Temesgen H，Gadow K V.Generalized height-diameter modelsan application for major tree species in complex stands of interior british columbia［J］.European Journal of Forest Research，2004，123(1):45 -51.

［17］ Crecente-Campo F，Tomé M，Soares P，et al.A generalized nonlinear mixed-effects height-diameter model for Eucalyptus globulus L.in northwestern spain［J］.Forest Ecology and Management，2010，259(5):943 -952.

［18］ Jones A，Huddleston E.SAS/STAT 9.2 user’s guide［M］.SAS Institute Inc Cary NC，2009.

［19］ Grégoire T G，Schabenberger O，Barrett J P.Linear modelling of irregularly spaced，unbalanced，longitudinal data from permanent-pBot measurements［J］.Canadian Journal of Forest Research，1995，25(1):137 -156.

［20］ Lindstrom M L，Bates D M.Nonlinear mixed effects models for repeated measures data［J］.Biometrics，1990，46(3):673 -687.