童建明

(國家林業(yè)和草原局中南調(diào)查規(guī)劃設計院，長沙 410014)

枝下高(Heighttocrownbase)一般指的是直立樹干上第一根活樹枝到地面的高度，它是單木樹冠特征的一個重要指標[1]，不僅能反映樹木的生長活力和立地生產(chǎn)力，還是林分內(nèi)競爭水平等的外在直觀體現(xiàn)[2]。大量研究表明，枝下高常被作為輸入變量應用于冠形模型、冠幅模型以及生物量模型等，此外也可用于計算單株或林分生長收獲模型的關鍵預測變量，如冠長、冠長率和其它等。但是，在實際生產(chǎn)經(jīng)營活動中對立木枝下高的獲取較為困難，尤其是在林況通透性差、郁閉度高的林分中，其測量誤差較大，造成了大量人力和物質(zhì)資源的浪費。同時，因為工作人員的自我意識，在測量時對第一活枝的判斷會存在偏差，從而影響數(shù)據(jù)的準確性；另外，同一樹種由于立地類型的不同，其枝下高也會存在較大的差異。因此在枝下高的廣義模型中，僅僅以樹高和胸徑作為自變量建立枝下高的單一模型，不足以反映所有可能的枝下高與相關測樹因子的關系，而基于立地類型啞變量為枝下高模型的研究，嘗試提供了一種有效的思路與途徑。

啞變量(dummyvariable)，也叫做指示變量、假變量。在啞變量基礎上的回歸分析方法，可參照李希非等[3]的探索。目前大部分的探索成果表明，啞變量是處理分類變量和定性因子經(jīng)常使用的一種方法，統(tǒng)計學中的每一種量化方法全部都有涉及啞變量處理的問題[4-6]。同時啞變量模型廣泛運用于不同的建模實踐與回歸分析中[7]，該模型能有效表示調(diào)查中的定性因子，模型內(nèi)既能體現(xiàn)復雜的林分結構，也能反映海拔、坡度等立地因子的不同水平。

1 研究區(qū)概況

湖南省位于中國中南部，24°38′—30°08′N，108°47′—114°15′E，在長江中游，省會長沙。陸地總面積大概是3 174.35萬km2，其中：山地占51%、盆地占7%、平原占13%、丘陵占29%。東部、西部和南部的地形很高，形成了一個東北有開口的馬蹄形形狀。屬亞熱帶季風潮濕的大陸氣候，年均氣溫16～18 ℃度至攝氏度，年日照時間1 300～1 800 h；年均降水量為1 200～1 700 mm，降雨量豐富，是我們國家降雨量較多的一個省份[8]；土壤重點是黃色土壤與紅色土壤，沖積土壤與紫色土壤其次。

馬尾松(Pinusmassoniana)，是松科樹種，廣泛分布于我國中南部，北自河南、魯南，南自廣東、廣西，東自湖南、臺灣，西自四川中部、貴州[9]，是中國南部主要的木材品種，具有很高的經(jīng)濟價值。

2 研究方法

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

數(shù)據(jù)采集于2020年9—11月，在湖南省懷化、邵陽和益陽三地共調(diào)查了42塊馬尾松人工林臨時樣地，測量樣地內(nèi)胸徑大于5 cm的每株立木。樣地調(diào)查內(nèi)容主要涉及樹高、枝下高、胸徑和其它測樹因子，以及坡面、海拔、坡度、土壤類型、土壤厚度和其它立地因素。優(yōu)勢樹種(組)由實測數(shù)據(jù)計算，選取優(yōu)勢木高(HD)等變量指標。

2.2 優(yōu)勢木選取

以樣地內(nèi)所有馬尾松胸高斷面積(ΣBAi)與所有樹種胸高斷面積(BA)之和為商，計算出各樣地馬尾松的組成系數(shù)(XSi)，也就是XSi=BAi/BA。為劃分優(yōu)勢樹種(組)，若ΣXSi≥0.65，則認為該臨時樣地優(yōu)勢樹種(組)為馬尾松，樣地保留，否則舍棄。

在優(yōu)勢樹種(組)為馬尾松的樣地中，在優(yōu)勢樹種(組)中選取3棵優(yōu)勢樹(包括次優(yōu)勢樹)，分別統(tǒng)計樹高和胸徑，取其平均值，得到樣地平均優(yōu)勢木[10]樹高和胸徑，根據(jù)這兩個數(shù)值在樣地內(nèi)選取與之最接近的一株立木，該立木數(shù)據(jù)即為構建模型的原始數(shù)據(jù)。所有樣地構建模型數(shù)據(jù)，詳見表1。

表1 建模數(shù)據(jù)基本統(tǒng)計量項目枝下高(HCB)/m樹高(HD)/m胸徑(D)/cm高徑比(RHD)平均值 6.915.617.10.92最小值 18.5110.69最大值 13.52224.41.29標準差 3.23.23.50.14變異系數(shù)0.50.20.20.15

2.3 模型選擇

本研究中，選擇5種常見的枝下高模型[11-14](表2)，作為研究枝下高的基礎模型，討論立地類型啞變量馬尾松優(yōu)勢木枝下高模型的最佳參數(shù)形式。

表2 枝下高基礎模型模型 表達式M1HCB= HD(1+exp(X))M2HCB= HD(1+exp(X))M3HCB= HD6(1+exp(X))M4HCB= HD×(1-exp(X))M5HCB= HD×(a+exp(X)) 注: 式中,X=b0+b1?D;其中 a ,b0,b1為模型參數(shù)。

2.4 模型及參數(shù)檢驗

模型用確定系數(shù)(R2)、均方根誤差(RMSE)和平均系統(tǒng)誤差(MSE)3個標準進行評估和比較。其中R2反映模型的適用性，R2越接近1，模型就越準確；RMSE值越小，模型就越準確；MSE是反映一定范圍(3%或5%)內(nèi)擬合效果的關鍵指標，接近0時，可以得到最佳效果。具體的公式為：

2.5 啞變量模型構建

立地類型對林木枝下高影響較顯著，為探索這種顯著性，首先將初始立地類型作為啞變量應用在基礎模型中擬合，再采用R語言k均值算法(k-means)對初始立地類型分級處理，其分類數(shù)標準為聚類精度≥0.99[15]，基于聚類結果最終構建含聚類后的立地類型啞變量模型。

本文在構建馬尾松優(yōu)勢木枝下高模型時，在確保預測準確性的基礎上，盡可能簡化了模型，建立了具有立地類型啞變量的優(yōu)勢木枝下高模型，以確定模型的最佳模式。

那么基于模型M1的啞變量模型的形式可表示為：

式中:Zi為啞變量，bi為對應的具體參數(shù)或局部參數(shù)。

按照立地類型的分類、聚類，立地類型涉及6個啞變量；即Ⅰ，Ⅱ，Ⅲ，Ⅳ，Ⅴ和Ⅵ。如果立地類型是Z1的時候，取Z1=1；那么Z2，Z3，…，Z6均為0，以此類推。

2.5.1 初始立地類型啞變量

對于初始立地類型的劃分，通過數(shù)量化方法I得出6個立地因子(海拔、坡度、坡向、坡位、土層厚度、土壤類型)對枝下高的顯著影響，并根據(jù)《湖南省森林資源規(guī)劃設計調(diào)查技術規(guī)程(2013年版)》對所有立地因子進行組合，每個立地因子水平組合為一個立地類型[16-17]，42塊樣地共劃分為21個初始立地類型。

2.5.2 聚類立地類型啞變量

基于最優(yōu)基礎模型考慮21個初始立地類型啞變量擬合的得分值，為了研究方便以及考慮模型精度，采用R語言k均值算法(k-means)對各因子分級處理，允許合并后因子水平信息損失≤1%。

3 結果與分析

3.1 基礎模型擬合與評價

基于樣地的平均優(yōu)勢木數(shù)據(jù)，對基本模型的枝下高模型進行擬合，結果詳見表3。

表3 基礎模型擬合結果模型R2RMSEMSE/%M10.320 42.738 3-2.594 0M20.320 22.738 8-2.632 5M30.319 82.739 6-2.678 2M40.321 52.736 2-2.430 5M5不收斂——

結果顯示，模型M1—M4的確定系數(shù)都較低。其中，模型M4的確定系數(shù)最大，為0.321 5；模型M3的確定系數(shù)最小，為0.319 8；模型M5擬合結果不收斂?？紤]到4種基礎模型(M1—M4)擬合的確定系數(shù)差值均小于0.001 7，并且該模型具有固有的可變性。所以，四個基本模型都作為構建初始立地類型啞變量枝下高的模型，其后再根據(jù)評價指標選出最優(yōu)模型。

3.2 啞變量模型擬合與評價

利用Forstat軟件中非線性混合效應模塊，以初始立地類型為啞變量加在4種候選枝下高模型不同參數(shù)上擬合；依據(jù)R2，AIC和BIC進行評價。再分析啞變量在不同模型上的擬合效果，選出最優(yōu)模型來構建含聚類后立地類型啞變量的候選模型。

3.2.1 初始立地類型啞變量擬合結果

根據(jù)數(shù)量化方法I結果得出6個立地因子對枝下高均顯著，不同因子水平組合為21個初始立地類型。利用初始立地類型作為啞變量加在4種基礎模型不同參數(shù)上擬合，所有模型擬合指標結果詳見表4。

表4 初始立地類型啞變量擬合結果模型R2AICBICM1-b00.773 8198.135 4221.114 9M1-b10.772 4313.831 5 336.733 4 M2-b00.773 8183.764 2210.505 6M2-b10.770 2301.255 2324.234 7M3-b00.773 8147.386 9170.366 4M3-b10.767 7265.059 7288.039 2M4-b00.773 8237.653 1260.632 6M4-b10.771 8354.972 9377.952 4

由表4可知，加入初始立地類型啞變量后，模型確定系數(shù)從0.319 8～0.321 5提高到0.767 7～0.773 8，且收斂度均較好，4種模型的確定系數(shù)差值均小于0.006 1，都是啞變量加在模型參數(shù)b0上最大，均為0.773 8。另外啞變量加在同一模型的不同參數(shù)上時，其確定系數(shù)差值很小。

基于4種模型的擬合結果，綜合考慮評價指標，最終選擇模型M1作為后期研究的候選模型。同時，由于初始立地類型啞變量加在模型M1上的兩個不同參數(shù)(b0和b1)得到擬合的確定系數(shù)很接近，分別為0.773 8和0.772 4，考慮到模型內(nèi)在的相容性，將分別以模型M1兩個不同參數(shù)的擬合結果來構建含聚類后立地類型啞變量模型，最終選擇最優(yōu)模型形式。

3.2.2 聚類立地類型啞變量擬合結果

根據(jù)42塊樣地初始劃分的21個立地類型，分別以模型M1兩個參數(shù)形式擬合的初始立地類型得分值聚類，聚類后分類數(shù)對應的立地編號及數(shù)量如表5。

表5 立地類型聚類模型M1-b0模型M1-b1立地類型編號數(shù)量立地類型編號數(shù)量Ⅰ16Ⅰ17Ⅱ4Ⅱ6Ⅲ12Ⅲ10Ⅳ5Ⅳ4Ⅴ1Ⅴ1Ⅵ4Ⅵ4

每個樣本的初始立地類型被轉換成相應的類，并作為啞變量添加到模型M1的相應參數(shù)中進行擬合，也就是說模型M1的參數(shù)b0上的聚類立地類型啞變量僅被添加到參數(shù)b0中進行擬合，而參數(shù)b1上的聚類立地類型啞變量僅被添加到參數(shù)b1中進行擬合，分析并比較兩種擬合結果。

基于聚類后立地類型啞變量，獲得在模型M1的兩個不同參數(shù)下模擬的RMSE，R2，MSE，AIC和BIC，詳見表6。

表6 模型不同參數(shù)擬合結果模型R2RMSEMSE/ %AICBICM1-b00.769 21.595 8-1.031 3184.413 1196.855 8M1-b10.767 01.603 6-1.140 5211.969 1224.411 9

結果表明，在聚類后添加立地類型啞變量后，模型的擬合結果要比基本模型好得多，但相比初始立地類型啞變量擬合的確定系數(shù)有所降低。其中，啞變量加在模型M1參數(shù)b0上的R2最大，為0.769 2；RMSE和MSE最小，分別為1.595 8和-1.0313；AIC和BIC值也優(yōu)于參數(shù)b1上的擬合結果。因此，聚類立地類型啞變量中，模型M1參數(shù)b0上的形式擬合最優(yōu)。

3.2.3 模型參數(shù)估計

基礎模型參數(shù)值和立地類型啞變量模型參數(shù)估計，詳見表7和表8。

表7 基礎模型參數(shù)值參數(shù) 估計值漸近標準差參數(shù)下限參數(shù)上限b00.564 30.596 4-0.641 11.769 6b1-0.019 80.032 4-0.085 20.045 6

在模型M1中，基礎模型和啞變量模型的全部參數(shù)均不錯。立地類型啞變量模型中啞變量的估計值有很大不同，即此項探索中，聚類后的立地類型啞變量模型較好地體現(xiàn)了不同立地類型之間的不同，將模型形式進行了簡化。另外，啞變量模型中的形式參數(shù)都具有較好的穩(wěn)定性。

表8 立地類型啞變量模型參數(shù)估計參數(shù)估計值SDtPⅠb00.491 40.539 10.911 50.368 2Ⅱb0-0.852 70.519 1-1.642 80.109 4Ⅲb0-0.395 00.570 0-0.692 90.492 9Ⅳb0-0.164 90.457 8-0.360 30.720 8Ⅴb01.748 30.979 41.785 00.082 9Ⅵb01.439 60.659 22.183 90.035 8b10.007 60.028 20.269 70.789 0

3.2.4 最優(yōu)模型殘差圖

按照擬合的成果，將選擇一種最理想的模型M1最好的啞變量參數(shù)形式。用枝下高實測值減去枝下高預測值得到殘差，殘差圖以枝下高預測值為橫軸，殘差為縱軸進行描繪，具體結果見模型M1基礎的殘差圖(圖1)和聚類后立地類型啞變量最優(yōu)參數(shù)形式的殘差圖(圖2)。

圖1 基礎模型殘差圖

圖2 聚類立地類型啞變量殘差圖

圖1和圖2直觀地反映了啞變量模型方法的優(yōu)勢，以聚類立地類型啞變量進行估計后的枝下高殘差圖要優(yōu)于基礎模型的擬合；另外，啞變量模型的枝下高估計值與實測值相差程度都較小，也就是說，各對應的立地類型啞變量模型得到的估計值與基礎模型估計值差異較小，這驗證了模型之間的相容性。

3.3 模型檢驗

利用建模樣本對構建的枝下高模型進行檢驗。由于調(diào)查的樣地數(shù)量有限，在全部樣本數(shù)用于建模的情況下，本文模型檢驗采取分徑階檢驗的方式，即將樣本數(shù)以胸徑分組，最終分別以徑階12、14、16、18、20和22～24六組數(shù)據(jù)檢驗模型，檢驗指標選擇確定系數(shù)(R2)、均方根誤差(RMSE)和平均系統(tǒng)誤差(MSE)3個評價指標進行評價和比較。結果見表9。

表9 模型分徑階擬合結果徑階/cm樣本數(shù)R2RMSEMSE/%1250.999 00.060 2-0.221 5 1480.694 71.886 50.134 7 1690.960 20.509 1-0.124 9 1870.855 51.301 9-1.521 4 2080.751 71.046 7-0.318 0 22～2450.669 92.491 2-3.045 9

從表9中可知，建模樣本檢驗中，各徑階擬合結果均較好，但差異較大。在各徑階上，聚類后立地類型啞變量模型確定系數(shù)均優(yōu)于基礎模型，均方根誤差(RMSE)和平均系統(tǒng)誤差(MSE)也均表明良好。檢驗結果較好地顯示了模型M1所建的馬尾松枝下高模型具有良好的全面切合性能。

4 結論與討論

因為此項探索的數(shù)據(jù)收集僅限于湖南三個地方，所以模型的應用范圍具有局限性。在優(yōu)勢樹的選擇中，選擇優(yōu)勢樹種(組)中的三棵優(yōu)勢木，沒有與其他選取方法數(shù)據(jù)作對比，在以后的研究中可嘗試以兩種以上不同優(yōu)勢木來分析，如最高優(yōu)勢木(優(yōu)勢樹種組中最高的一株立木)，平均優(yōu)勢木等。除此之外，關于立地類型啞變量，立地類型啞變量模型法的應用受到限制，因為在此項探索中全部的樣地初始的立地類型并不包括湖南地區(qū)所有的立地類型。此篇論文運用R語言的k-means算法將初始立地類型聚類成若干類作為啞變量，簡化啞變量類型是為了說明此方法的應用，因為不包括這一區(qū)域的全部立地類型，其實用性需進一步探討。研究考慮了立地類型啞變量，為區(qū)域性馬尾松人工林枝下高模型的研究提供了一種有效的思路與途徑。