嚴(yán)銘海 方靜儀 李睿宇 江希鈿

(福建農(nóng)林大學(xué)林學(xué)院，福建 福州 350002)

樹皮厚度受每個(gè)樹種的遺傳因子控制，導(dǎo)致不同樹種的樹皮厚度往往有差異[1]。Zeibig-Kichas等[2]對(duì)加利福尼亞針闊混交林中8種針葉樹的樹皮厚度進(jìn)行了研究，發(fā)現(xiàn)不同針葉樹種之間差異大，8種針葉樹中樹皮厚度最大的是北美翠柏 (Calocedrusdecurrens)，最小的是美國(guó)黑松 (Pinuscontorta)。樹皮厚度在樹種間變化很大，但一般隨直徑的增加而增加[3]。有研究結(jié)果顯示，直徑每增加10 cm，樹皮厚度約增加1.2 mm[4]。同一樹干不同高度對(duì)應(yīng)的直徑往往不同，導(dǎo)致樹皮厚度也隨之變化。Laasasenaho等[5]也認(rèn)為樹皮厚度不僅取決于樹種，同時(shí)也和樹干的不同高度位置密切相關(guān)。還有研究表明，針葉樹比常綠闊葉樹種具有更大的樹皮厚度[6]。

國(guó)內(nèi)外學(xué)者對(duì)樹皮厚度進(jìn)行研究，得出了許多計(jì)算樹皮厚度的數(shù)學(xué)模型。其中，胸徑通常作為一個(gè)很重要的解釋變量，此外樹高和年齡等變量的加入，也取得了較好的擬合效果[7-10]。但同一樹干的不同位置，樹皮厚度往往不一樣，因此任意樹干高度的樹皮厚度方程被建立[11-13]。除此之外，相對(duì)樹皮厚度的方程也被一些學(xué)者用來計(jì)算樹皮厚度[12-14]。較常見的還有去皮直徑模型，在求不同高度的樹皮厚度時(shí)也很方便[15-17]。

濕地松 (Pinuselliottii) 原產(chǎn)于美國(guó)東南部，我國(guó)從20世紀(jì)30年代開始引種，因其具有速生、適應(yīng)性強(qiáng)、木材用途廣和松脂產(chǎn)量高等優(yōu)良特性，現(xiàn)已成為亞熱帶地區(qū)大面積造林樹種，是我國(guó)引種的優(yōu)良速生樹種之一[18-19]。濕地松的樹皮有不同于松脂的濃厚香味，其樹皮提取物不僅可以用于制備改性酚醛樹脂膠和膠合劑，還具有治療平喘的藥理作用[20-21]。目前我國(guó)濕地松人工林種植面積超過300萬hm2[22]，木材和樹皮資源量巨大。研究濕地松樹皮厚度模型對(duì)于木材和樹皮資源的生產(chǎn)利用都非常重要，但目前還缺少合適的模型來預(yù)估濕地松樹皮厚度?；诖?，本研究通過收集前人構(gòu)建的模型和自建模型，應(yīng)用前期濕地松樣木獲取的帶皮直徑 (dob)、去皮直徑 (dib)、樹皮厚度等數(shù)據(jù)對(duì)模型進(jìn)行擬合，對(duì)比分析不同模型的擬合效果，并篩選出各類樹皮厚度預(yù)測(cè)的合理模型，從而為濕地松材積、出材率和樹皮生長(zhǎng)量的估測(cè)提供參考。

1 研究區(qū)概況

研究區(qū)為福建省福州市，位于福建省中部東端、閩江下游，全市面積11 968 km2。地處北緯25°15′～26°39′，東經(jīng)118°08′～120°31′，海拔600～1 000 m。該市屬于典型的亞熱帶季風(fēng)氣候，溫暖濕潤(rùn)，雨量充沛，陽光充足，無霜期326 d，年均日照數(shù)1 700～1 980 h，年均溫20～25 ℃，年均降雨量900～2 100 mm。濕地松是當(dāng)?shù)刂匾脑炝謽浞N。

2 研究方法

2.1 數(shù)據(jù)采集

結(jié)合伐區(qū)木材生產(chǎn)，在福州市濕地松人工林中設(shè)置標(biāo)準(zhǔn)地共87塊，考慮不同的海拔 (h)、坡度 (s)、郁閉度 (c)、林分年齡 (t)、胸徑 (DBH)、樹高 (H) 等因素，從標(biāo)準(zhǔn)地中隨機(jī)抽取217株樣木，其所在標(biāo)準(zhǔn)地基本信息詳見表1。伐倒樣木之后，用皮尺 (精度為0.01 m) 測(cè)量全樹高和冠長(zhǎng)，然后再以1 m為間隔，于樹干上0.3、0.5、1.3、1.5、2.5 m等位置處，用鋼卷尺 (精度為0.1 cm) 依次讀取帶皮直徑和去皮直徑 (剝除樹皮后讀取)，并計(jì)算樹皮厚度。

表1 樣木所在標(biāo)準(zhǔn)地基本信息Table 1 Information of sample plots

由于需要考慮模型的擬合和驗(yàn)證，故將217株樣木按75%和25%的比例隨機(jī)分為兩部分，一部分用于建模 (樣木163 株)，一部分用于模型檢驗(yàn) (樣木54 株)，樣木數(shù)據(jù)的描述性統(tǒng)計(jì)結(jié)果見表2。模型變量包括林分年齡、樹高、胸徑、胸高處樹皮厚度 (BBT)、任意樹高處樹皮率 (z)、相對(duì)樹皮厚度 (RBT)、樹冠率 (CR)、樹冠長(zhǎng) (CL)，其中任意樹高處樹皮率為該處去皮直徑與該處帶皮直徑之比。

2.2 模型形式

通過查閱資料，參考前人建立的14個(gè)模型，再加上自行構(gòu)建的4個(gè)模型，共獲取18個(gè)模型用于本研究，具體模型見表3。樹皮模型的分類采用唐誠(chéng)等[13]的研究結(jié)果，依據(jù)模型因變量將其歸納為4種。

表2 濕地松人工林樹皮厚度建模和檢驗(yàn)數(shù)據(jù)Table 2 Fitting and validating data for bark thickness models of P.elliottii plantations

表3 模型形式Table 3 The models form

注：各模型中，BT表示任意高度處樹皮厚度；RD表示相對(duì)直徑；RH表示相對(duì)樹高；a表示常數(shù)；B1、B2、B3表示模型參數(shù)。

2.3 模型質(zhì)量評(píng)價(jià)指標(biāo)

4類樹皮厚度模型均采用最小二乘法進(jìn)行擬合，并對(duì)模型參數(shù)進(jìn)行顯著性檢驗(yàn)，滿足所有參數(shù)均顯著的模型，進(jìn)一步計(jì)算赤池信息準(zhǔn)則 (AIC)、交叉驗(yàn)證誤差 (CVE)、偏差 (B)、絕對(duì)偏差 (AB)、決定系數(shù) (R2) 作為模型質(zhì)量評(píng)價(jià)的5個(gè)指標(biāo)，其計(jì)算方法分別見公式 (1)～(5)。

(1)

(2)

(3)

AIC=-2lnζ+2p

(4)

(5)

在計(jì)算交叉驗(yàn)證誤差時(shí)，將建模數(shù)據(jù)分為沒有重疊的k份 (k=5)[23]。用k-1份去擬合模型，用第k份進(jìn)行檢驗(yàn)，這個(gè)過程重復(fù)k次，直到每一部分都被抽出做一次模型檢驗(yàn)，做k-1次模型擬合，計(jì)算k次均方誤差 (MSE) 的平均值即得到交叉驗(yàn)證誤差值。在R軟件中，CVE通過cv.glm函數(shù)計(jì)算得到。

2.4 模型優(yōu)選方法

應(yīng)用相對(duì)排序法[24]進(jìn)行模型優(yōu)選，分別計(jì)算各模型AIC、B、AB、CVE和R2的相對(duì)排序值，其計(jì)算方法見公式 (6)。

(6)

式中：Ranki為模型i的相對(duì)排序值 (i=1, 2, …,w)，其值越小說明模型擬合效果越好；w為參與排序的模型數(shù)目;Si為模型i的統(tǒng)計(jì)指標(biāo)值 (分別為AIC、B、AB、CVE和R2)；Smin和Smax分別為Si的最小值和最大值。

偏差均值和AIC都存在有負(fù)值，但偏差需取絕對(duì)值來計(jì)算排序值，而AIC是無論正負(fù)，越小越好，故可直接參與排序。對(duì)于R2則是越大越好，為了統(tǒng)一越小越好的標(biāo)準(zhǔn)，將1-R2作為排序的統(tǒng)計(jì)指標(biāo)值。最后求出5個(gè)統(tǒng)計(jì)指標(biāo)的平均相對(duì)排序值，通過比較其大小可判斷模型質(zhì)量的優(yōu)劣，并篩選出最優(yōu)模型。

2.5 模型檢驗(yàn)方法

對(duì)最優(yōu)模型進(jìn)行配對(duì)t檢驗(yàn)來判斷預(yù)估值和實(shí)際值之間差異是否顯著，若差異不顯著，則說明模型有效性良好。進(jìn)一步通過方差擴(kuò)大因子法對(duì)其進(jìn)行多重共線性的診斷，再通過懷特檢測(cè)其是否存在異方差，如果有明顯異方差的模型，則可通過變量變換進(jìn)行修正[25]。

3 結(jié)果與分析

3.1 模型擬合參數(shù)檢驗(yàn)

以217株樣木的相關(guān)數(shù)據(jù)為依托，對(duì)18個(gè)模型進(jìn)行擬合，其參數(shù)顯著性檢驗(yàn)結(jié)果 (表4) 表明：模型2、5和11均有一個(gè)參數(shù)預(yù)估值與零差異不顯著，將不參與進(jìn)一步分析。分析結(jié)果可以看出模型2、5、11中不顯著的都是H或lnH的參數(shù)，這與唐誠(chéng)等[13]的研究結(jié)果一致。

3.2 模型擬合質(zhì)量評(píng)價(jià)

模型擬合質(zhì)量好壞很大程度上影響著預(yù)估的精度，模型統(tǒng)計(jì)指標(biāo)計(jì)算及排序見表5。

表4 模型參數(shù)及其顯著性檢驗(yàn)Table 4 Model parameters and its significance test

注：*和ns分別表示模型參數(shù)預(yù)估值在5%水平與零差異顯著和不顯著。

表5 模型統(tǒng)計(jì)指標(biāo)計(jì)算及排序Table 5 Model statistical index calculation and ranking

模型B值相對(duì)排序值A(chǔ)B值相對(duì)排序值CVE值相對(duì)排序值A(chǔ)IC值相對(duì)排序值1-R2值相對(duì)排序值相對(duì)排序值均值10.000 1.001 0.437 11.408 0.331 1.005 249.161 4.402 0.640 14.000 6.363 30.094 12.293 0.442 11.574 0.208 1.003 186.095 4.379 0.552 12.211 8.292 40.091 11.989 0.443 11.602 0.204 1.003 183.141 4.378 0.534 11.836 8.161 60.108 14.000 0.285 6.445 0.126 1.002 2055.314 5.078 0.375 8.591 7.023 70.093 12.194 0.281 6.320 0.107 1.001 1619.298 4.915 0.318 7.436 6.373 8-0.0001.000 0.229 4.636 0.093 1.001 1246.407 4.775 0.277 6.596 3.60290.000 1.000 0.237 4.884 0.096 1.001 1338.884 4.810 0.286 6.794 3.698 10-0.0001.000 0.227 4.567 0.091 1.001 1174.355 4.748 0.269 6.445 3.552 120.000 1.000 0.108 0.661 0.025 1.000 -2315.266 3.443 0.226 5.563 2.333 130.0001.000 0.118 1.000 0.027 1.000 -2082.532 3.530 0.246 5.981 2.502 140.000 1.000 0.125 1.222 0.032 1.000 -1659.435 3.688 0.288 6.835 2.749 150.021 3.524 0.516 14.000 0.505 1.007 5781.896 6.472 0.002 1.000 5.201 16-0.0001.000 0.508 13.718 0.504 1.007 5771.969 6.468 0.011 1.181 4.675 17-0.043-4.159 0.509 13.773 928.422 14.000 25908.153 14.000 0.008 1.122 7.747 180.003 1.414 0.508 13.736 0.002 1.000 -8846.419 1.000 0.007 1.096 3.649

由表5可知，14個(gè)模型中6、7、3和4的偏差值最大，接近于0.1，模型17和15的偏差值為-0.043和0.021，其余9個(gè)模型的偏差值均接近于0。擬合去皮直徑的4個(gè)模型絕對(duì)偏差最大，約為0.51；其次是擬合胸高處樹皮厚度的3個(gè)模型，約為0.44；擬合任意樹高處樹皮厚度的模型中，模型6和7的絕對(duì)偏差最大，其余的模型8、9、10約為0.23；擬合相對(duì)樹皮厚度的模型絕對(duì)偏差最小。模型18的CVE值最小，模型17的CVE最大，剩余模型中CVE最大的是去皮直徑模型；相對(duì)樹皮厚度CVE模型最小，任意樹高處樹皮厚度模型的CVE值在0.1附近，胸高處樹皮厚度CVE值約在0.2～0.3。AIC的排序情況和CVE基本一致。R2以模型15、16、17、18為最高，均接近1。而模型1、3、4的R2則在0.5以下。

從5項(xiàng)指標(biāo)的綜合排序值可以看出，模型3、4、6、7、17的排序值較高，模型12、13、14的綜合排序值較低，其他模型居中。同形式模型進(jìn)行比較，胸高處樹皮厚度的3個(gè)模型中，模型1擬合質(zhì)量最好，模型4最差。任意樹高處樹皮厚度的5個(gè)模型中，模型10擬合效果最好，模型6最差。相對(duì)樹皮厚度的3個(gè)模型中，模型12的擬合質(zhì)量最好，模型15的擬合質(zhì)量最差。去皮直徑的4個(gè)模型中，模型18的擬合質(zhì)量最優(yōu)，模型17最差。由此篩選出較優(yōu)模型為1、10、12、18。

3.3 模型檢驗(yàn)

運(yùn)用54株樣木的相關(guān)數(shù)據(jù)對(duì)最優(yōu)模型進(jìn)行配對(duì)t檢驗(yàn)，模型1、10、12、18的預(yù)估值和實(shí)際值之間沒有顯著差異 (表6)。進(jìn)一步對(duì)其進(jìn)行多重共線性診斷和異方差分析 (圖1)。

表6 模型配對(duì)t檢驗(yàn)Table 6 t test for model paired

圖14個(gè)模型殘差
Fig.1 Residuals of 4 selected models

4個(gè)最優(yōu)模型中的模型1、12、18均只有1個(gè)自變量，故不需要考慮共線性問題。而模型10則用方差擴(kuò)大因子法分析其共線性，其3個(gè)自變量的VIF值分別為1.079、2.414、2.361，均小于10，認(rèn)為模型10的各自變量不存在嚴(yán)重的多重共線性。

從4個(gè)模型的殘差圖 (圖1) 可以看出，模型1和12的殘差分布均勻且隨機(jī)，經(jīng)過懷特檢驗(yàn)表明模型1、12不存在異方差，模型10和18則有異方差問題，故通過取對(duì)數(shù)后回歸。變換后模型的殘差圖 (圖2) 和評(píng)價(jià)指標(biāo) (表7) 均有改善。

圖2異方差修正后模型殘差
Fig.2 Residuals against fitted values for 5 models after correction of heteroskedasticity

表7 原模型與變換后模型統(tǒng)計(jì)指標(biāo)對(duì)比Table 7 Statistics indices comparison between the original models and changed models

4 結(jié)論與討論

本研究通過5個(gè)評(píng)價(jià)指標(biāo)篩選出模型1、10、12和18，分別用于擬合胸高處樹皮厚度、任意樹高處樹皮厚度、相對(duì)樹皮厚度、任意高度去皮直徑，4個(gè)模型擬合精度高，不存在多重共線性問題，且模型沒有異方差或經(jīng)過變量變換得到較好的修正。擬合模型所需的數(shù)據(jù)除了冠長(zhǎng)率之外，其余數(shù)據(jù)在林業(yè)生產(chǎn)實(shí)踐中具有易獲取性。在樹皮厚度模型中加入冠長(zhǎng)率，對(duì)于濕地松而言，擬合效果較好。

胸徑是擬合胸徑處樹皮厚度模型的重要自變量，此外，再通過增加一些自變量也可以提高模型的決定系數(shù)。本研究中的模型2、4、5增加了自變量樹高或年齡，然而因?yàn)樾馗咛帢淦ず穸扰c樹高不相關(guān)，導(dǎo)致模型2、5的參數(shù)檢驗(yàn)不顯著。模型3、4是在模型1的基礎(chǔ)上通過變量的對(duì)數(shù)變換后得到，可以看出其決定系數(shù)明顯提高，但模型3的實(shí)際值和預(yù)測(cè)值差異顯著。剩余的模型中，從模型質(zhì)量指標(biāo)綜合排序來看，模型1的擬合效果優(yōu)于模型4。

王曉林等[9]對(duì)落葉松的研究中發(fā)現(xiàn)，樹皮厚度與冠長(zhǎng)和冠長(zhǎng)率存在著一定相關(guān)性，但無法確保其對(duì)樹皮厚度有顯著影響。因此本研究在任意高度處樹皮厚度的模型中考慮了冠長(zhǎng)和冠長(zhǎng)率2個(gè)自變量，在模型8中新加入了自變量冠長(zhǎng)，模型9和10中新加入了自變量冠長(zhǎng)率，結(jié)果發(fā)現(xiàn)模型參數(shù)檢驗(yàn)均顯著并且決定系數(shù)均提高至0.7以上。任意高度的帶皮直徑是模擬任意高度處樹皮厚度的重要解釋變量，模型7、8、9、10在模型6的基礎(chǔ)上陸續(xù)加入RH、CL、CR和RD等自變量，擬合效果均有改善。經(jīng)過比較模型質(zhì)量的綜合排序值，認(rèn)為模型10擬合效果優(yōu)于模型6、7、8和9。

鑒于冠長(zhǎng)率在模型10中取得良好的擬合效果，因此嘗試在模型14的自變量中增加冠長(zhǎng)率，與其他相對(duì)樹皮厚度模型相比，模型14雖然參數(shù)檢驗(yàn)顯著，但相比較而言決定系數(shù)較低。模型11也因?yàn)樽宰兞坑袠涓?，?dǎo)致模型參數(shù)與零相比不顯著。在擬合相對(duì)樹皮厚度時(shí)，有學(xué)者認(rèn)為模型12中的B2會(huì)因樹種差異而不固定[16,25]。Johnson等[14]對(duì)輻射松 (Pinusradiata) 的研究發(fā)現(xiàn)，B2賦值4.0時(shí)擬合效果良好；唐誠(chéng)等[13]對(duì)西南樺 (Betulaalnoides) 的研究表明，B2賦值約1.0時(shí)效果最好；對(duì)于濕地松而言，B2約1.7時(shí)擬合效果最佳。比較同類模型的綜合排序值，模型 12的擬合效果略優(yōu)于模型13、14。

任意樹高的去皮直徑在森林經(jīng)營(yíng)和木材加工過程中非常重要，Marshall等[26]研究發(fā)現(xiàn)去皮直徑的預(yù)測(cè)誤差會(huì)導(dǎo)致高達(dá)11%的木材價(jià)值損失,同時(shí)去皮直徑也是木材企業(yè)生產(chǎn)時(shí)關(guān)心的問題。去皮直徑主要通過帶皮直徑來擬合，本研究中用來擬合去皮直徑模型的R2都接近1。模型15的B1常稱為樹皮因子，模型17是在模型15的基礎(chǔ)上對(duì)變量進(jìn)行平方變換來的。模型16和模型18分別在模型15的基礎(chǔ)上增加常數(shù)項(xiàng)和參數(shù)B2。綜合比較各統(tǒng)計(jì)指標(biāo)值，認(rèn)為模型18的擬合效果最佳。

樹皮厚度除受直徑、冠長(zhǎng)、冠幅、樹高、年齡等因子的影響外，還可能受生長(zhǎng)速率的影響。有學(xué)者的研究表明，不同樹種的生長(zhǎng)速率往往不同，生長(zhǎng)較快的樹種一般樹皮略薄[2]。St?ngle等[27]對(duì)德國(guó)西南部云杉 (Piceaasperata) 樹皮厚度估測(cè)模型的研究中，發(fā)現(xiàn)立地條件對(duì)樹皮厚度的影響很大，立地質(zhì)量較好，相對(duì)應(yīng)的樹皮厚度較小。立地質(zhì)量有時(shí)也會(huì)影響樹木的生長(zhǎng)速率，同時(shí)也是影響樹皮厚度的重要因素。但本研究中未考慮立地因子，在往后的研究可以嘗試加入立地因子來提高模型精度。目前，樹皮厚度方面的林業(yè)數(shù)表很少，如果有條件可基于樹皮模型的研究構(gòu)建相對(duì)應(yīng)的林業(yè)數(shù)表，也是對(duì)林業(yè)數(shù)表的一種補(bǔ)充。此外，樹皮厚度在森林動(dòng)態(tài)模擬中的其他方面應(yīng)用，還有待深入研究。