基于啞變量的木荷人工林樹皮厚度模型研究

張偉志 吳宏煒 黃光燦 田意 莊崇洋 江希鈿

摘? 要：在木荷人工林中，基于啞變量模型法構建樹皮厚度模型，以提高預測模型參數的穩定性，為木荷人工林出材率計算和經濟價值評估方面提供參考依據。基于福建省南平市213塊木荷人工林固定樣地的調查數據，選取19個包含胸高處、任意高度處、相對樹皮厚度和去皮直徑的理論方程作為構建木荷人工林樹皮厚度的基礎模型，使用R軟件進行模型擬合分析，運用決定系數（R2）、均方根誤差（Erms）、和方差（Ess）、赤池信息準則（AIC）、貝葉斯準則（BIC）等模型評價指標，結合含熵值的TOPSIS法篩選出擬合度較高的4個基礎模型（M1、M10、M12、M16），進一步構建含齡組和立地質量啞變量的木荷人工林樹皮厚度模型。用于建模的胸高處樣本數為130個，任意高度處、相對樹皮厚度和去皮直徑樣本數為2386個，用于檢驗模型擬合效果相對應的樣本數分別為55個和1013個。結果表明：含齡組啞變量胸高處、任意高度處、相對樹皮厚度模型（M20、M22、M25）的R2分別為0.9769、0.9214、0.9111，比基礎模型（M1、M10、M12）的R2（分別為0.6981、0.5540、0.5056）提高了39.9%、66.3%、80.2%，有了明顯提升;含齡組啞變量去皮直徑模型M27的AIC為21.62，BIC為76.44，比其基礎模型M16分別降低了95.8%、85.46%。通過含熵值的TOPSIS法對模型指標綜合評價，并經模型配對t檢驗，這4個含齡組啞變量模型（M20、M22、M25、M27）擬合效果更佳，適合于福建木荷人工林樹皮厚度預測。

關鍵詞：樹皮厚度;啞變量;齡組;立地質量;木荷

中圖分類號：S758;S79? ? ? 文獻標識碼：A

Abstract： In the artificial forest of Schima superba， a bark thickness model was constructed based on the dummy variable model method， in order to improve the stability of the prediction model parameters， and provide reference for the calculation of the timber yield and the evaluation of the economic value of the artificial forest of S. superba. Based on the survey data of 213 sample plots of S. superba plantation in Nanping City， Fujian Province， 19 theoretical equations including breast height， any height， relative bark thickness and peeling diameter were selected as the basic models for building the bark thickness of S. superba plantation. Based on R software model fitting analysis， determination coefficient （R2）， root mean square error （Erms）， the sum of squares due to error （Ess）， akaike information criterion （AIC） and bayesian information criterion （BIC） models as the evaluation index， and the entropy TOPSIS method to screen high fitting degree of four basic models （M1， M10， M12， M16） model， an age group and the site quality of S. superba plantation bark thickness of dummy variable model was further built. The number of samples at the breast height used for modeling was 130， the number of samples at any height， the relative bark thickness， and the peeling diameter was 2386， and the corresponding numbers for testing the model fitting effect was 55 and 1013 respectively. The results showed that R2 of the model with age dummy variable at breast height， any height and relative bark thickness （M20， M22 and M25） was 0.9769， 0.9214 and 0.9111， respectively， which was 39.9%， 66.3% and 80.2% higher than that of the basic model （M1， M10 and M12）， and AIC and BIC of M27 with age dummy variable was 21.62 and 76.44， respectively， it was 95.8% and 85.46% lower than that of the basic model M16. Through the comprehensive evaluation of the model indexes by the TOPSIS method with entropy value and t-test of model pairing， the fitting effect of these four age group dummy variable models （M20， M22， M25， M27） was better， which is suitable for the prediction of bark thickness of S. superba plantation in Fujian.

Keywords： bark thickness; dummy variable; age groups; site quality; Schima superba

DOI： 10.3969/j.issn.1000-2561.2020.12.027

樹皮不僅對樹木自身成長起到保護作用，還可成為生產及藥用材料等重要生物資源，若能合理對樹皮進行基質化利用還可進一步提高其經濟價值[1]。樹皮厚度是指樹干某處帶皮直徑（diameter outside bark， DOB）與該處去皮直徑（diameter inside bark， DIB）差的一半[2]。樹皮大約占到整個木材體積的12%～20%，不同樹種的樹皮總量變化較大[3]，同一樹種在不同位置的樹皮厚度也不同[4]。因此，樹皮厚度合理測算將直接影響到樹皮蓄積量準確估算。

國外學者對樹皮厚度的研究，早期主要集中在胸徑處，如Grosenhbaugh[5]建立比率方程式，通過預測胸徑高度處去皮直徑來測算樹皮厚度。Wingerd等[6]利用多元回歸分析方法，構建美國阿巴拉契亞山脈7個樹種的去皮直徑回歸模型，Maguire等[7]則利用削度方程來測算樹干不同高度處去皮直徑，這些方法都較大提高了樹皮厚度計算的準確度。國內研究者通過加入胸徑、樹高、任意樹高、年齡等解釋變量，對白樺、落葉松、西南樺、濕地松、厚樸、云杉等樹種的樹皮厚度進行模擬，發現不同變量對樹皮厚度影響也不同[8-13]。但目前研究還較少考慮到立地質量、齡組等因子對樹皮厚度影響。賈煒瑋等[14]建立含齡組、區域啞變量碳儲量預測模型系統，結果顯示含有啞變量的模型精度較高。朱光玉等[15]則構建含林分和立地類型啞變量的林分斷面積生長模型，通過比較發現啞變量模型精度均優于基礎模型。所以筆者將通過使用啞變量的方法將這些因子納入到研究中。

木荷（Schima superba），也稱為荷樹、荷木，屬于山茶科（Theaceae）常綠喬木[16-17]。木荷不僅是珍貴高級用材樹種[18]，也是重要的防火樹種[19]。近年來，福建和廣東兩省種植了大量的木荷人工林[20]，已有學者在木荷培育[21]、病蟲害防治[22-23]等方面有研究，但還較少對木荷樹皮厚度進行探討。為此，本研究以福建省南平市木荷人工林為研究對象，通過對樣木及樣地調查獲取模型所需因子，分別用胸高處樹皮厚度（bark thickness at breast height， BBT）、任意高度處樹皮厚度（bark thickness， BT）、相對樹皮厚度（relative bark thickness， RBT）和去皮直徑（diameter inside bark， DIB）4類基礎模型[10]，擬合篩選出各自最優模型后，加入包含齡組和立地質量的啞變量，通過含熵值的TOPSIS綜合分析法選出相容且統一的木荷樹皮厚度模型，為計算木荷材積、出材量及生產經營措施制訂提供參考依據，對木荷人工林經濟價值精準評估具有重要現實意義。

1? 材料與方法

1.1? 材料

1.1.1? 研究區概況? 研究區為福建省南平市，位于福建北部，武夷山脈北段東南側，全市面積約2.63萬km2。南平地處北緯26°15～28°19，東經117°00～119°17之間，海拔50～2158 m。該地區屬于中亞熱帶季風氣候，年均氣溫19～22 ℃，無霜期300 d以上，年均日照時數1670～1980 h，年均降水量1400～1800 mm，相對濕度為79%。

1.1.2? 數據獲取? 在福建省南平市木荷人工林中設置213塊固定樣地（含標準地），根據福建省森林資源清查中相關技術標準和實際林業生產情況，考慮不同海拔、地貌類型、林分年齡、胸徑、樹高、立地質量、郁閉度等因子，采取隨機抽樣的方式，選取189株樣木，其中幼齡林19株、中齡林67株、近熟林49株、成熟林45株和過熟林9株;這些樣木所在樣地的立地質量等級分別為：（Ⅰ）肥沃15株，（Ⅱ）較肥沃68株，（Ⅲ）中等肥沃88株，（Ⅳ）貧瘠18株。根據國家氣象科學數據中心提供的數據顯示（http：//www.nmic. cn/），樣木所在區域1981—2010年的極端氣溫、年均降水量、年均氣溫及年均濕度基本趨于一致，具體解析木所在標準地基本信息見表1。對樣木伐倒后用精度為0.01 m的皮尺測量其全樹高，打枝后以2 m為區分段，在樹干的0.3、1.0、1.3、3.0、5.0 m等位置處用精度為0.01 cm的鋼尺分別測量出帶皮直徑（DOB）和去皮直徑（DIB）。

使用R語言簡單隨機抽樣功能，從這189株樣木中隨機抽取兩組數據，其中用于模型構建133株，其余部分用于模型檢驗（56株），進行本次研究模型的變量有：胸徑（diameter at breast height， DBH）、樹高（tree height， H）、林分年齡（stand age， T）、胸高處樹皮厚度（bark thickness at breast height， BBT）、相對樹皮厚度（relative bark thickness， RBT）、任意樹高處樹皮率（該處去皮直徑與帶皮直徑之比，DIB/DOB）、郁閉度（canopy closure， CC），樣木基本信息見表2。

1.2? 方法

1.2.1? 基礎模型確立? 本研究通過收集他人建立的14個模型，加上自建5個模型，共19個模型作為研究木荷樹皮厚度基礎模型（表3），探討木荷最優樹皮厚度模型。其中，樹皮模型分類方法是根據唐誠等[10]的研究結果。

1.2.2? 模型檢驗和評價? 使用最小二乘法對以上19個基礎模型進行擬合，計算各參數值，并進行顯著性檢驗。經過篩選，對于滿足所有參數均具有顯著性的基礎模型，選取數學統計量指標：決定系數（determination coefficient， R2）、均方根誤差（root mean square error， Erms）、和方差（the sum of squares due to error， Ess）、赤池信息準則（akaike information criterion， AIC）、貝葉斯準則（bayesian information criterion， BIC）等5個統計指標來評估基礎模型擬合效果。判別標準為：R2越大，精度越高;Erms、Ess、AIC、BIC越小，精度越高;模型擬合效果也就越好。具體計算公式如下：

式中：yi為第i個實測值，i為第i個預測值，i為實測值的平均值，n為用于擬合模型的實測值數目，q為模型參數數目，為對數似然函數的最大值。

1.2.3? 模型優選方法? TOPSIS法又可稱為優劣距離法，1981年由Hwang和Yoon等[24]提出，它是一種有效的對多目標決策方法。但TOPSIS法在按歐氏距離進行優劣排序時，常出現接近于理想解，同時也可能出現貼近于負理想解的結果[25-26]。因此，需要用含熵值的TOPSIS法來彌補不足。熵值法是用來判斷評價指標離散程度，指標離散程度越大則對綜合評價影響也就越大[27-28]。具體計算步驟如下：

（1）根據模型求解結果構建m個評價對象、n個模型評價指標原始矩陣R。

（2）將原始矩陣進行歸一化處理，得到無量綱化矩陣B，指標范圍在0～1之間，值越大越優，最優值等于1，最差值為0。

xqmin為指標最小值，xqmax為指標最大值，rpq為指標歸一化值，fpq為一致性指標，Hq為指標信息熵，ωq即為第q個指標的熵值，其中ωq不能大于1，lnn必大于0。

1.2.4? 啞變量設置? 啞變量，又稱指示變量或虛擬變量，是對分類變量或定性因子進行處理的一種常用方法[29]。啞變量一般將定性因子轉為定量因子，通常取值為0或1。在林業相關回歸分析和數量化方法上經常引入啞變量，這不僅可以減少工作量，還可以提高計算準確度和模型精度[30-32]。本研究將通過優選基礎模型，在不同類型樹皮厚度模型中引入齡組和立地質量2個啞變量，用定性代碼來表示，以構建具有相容性的模型[33]。

以齡組（i）做為啞變量為例，不同齡組使用定量變量Ji，表示如下：

當樣木齡組為幼齡林時，取J1=1，J2=0，J3=0，J4=0，J5=0;當樣木齡組為中齡林時，取J2=1，J1=0，J3=0，J4=0，J5=0;依此類推。同理，使用Ks表示立地質量，s=1，2，3，4，當樣木來自立地質量為Ⅰ肥沃時，取K1=1，K2=0，K3=0，K4=0;當樣木來自立地質量為Ⅱ較肥沃時，取K2=1，K1=0，K3=0，K4=0;依此類推。將Ji和Ks加入到各個模型參數中，以達到預測不同齡組和立地質量所包含的林分因子的目的。

1.3? 數據處理

采用Excel 2010、SPSS 21.0、R軟件（ver 3.6.1）進行各項數據處理分析。

2? 結果與分析

2.1? 基礎模型擬合參數檢驗

將189株樣木進行解析后進行模型擬合、參數求解并進行顯著性檢驗后發現（表4）：在胸高處樹皮厚度模型M2、M4、M5、M6中，均有1～3個參數預估值在5%水平與零差異不顯著，且4個模型中都含有T、H等自變量，而最優模型M1中只含有胸徑（DBH）自變量，說明胸高處樹皮厚度與胸徑更明顯相關;在任意高度處樹皮厚度模型M8、M9中，各有1個參數預估值在5%水平與零差異不顯著，分別對應自變量為RH、CC、T;在相對樹皮厚度模型M14、M15中，有1～2個參數預估值在5%水平與零差異不顯著，分別對應自變量為RH、CC;而在去皮直徑模型中未發現有參數不顯著因子。因此這8個模型不再進一步參與分析，將剩下參數均顯著的11個模型做進一步評價。

2.2? 基礎模型評價指標及優選

基于樣木數據，對各基礎模型進行擬合，結果顯示：胸高處樹皮厚度模型M1決定系數（R2 = 0.6981）在同類模型中最高;任意高度處樹皮厚度模型M10決定系數（R2 = 0.5540）大于模型M7、M8;相對樹皮厚度模型M12決定系數（R2 = 0.5056）大于模型M13;同時模型M1、M10、M12的Erms、Ess、AIC、BIC均為各類模型中最小（表5）。因此，模型M1、M10、M12分別為胸高處、任意高度處、相對樹皮厚度最優基礎模型。

由于去皮直徑基礎模型（M16、M17、M18、M19）的各評價指標互有優劣，因此，在對模型優選時，先對各評價指標進行無量綱化，再求解出各無量綱化值;然后利用無量綱化值計算各指標的熵值，其中Ess熵值最高為0.213102，R2熵值最低為0.161043;最后綜合計算得出模型M16最優解總距離是0.000306，距離最小，為去皮直徑最優基礎模型（表6）。

2.3? 含啞變量模型擬合參數檢驗和評價指標

為進一步研究模型，將齡組和立地質量啞變量分別同時加在基礎模型的不同參數上，運用R軟件進行分析與比較。結果顯示，含齡組和立地質量啞變量模型的精度均要優于同時含有齡組和立地質量的啞變量模型。因此選出所有參數預估值與零差異顯著模型進行評價。

由表7可知，在胸高處樹皮厚度啞變量模型M20、M21的參數a中，分別加入齡組和立地質量啞變量后，R2從0.6981提高到0.9769和0.9758，但在Erms、Ess、AIC、BIC等指標上有略微增大，需進一步評價;在任意高度處樹皮厚度啞變量模型M22、M23的參數a中加入齡組和立地質量啞變量后，其中模型M22的R2為0.9214，高于模型M23、M10，并且模型M22的Erms、Ess、AIC、BIC均小于同類其他模型，因此模型M22可直接進行模型t檢驗;在相對樹皮厚度啞變量模型M24的參數a、b、c中加入齡組啞變量，以及在模型M25和M26的參數a中分別加入齡組和立地質量啞變量后，其中模型M25的R2最高為0.9111，但在Erms、Ess、AIC、BIC上3個啞變量模型和基礎模型互有高低，需要進一步篩選;去皮直徑啞變量模型M27、M28的參數b中，分別加入齡組和立地質量啞變量后，R2均為0.9997，但在Erms、Ess、AIC、BIC等指標上無法直接選出最優，需進一步分析。

Wingerd D E， Wiant H V. Variables for predicting inside-bark diameters of upper stems of appalachian hardwoods[J]. Journal of Forestry-Washington， 1982， 80（12）： 791-792.

Maguire D A， Hann D W. Bark thickness and bark volume in southwestern oregon douglas-fir[J]. Western Journal of Applied Forestry， 1990， 5（1）： 5-8.

王曉明. 興安落葉松、白樺樹皮厚度與胸徑關系研究[J]. 林業科技情報， 2017， 49（4）： 45-47.

張如義， 鄧小東， 楊小建， 等. 樹高、胸徑、樹皮厚對厚樸產量的影響[J]. 經濟林研究， 2017， 35（3）： 239-243.

唐? 誠， 王春勝， 龐圣江， 等. 西南樺人工林樹皮厚度模型模擬[J]. 林業科學， 2017， 53（7）： 85-93.

嚴銘海， 方靜儀， 李睿宇， 等. 福州市濕地松人工林樹皮厚度模型研究[J]. 西南林業大學學報（自然科學）， 2018， 38（5）： 124-131.

王曉林， 蔡可旺， 姜立春. 落葉松樹皮厚度變化規律的研究[J]. 森林工程， 2011， 27（2）： 8-11.

郭孝玉， 孫玉軍， 馬? 煒， 等. 適于FVS的長白落葉松樹皮因子[J]. 東北林業大學學報， 2011， 39（10）： 28-31.

賈煒瑋， 孫赫明， 李鳳日. 包含啞變量的黑龍江省落葉松人工林碳儲量預測模型系統[J]. 應用生態學報， 2019， 30（3）： 814-822.

朱光玉， 胡? 松， 符利勇. 基于啞變量的湖南櫟類天然林林分斷面積生長模型[J]. 南京林業大學學報（自然科學版）， 2018， 42（2）： 155-162.

陳? 林， 董安強， 王發國， 等. 廣東南嶺國家級自然保護區疏齒木荷+福建柏群落結構與物種多樣性研究[J]. 熱帶亞熱帶植物學報， 2010， 18（1）： 59-67.

劉發林. 南方集體林區木荷次生林生長規律及經營技術研究[D]. 長沙： 中南林業科技大學， 2013.

阮傳成， 李振問， 陳誠和， 等. 木荷生物工程防火機理及應用[M]. 成都： 電子科技大學出版社， 1995.

林? 艷， 何紫迪， 毛積鵬， 等. 木荷基因組SSR位點開發及初步分析[J]. 熱帶作物學報， 2018， 39（9）： 1766-1771.

王秀花， 馬麗珍， 馬雪紅， 等. 木荷人工林生長和木材基本密度[J]. 林業科學， 2011， 47（7）： 138-144.

黃? 宇. 木荷無性系種子園生長與結實性狀的初步研究[J]. 熱帶作物學報， 2017， 38（2）： 213-217.

鄧志勇. 木荷提取物對小菜蛾和菜粉蝶生長發育的抑制效果[J]. 江蘇農業科學， 2013， 41（4）： 129-130.

黃光強. 用木荷樹皮粉末作為非殺傷性驅散武器填充劑的設想[J]. 森林公安， 2002（2）： 20.

Hwang C L， Yoon K. Multiple attribute decision making and applications[M]. New York： Springer-Verlag， 1981.

徐小萬， 雷建軍， 李? 穎， 等. 現蕾期辣椒耐高溫多濕性CA-TOPSIS綜合評定[J]. 熱帶作物學報， 2013， 34（9）： 1747-1751.

Suthar H A， Gadit J J. Multiobjective optimization of 2dof controller using evolutionary and swarm intelligence enhanced with TOPSIS[J]. Heliyon， 2019， 5（4）： 321-332.

張一達， 劉學錄， 范亞紅， 等. 基于改進TOPSIS法的蘭州市土地利用多功能性評價[J]. 干旱區地理， 2019， 42（2）： 444-451.

張浩渺. 基于熵權改進的TOPSIS模型在蘇州河水質綜合評價中的應用[J]. 西北水電， 2019（3）： 12-15.

曾偉生， 唐守正， 夏忠勝， 等. 利用線性混合模型和啞變量模型方法建立貴州省通用性生物量方程[J]. 林業科學研究， 2011， 24（3）： 285-291.

金? 龍， 吳志祥， 楊? 川， 等. 不同林齡橡膠凋落物葉分解特性與有機碳動態研究[J]. 熱帶作物學報， 2015， 36（4）： 698-705.

高東啟， 鄧華鋒， 王海賓， 等. 基于啞變量的蒙古櫟林分生長模型[J]. 東北林業大學學報， 2014， 42（1）： 61-64.

鄭冬梅， 曾偉生. 用啞變量方法構建東北落葉松和櫟類分段地上生物量模型[J]. 北京林業大學學報， 2013， 35（6）： 23-27.

唐守正， 郎奎建， 李海奎. 統計和生物數學模型計算： ForStat教程[M]. 北京： 科學出版社， 2009.

陳東來， 秦淑英. 樹皮厚度、樹皮材積與直徑和樹高相關關系的研究[J]. 河北林學院學報， 1994（3）： 248-250.

于世川， 張文輝， 李? 罡， 等. 黃龍山林區不同郁閉度對遼東櫟種群結構的影響[J]. 生態學報， 2017， 37（5）： 1537-1548.

Li R， Weiskittel A R . Estimating and predicting bark thickness for seven conifer species in the Acadian Region of North America using a mixed-effects modeling approach： comparison of model forms and subsampling strategies[J]. European Journal of Forest Research， 2011， 130（2）： 219-233.

池上評， 陳金章， 江傳陽， 等. 基于間隔期的福建柏人工林動態生長模型及應用[J]. 福建林學院學報， 2014， 34（4）： 304-308.