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

郭強 劉蔚漪 輝朝茂 官鳳英 鄒學明

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

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

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

1 研究區概況

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

2 研究方法

2.1 樣地選擇與數據調查

以研究區班洪鄉處連續分布的整個巨龍竹林分為試驗對象(見表1)。2019年12月份，篩選得到健康且非林緣處生長的巨龍竹共計93叢，對其竹高達到主林層及以上的2 090株巨龍竹進行直徑和年齡檢尺(見表2)。

表1 研究樣地基本概況

表2 研究樣地的巨龍竹直徑與年齡分布統計

2.2 分布模型選擇

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

一元3參數威布爾概率密度函數分布模型，即式(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為位置參數，與林木直徑分布最小值有關；b為尺度參數，與林木直徑分布范圍有關；c為形狀參數，與林木直徑分布形狀有關。

2.2.2 二元威布爾生存函數模型

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

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

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

(2)

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

2.2.3 改進的二元威布爾概率密度函數模型

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

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

g7+g8·T+g9·D；

將相關回歸式對應代入式(1)后，導出2種改進的二元威布爾概率密度函數分布模型，即式(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為改進的二元威布爾分布模型Ⅰ，即式(3)的參數；q1、…、q10為改進的二元威布爾分布模型Ⅱ，即式(4)的參數。

利用Office Excel 2007和SPSS 19.0進行數據整理和分析。

3 結果與分析

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

3.1.1 巨龍竹直徑分布服從威布爾分布假設檢驗

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

圖1 巨龍竹直徑分布概率散點圖

圖2 巨龍竹直徑分布無趨勢概率散點圖

同時，利用柯爾莫可洛夫-斯米洛夫檢驗(K-S)檢驗法對巨龍竹直徑分布進行假設檢驗。首先假設其直徑分布服從威布爾分布，計算其統計量Dn=sup[Fn(x)-F(x)]=0.010 959，經查表(0.05水平)得其臨界值D=0.029 749，通過檢驗可知Dn

3.1.2一元3參數威布爾概率密度函數模型求解與評價

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

表3 最后5次迭代過程

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。

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

圖3 巨龍竹直徑分布觀測數據與威布爾分布擬合

3.2 巨龍竹直徑與年齡的二元威布爾生存函數模型求解與評價

3.2.1 二元威布爾生存函數模型求解

3.2.2 二元威布爾生存函數模型評價

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

表4 巨龍竹生存干數的實際觀測值與理論分布值

3.3 巨龍竹直徑與年齡的二元威布爾概率密度函數模型求解與評價

3.3.1 改進二元威布爾概率密度函數模型求解

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

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

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

3.3.2 2種改進二元威布爾概率密度函數模型評價

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

表5 巨龍竹直徑與年齡的實測概率值和理論概率值

4 討論

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

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

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

4.2 林木直徑與年齡的二元威布爾分布模型適用性評價

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

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

5 結論

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

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