童建明

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

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

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

1 研究區(qū)概況

湖南省位于中國中南部，24°38′—30°08′N，108°47′—114°15′E，在長江中游，省會長沙。陸地總面積大概是3 174.35萬km2，其中：山地占51%、盆地占7%、平原占13%、丘陵占29%。東部、西部和南部的地形很高，形成了一個東北有開口的馬蹄形形狀。屬亞熱帶季風(fēng)潮濕的大陸氣候，年均氣溫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)等變量指標(biāo)。

2.2 優(yōu)勢木選取

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

在優(yōu)勢樹種(組)為馬尾松的樣地中，在優(yōu)勢樹種(組)中選取3棵優(yōu)勢樹(包括次優(yōu)勢樹)，分別統(tǒng)計樹高和胸徑，取其平均值，得到樣地平均優(yōu)勢木[10]樹高和胸徑，根據(jù)這兩個數(shù)值在樣地內(nèi)選取與之最接近的一株立木，該立木數(shù)據(jù)即為構(gòu)建模型的原始數(shù)據(jù)。所有樣地構(gòu)建模型數(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標(biāo)準(zhǔn)差 3.23.23.50.14變異系數(shù)0.50.20.20.15

2.3 模型選擇

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

表2 枝下高基礎(chǔ)模型模型 表達式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個標(biāo)準(zhǔn)進行評估和比較。其中R2反映模型的適用性，R2越接近1，模型就越準(zhǔn)確；RMSE值越小，模型就越準(zhǔn)確；MSE是反映一定范圍(3%或5%)內(nèi)擬合效果的關(guān)鍵指標(biāo)，接近0時，可以得到最佳效果。具體的公式為：

2.5 啞變量模型構(gòu)建

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

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

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

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

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

2.5.1 初始立地類型啞變量

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

2.5.2 聚類立地類型啞變量

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

3 結(jié)果與分析

3.1 基礎(chǔ)模型擬合與評價

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

表3 基礎(chǔ)模型擬合結(jié)果模型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不收斂——

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

3.2 啞變量模型擬合與評價

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

3.2.1 初始立地類型啞變量擬合結(jié)果

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

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

3.2.2 聚類立地類型啞變量擬合結(jié)果

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

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

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

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

表6 模型不同參數(shù)擬合結(jié)果模型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

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

3.2.3 模型參數(shù)估計

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

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

在模型M1中，基礎(chǔ)模型和啞變量模型的全部參數(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ù)形式。用枝下高實測值減去枝下高預(yù)測值得到殘差，殘差圖以枝下高預(yù)測值為橫軸，殘差為縱軸進行描繪，具體結(jié)果見模型M1基礎(chǔ)的殘差圖(圖1)和聚類后立地類型啞變量最優(yōu)參數(shù)形式的殘差圖(圖2)。

圖1 基礎(chǔ)模型殘差圖

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

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

3.3 模型檢驗

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

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

4 結(jié)論與討論

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