威布爾分布模型在巨龍竹林直徑與年齡分布特征估測(cè)中的應(yīng)用1)

郭強(qiáng) 劉蔚漪 輝朝茂 官鳳英 鄒學(xué)明

(國(guó)際竹藤中心，北京，100102) (西南林業(yè)大學(xué)) (國(guó)際竹藤中心) (滄源縣林業(yè)和草原局)

林分直徑結(jié)構(gòu)反映了林木個(gè)體隨直徑變化的分布狀況，與林木的生長(zhǎng)發(fā)育、生產(chǎn)力、生物量和林分的生態(tài)穩(wěn)定性等密切相關(guān)，是重要且關(guān)鍵的生長(zhǎng)特征因子和參數(shù)推導(dǎo)因子[1]。利用相關(guān)分布模型，如正態(tài)(Normal)分布、威布爾(Weibull)分布、邏輯斯蒂(Logistic)分布等，可以有效地表征林木直徑分布規(guī)律，其中威布爾分布模型因其靈活性大、適應(yīng)性強(qiáng)、參數(shù)易于求解等優(yōu)勢(shì)而被廣泛用于林木直徑分布擬合和預(yù)估[2-4]，如蒙古櫟(Quercusmongolica)、油松(Pinustabulaeformis)、杉木(Cunninghamialanceolata)等[5-7]。通常林木直徑威布爾分布模型，主要指一元3參數(shù)威布爾分布模型，其位置參數(shù)(a)、尺度參數(shù)(b)、形狀參數(shù)(c)與林木直徑分布特征密切相關(guān)(當(dāng)a=0時(shí)，即為2參數(shù)分布)。但綜合看，分析直徑與其它結(jié)構(gòu)因子的二元威布爾分布模型比一元威布爾分布模型能更全面地反映林木直徑分布的實(shí)際情況。為此，一些學(xué)者對(duì)林木直徑與其它因子的二元威布爾分布模型進(jìn)行了探索[8]，在一定程度上為二元威布爾分布模型的理論研究和林業(yè)應(yīng)用提供了參考[9]；但其總體上還不成熟，相關(guān)研究成果也較少，還需進(jìn)一步研究和完善。

巨龍竹(Dendrocalamussinicus)為禾本科竹亞科牡竹屬特大型合軸叢生竹種，是中國(guó)西南地區(qū)特有珍惜竹種，同時(shí)也是迄今發(fā)現(xiàn)徑級(jí)最大的竹種，有“世界竹王”之美譽(yù)[10]。但是長(zhǎng)期以來(lái)，由于巨龍竹形態(tài)巨大、分布偏遠(yuǎn)，且分布山區(qū)地形與氣候復(fù)雜等原因，對(duì)其相關(guān)研究較少，且較為局限，對(duì)其的保護(hù)和研究具有一定的緊迫性[11-12]。由于廣泛適用于一般林木樹(shù)種的威布爾分布模型，對(duì)這類特大型叢生竹種的適用性還不確定，這一珍稀竹種的直徑分布規(guī)律還不清晰，為此，本研究以滇西南連續(xù)分布的整個(gè)巨龍竹林分為試驗(yàn)對(duì)象，于2019年12月份，篩選得到健康且非林緣處生長(zhǎng)的巨龍竹共計(jì)93叢，對(duì)其竹高達(dá)到主林層及以上的2090株巨龍竹進(jìn)行直徑和年齡檢尺統(tǒng)計(jì)；采用一元3參數(shù)威布爾概率密度函數(shù)模型、二元威布爾生存函數(shù)模型、導(dǎo)出的2種改進(jìn)二元威布爾概率密度函數(shù)模型，分析滇西南巨龍竹林直徑與年齡分布特征。旨在為巨龍竹林業(yè)工作開(kāi)展和林用二元威布爾分布模型研建提供參考。

1 研究區(qū)概況

研究區(qū)位于云南省臨滄市滄源佤族自治縣(98°52′～99°43′E、23°4′～23°40′N)，屬亞熱帶季風(fēng)氣候類型，干濕季分明，一般夏秋季為雨季、春冬季為旱季；年平均氣溫21 ℃，年積溫約7 500 ℃，年降水量約2 100 mm，月均空氣相對(duì)濕度約80%；土壤類型主要為紅壤、赤紅壤、磚紅壤等；林下分布植被主要為小葉肉實(shí)樹(shù)(Sarcospermagriffithii)、光葉火筒樹(shù)(Leeaglabra)、豬肚木(Canthiumhorridum)、多羽鳳尾蕨(Pterisdecrescens)、破壞草(Ageratinaadenophora)等。

2 研究方法

2.1 樣地選擇與數(shù)據(jù)調(diào)查

以研究區(qū)班洪鄉(xiāng)處連續(xù)分布的整個(gè)巨龍竹林分為試驗(yàn)對(duì)象(見(jiàn)表1)。2019年12月份，篩選得到健康且非林緣處生長(zhǎng)的巨龍竹共計(jì)93叢，對(duì)其竹高達(dá)到主林層及以上的2 090株巨龍竹進(jìn)行直徑和年齡檢尺(見(jiàn)表2)。

表1 研究樣地基本概況

表2 研究樣地的巨龍竹直徑與年齡分布統(tǒng)計(jì)

2.2 分布模型選擇

2.2.1 一元3參數(shù)威布爾概率密度函數(shù)模型

一元3參數(shù)威布爾概率密度函數(shù)分布模型，即式(1)：

f(D)=[(D-a)/b]c·[c/(D-a)]·exp{-[(D-

a)/b]c}=(c/b)·[(D-a)/b]c-1·exp{-

[(D-a)/b]c}。

(1)

式中：a為位置參數(shù)，與林木直徑分布最小值有關(guān)；b為尺度參數(shù)，與林木直徑分布范圍有關(guān)；c為形狀參數(shù)，與林木直徑分布形狀有關(guān)。

2.2.2 二元威布爾生存函數(shù)模型

二元威布爾生存函數(shù)分布模型[9]，即式(2)：

N(D>d,T>t)=n·P(D>d,T>t)=n·exp{-[(d-

a1)/b1]c1+(t/b2)c2}。

(2)

式中：a1為位置參數(shù)，與林木直徑分布最小值有關(guān)；b1、b2為尺度參數(shù)，分別與林木直徑和年齡分布范圍有關(guān)；c1、c2為形狀參數(shù)，分別與林木直徑和年齡分布形狀有關(guān)；n為數(shù)量參數(shù)，與林木分布總數(shù)量有關(guān)；P(D>d,T>t)為林木的生存概率。

2.2.3 改進(jìn)的二元威布爾概率密度函數(shù)模型

研究表明[8，13-14]，一元3參數(shù)威布爾分布模型的3個(gè)參數(shù)(a、b、c)與部分林分結(jié)構(gòu)因子具有多元線性或指數(shù)等回歸關(guān)系。因此，假設(shè)其分布模型的3個(gè)參數(shù)(a、b、c)和林木直徑、年齡分別具多元線性及指數(shù)回歸關(guān)系，則有：

c1=g1+g2·T+g3·D、b1=g4+g5·T+g6·D、a1=

g7+g8·T+g9·D；

將相關(guān)回歸式對(duì)應(yīng)代入式(1)后，導(dǎo)出2種改進(jìn)的二元威布爾概率密度函數(shù)分布模型，即式(3)、式(4)：

f1(T,D)=g10·[(g1+g2·T+g3·D)/(g4+g5·T+

g6·D)]·{[D-(g7+g8·T+g9·D)]/(g4+g5·T+g6·D)}(g1+g2·T+g3·D-1)·exp{-[D-

(g7+g8·T+g9·D)]/(g4+g5·T+g6·

D)}(g1+g2·T+g3·D)；

(3)

(4)

式中：g1、…、g10為改進(jìn)的二元威布爾分布模型Ⅰ，即式(3)的參數(shù)；q1、…、q10為改進(jìn)的二元威布爾分布模型Ⅱ，即式(4)的參數(shù)。

利用Office Excel 2007和SPSS 19.0進(jìn)行數(shù)據(jù)整理和分析。

3 結(jié)果與分析

3.1 巨龍竹直徑的一元威布爾概率密度函數(shù)模型擬合及求解

3.1.1 巨龍竹直徑分布服從威布爾分布假設(shè)檢驗(yàn)

概率散點(diǎn)圖可檢驗(yàn)樣本數(shù)據(jù)是否服從某種理論分布，即當(dāng)符合既定理論分布時(shí)，各數(shù)據(jù)點(diǎn)在概率散點(diǎn)圖中近似呈一條直線地分布于右斜對(duì)角線位置，而在無(wú)趨勢(shì)概率散點(diǎn)圖中則呈離散分布。對(duì)現(xiàn)實(shí)巨龍竹林分直徑分布進(jìn)行威布爾分布檢驗(yàn)(見(jiàn)圖1、圖2)，由圖1、圖2可見(jiàn)，概率散點(diǎn)圖中的各數(shù)據(jù)點(diǎn)近乎一條直線地分布于“Y=X”線處，而無(wú)趨勢(shì)概率散點(diǎn)圖中的各數(shù)據(jù)點(diǎn)則呈離散狀分布于“Y=0”線的上下側(cè)。因此，可以認(rèn)為現(xiàn)實(shí)巨龍竹林分直徑分布服從威布爾分布。

圖1 巨龍竹直徑分布概率散點(diǎn)圖

圖2 巨龍竹直徑分布無(wú)趨勢(shì)概率散點(diǎn)圖

同時(shí)，利用柯?tīng)柲陕宸?斯米洛夫檢驗(yàn)(K-S)檢驗(yàn)法對(duì)巨龍竹直徑分布進(jìn)行假設(shè)檢驗(yàn)。首先假設(shè)其直徑分布服從威布爾分布，計(jì)算其統(tǒng)計(jì)量Dn=sup[Fn(x)-F(x)]=0.010 959，經(jīng)查表(0.05水平)得其臨界值D=0.029 749，通過(guò)檢驗(yàn)可知Dn

3.1.2一元3參數(shù)威布爾概率密度函數(shù)模型求解與評(píng)價(jià)

應(yīng)用牛頓迭代法求解一元3參數(shù)威布爾分布模型的3個(gè)參數(shù)值(a、b、c)。根據(jù)一元威布爾分布模型參數(shù)值與其分布形狀的關(guān)系，并結(jié)合巨龍竹直徑分布的散點(diǎn)圖，分別選取其初始值為5.3、12.0、2.5。經(jīng)過(guò)12次迭代后穩(wěn)定(表3為最后5次迭代結(jié)果)，并得最優(yōu)解：a=5.330 583、b=13.409 918、c=5.154 559。因此，可以確定，巨龍竹直徑分布的一元3參數(shù)威布爾概率密度函數(shù)最優(yōu)預(yù)估模型為：

表3 最后5次迭代過(guò)程

f(D)=[(D-5.330 583)/13.409 918]5.154 559·

[5.154 559/(D-5.330 583)]·exp[-(D-

5.330 583)/13.409 918]5.154 559，

R2=0.993 436。

根據(jù)巨龍竹直徑分布的實(shí)際觀測(cè)和分布模型擬合結(jié)果(見(jiàn)表3、圖3)分析可知，巨龍竹直徑一元3參數(shù)威布爾分布模型的擬合效果優(yōu)良，其擬合的決定系數(shù)R2=0.993 436、平均絕對(duì)誤差值ε=0.003 183。巨龍竹林分直徑分布呈單峰狀，主要分布于16～20 cm，其它徑階的巨龍竹數(shù)量則相對(duì)較少；其直徑分布的偏度值和峰度值分別為0.282 982、-1.378 235。

圖3 巨龍竹直徑分布觀測(cè)數(shù)據(jù)與威布爾分布擬合

3.2 巨龍竹直徑與年齡的二元威布爾生存函數(shù)模型求解與評(píng)價(jià)

3.2.1 二元威布爾生存函數(shù)模型求解

3.2.2 二元威布爾生存函數(shù)模型評(píng)價(jià)

巨龍竹直徑與年齡的二元威布爾生存函數(shù)模型可適用于其林分直徑與年齡分布規(guī)律研究，其擬合的決定系數(shù)R2=0.997 341、平均絕對(duì)誤差值ε=64.013 527。利用二元威布爾生存函數(shù)預(yù)估模型對(duì)巨龍竹的生存株數(shù)進(jìn)行預(yù)測(cè)，并與實(shí)際觀測(cè)值對(duì)比(見(jiàn)表4)。由表4可知，巨龍竹生存株數(shù)隨直徑分布的下降幅度在不同年齡特征的變化趨勢(shì)基本類似，即慢-快-慢的降幅變化，且包含的年齡結(jié)構(gòu)信息量越大，其降幅程度也越大，這間接地反映出巨龍竹林直徑分布呈單峰狀分布的特點(diǎn)。同時(shí)，將生存株數(shù)的分布模型預(yù)估值分別除以n(即2 145.405 447)，可相應(yīng)地得到巨龍竹林直徑與年齡分布的生存概率預(yù)測(cè)值。

表4 巨龍竹生存干數(shù)的實(shí)際觀測(cè)值與理論分布值

3.3 巨龍竹直徑與年齡的二元威布爾概率密度函數(shù)模型求解與評(píng)價(jià)

3.3.1 改進(jìn)二元威布爾概率密度函數(shù)模型求解

應(yīng)用牛頓迭代法求解2種改進(jìn)的巨龍竹直徑與年齡二元威布爾概率密度函數(shù)模型的相關(guān)參數(shù)，并將其參數(shù)估計(jì)值對(duì)應(yīng)代入其分布函數(shù)模型表達(dá)式式(3)、式(4)式后，則有：

圖4和圖5反映了2種改進(jìn)的巨龍竹直徑與年齡二元威布爾概率密度函數(shù)模型的擬合特征。由圖4、圖5可見(jiàn)，其分布形狀在直徑和分布概率的投影上，均呈一元威布爾分布狀，各年齡竹的直徑分布較為穩(wěn)定，但彼此仍略微差異，其總體上表現(xiàn)為新生竹比老竹的直徑有逐漸增大的趨勢(shì)，這與林分實(shí)際情況相吻合。可見(jiàn)，考慮竹年齡特征的二元威布爾分布模型，在反映巨龍竹林分直徑分布特征方面比一元威布爾分布模型更為客觀、全面。

圖5 巨龍竹直徑與年齡的二元威布爾分布模型Ⅱ

3.3.2 2種改進(jìn)二元威布爾概率密度函數(shù)模型評(píng)價(jià)

綜合2種改進(jìn)的二元威布爾概率密度函數(shù)模型擬合巨龍竹林直徑與年齡分布結(jié)果，并與現(xiàn)實(shí)林分實(shí)際觀測(cè)值對(duì)比(見(jiàn)表5)，由表5可見(jiàn)，2種改進(jìn)的二元威布爾分布模型對(duì)巨龍竹林直徑與年齡分布的擬合效果良好，其擬合的決定系數(shù)分別為0.919 323、0.906 260，平均絕對(duì)誤差值分別為0.000 072、0.001 325，可滿足基本的精度要求。根據(jù)現(xiàn)實(shí)林分巨龍竹總干數(shù)乘以表5中對(duì)應(yīng)徑階與齡階的理論分布概率值，可得到其林分直徑和年齡分布的估測(cè)預(yù)期值。

表5 巨龍竹直徑與年齡的實(shí)測(cè)概率值和理論概率值

4 討論

4.1 林木直徑與年齡的二元威布爾分布模型及求解

葛宏立等[9]在對(duì)浙江省毛竹(Phyllostachysheterocycla)直徑與年齡二元威布爾分布模型研究中，提出了一種二元威布爾生存函數(shù)模型，并認(rèn)為其研究方法在理論上可應(yīng)用于毛竹林分。本研究中，利用二元威布爾生存函數(shù)模型可較好地反映巨龍竹林分直徑與年齡分布特征，說(shuō)明這類分布模型，不僅可以應(yīng)用于林分，還可適用于巨龍竹這樣的特大型叢生竹種。但是，該類分布模型的擬合結(jié)果具有累積式分布特征，難以快速且方便地反映具體徑階與齡階的林木直徑與年齡分布狀況，而二元威布爾概率密度函數(shù)模型在這方面則具有一定的優(yōu)勢(shì)。本研究通過(guò)一元3參數(shù)威布爾分布模型參數(shù)(a、b、c)對(duì)林木直徑與年齡進(jìn)行多元線性、指數(shù)回歸，導(dǎo)出了2種改進(jìn)二元威布爾概率密度函數(shù)模型，并均具有良好的擬合效果，可滿足基本的精度要求，其中由多元線性回歸導(dǎo)出的改進(jìn)二元威布爾分布模型，與周國(guó)模等[8]研究結(jié)果相似。

本研究中，在利用迭代法求解2類3種二元威布爾分布模型的多個(gè)模型參數(shù)時(shí)，均存在參數(shù)初始值選取困難等難題，這常常使得模型參數(shù)值計(jì)算工作量大、擬合效果不佳。因此，如何巧妙且恰當(dāng)?shù)剡x取二元分布模型參數(shù)初始值，以便于高效且準(zhǔn)確地得到其模型參數(shù)最優(yōu)解，還需進(jìn)一步研究。

4.2 林木直徑與年齡的二元威布爾分布模型適用性評(píng)價(jià)

本研究中，二元威布爾生存函數(shù)模型在對(duì)巨龍竹直徑與年齡分布擬合時(shí)，盡管難以直觀地反映各個(gè)徑階與齡階的竹株具體分布特征，但其擬合精度較高，且模型參數(shù)值的生態(tài)學(xué)意義明確。本研究中，對(duì)于嘗試改進(jìn)的2種二元威布爾概率密度函數(shù)模型，盡管較為真實(shí)且方便地反映了巨龍竹林直徑與年齡分布情況，但其模型多個(gè)參數(shù)值的生態(tài)學(xué)意義還不明確，甚至其分布模型本身也有待于從數(shù)學(xué)理論層面做進(jìn)一步考證。同時(shí)，本研究未能考量巨龍竹“叢生”這一特性，如何在巨龍竹或其它叢生竹種的二元或多元威布爾分布模型研究中引入“竹叢”這一特征因素進(jìn)行分析，還需深入研究。

此外，本研究通過(guò)最小二乘的牛頓迭代法求解一元3參數(shù)威布爾概率密度函數(shù)模型擬合巨龍竹直徑分布概率、二元Weibull生存函數(shù)模型擬合巨龍竹直徑與年齡生存株數(shù)、2種改進(jìn)二元威布爾概率密度函數(shù)模型擬合巨龍竹直徑與年齡分布概率的參數(shù)最優(yōu)解，而針對(duì)不同的分布模型類型、參數(shù)估計(jì)方法、變量含義等會(huì)對(duì)其模型擬合效果產(chǎn)生怎樣的變化，仍需進(jìn)一步探討。

5 結(jié)論

利用一元3參數(shù)威布爾分布模型可以較好地描述巨龍竹林直徑分布特征；而利用巨龍竹直徑與年齡的二元威布爾生存函數(shù)模型和2種改進(jìn)二元威布爾概率密度函數(shù)模型，則可以更加真實(shí)地反映考慮年齡特征的巨龍竹林直徑分布規(guī)律。綜合各類威布爾分布模型擬合特征認(rèn)為，巨龍竹林直徑分布呈單峰狀，主要分布于16～20 cm，且具有新生竹比老竹直徑逐漸增大的變化趨勢(shì)。

本研究分別建立了巨龍竹林直徑一元3參數(shù)威布爾概率密度函數(shù)最優(yōu)預(yù)估模型、直徑與年齡的二元威布爾生存函數(shù)最優(yōu)預(yù)估模型、直徑與年齡的2種改進(jìn)二元威布爾概率密度函數(shù)最優(yōu)預(yù)估模型，并在此基礎(chǔ)上得到了巨龍竹林直徑與年齡的生存數(shù)量(概率)和分布概率的理論估測(cè)數(shù)表，可滿足基本的精度要求。相關(guān)研究方法和結(jié)果可為今后巨龍竹林業(yè)工作開(kāi)展和林用二元威布爾分布模型研建提供參考。