基于立地類型啞變量的馬尾松優(yōu)勢(shì)木枝下高模型研究

童建明

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

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

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

1 研究區(qū)概況

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

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

2 研究方法

2.1 數(shù)據(jù)來(lái)源

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

2.2 優(yōu)勢(shì)木選取

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

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

表1 建模數(shù)據(jù)基本統(tǒng)計(jì)量項(xiàng)目枝下高(HCB)/m樹(shù)高(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種常見(jiàn)的枝下高模型[11-14](表2)，作為研究枝下高的基礎(chǔ)模型，討論立地類型啞變量馬尾松優(yōu)勢(shì)木枝下高模型的最佳參數(shù)形式。

表2 枝下高基礎(chǔ)模型模型 表達(dá)式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ù)檢驗(yàn)

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

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

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

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

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

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

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

2.5.1 初始立地類型啞變量

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

2.5.2 聚類立地類型啞變量

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

3 結(jié)果與分析

3.1 基礎(chǔ)模型擬合與評(píng)價(jià)

基于樣地的平均優(yōu)勢(shì)木數(shù)據(jù)，對(duì)基本模型的枝下高模型進(jìn)行擬合，結(jié)果詳見(jiàn)表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è)基本模型都作為構(gòu)建初始立地類型啞變量枝下高的模型，其后再根據(jù)評(píng)價(jià)指標(biāo)選出最優(yōu)模型。

3.2 啞變量模型擬合與評(píng)價(jià)

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

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

根據(jù)數(shù)量化方法I結(jié)果得出6個(gè)立地因子對(duì)枝下高均顯著，不同因子水平組合為21個(gè)初始立地類型。利用初始立地類型作為啞變量加在4種基礎(chǔ)模型不同參數(shù)上擬合，所有模型擬合指標(biāo)結(jié)果詳見(jiàn)表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í)，其確定系數(shù)差值很小。

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

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

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

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

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

基于聚類后立地類型啞變量，獲得在模型M1的兩個(gè)不同參數(shù)下模擬的RMSE，R2，MSE，AIC和BIC，詳見(jiàn)表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ù)估計(jì)

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

表7 基礎(chǔ)模型參數(shù)值參數(shù) 估計(jì)值漸近標(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ù)均不錯(cuò)。立地類型啞變量模型中啞變量的估計(jì)值有很大不同，即此項(xiàng)探索中，聚類后的立地類型啞變量模型較好地體現(xiàn)了不同立地類型之間的不同，將模型形式進(jìn)行了簡(jiǎn)化。另外，啞變量模型中的形式參數(shù)都具有較好的穩(wěn)定性。

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

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

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

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

3.3 模型檢驗(yàn)

利用建模樣本對(duì)構(gòu)建的枝下高模型進(jìn)行檢驗(yàn)。由于調(diào)查的樣地?cái)?shù)量有限，在全部樣本數(shù)用于建模的情況下，本文模型檢驗(yàn)采取分徑階檢驗(yàn)的方式，即將樣本數(shù)以胸徑分組，最終分別以徑階12、14、16、18、20和22～24六組數(shù)據(jù)檢驗(yàn)?zāi)Ｐ停瑱z驗(yàn)指標(biāo)選擇確定系數(shù)(R2)、均方根誤差(RMSE)和平均系統(tǒng)誤差(MSE)3個(gè)評(píng)價(jià)指標(biāo)進(jìn)行評(píng)價(jià)和比較。結(jié)果見(jiàn)表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中可知，建模樣本檢驗(yàn)中，各徑階擬合結(jié)果均較好，但差異較大。在各徑階上，聚類后立地類型啞變量模型確定系數(shù)均優(yōu)于基礎(chǔ)模型，均方根誤差(RMSE)和平均系統(tǒng)誤差(MSE)也均表明良好。檢驗(yàn)結(jié)果較好地顯示了模型M1所建的馬尾松枝下高模型具有良好的全面切合性能。

4 結(jié)論與討論

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