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

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

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

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

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

1 研究區(qū)概況

研究區(qū)為福建省福州市，位于福建省中部東端、閩江下游，全市面積11 968 km2。地處北緯25°15′～26°39′，東經(jīng)118°08′～120°31′，海拔600～1 000 m。該市屬于典型的亞熱帶季風(fēng)氣候，溫暖濕潤，雨量充沛，陽光充足，無霜期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)地中隨機抽取217株樣木，其所在標(biāo)準(zhǔn)地基本信息詳見表1。伐倒樣木之后，用皮尺 (精度為0.01 m) 測量全樹高和冠長，然后再以1 m為間隔，于樹干上0.3、0.5、1.3、1.5、2.5 m等位置處，用鋼卷尺 (精度為0.1 cm) 依次讀取帶皮直徑和去皮直徑 (剝除樹皮后讀取)，并計算樹皮厚度。

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

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

2.2 模型形式

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

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

表3 模型形式Table 3 The models form

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

2.3 模型質(zhì)量評價指標(biāo)

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

(1)

(2)

(3)

AIC=-2lnζ+2p

(4)

(5)

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

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

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

(6)

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

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

2.5 模型檢驗方法

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

3 結(jié)果與分析

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

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

3.2 模型擬合質(zhì)量評價

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

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

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

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

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

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

3.3 模型檢驗

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

表6 模型配對t檢驗Table 6 t test for model paired

圖14個模型殘差
Fig.1 Residuals of 4 selected models

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

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

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

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

4 結(jié)論與討論

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

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

王曉林等[9]對落葉松的研究中發(fā)現(xiàn)，樹皮厚度與冠長和冠長率存在著一定相關(guān)性，但無法確保其對樹皮厚度有顯著影響。因此本研究在任意高度處樹皮厚度的模型中考慮了冠長和冠長率2個自變量，在模型8中新加入了自變量冠長，模型9和10中新加入了自變量冠長率，結(jié)果發(fā)現(xiàn)模型參數(shù)檢驗均顯著并且決定系數(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。

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

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

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