福建杉木樹冠外輪廓和樹冠體積相容性模型

吳丹子， 王成德， 李 倞， 劉 敏

（1. 北京林業(yè)大學(xué) 園林學(xué)院， 北京100083； 2. 北京林業(yè)大學(xué) 信息學(xué)院， 北京100083； 3. 北京林業(yè)大學(xué)城鄉(xiāng)生態(tài)環(huán)境北京實驗室， 北京100083）

樹冠外輪廓模型是指以樹木樹冠任意位置處半徑為因變量， 以冠幅、 冠長、 樹高或胸徑等因子為自變量的函數(shù)， 用來模擬樹冠形狀變化規(guī)律的數(shù)學(xué)表達式［1］。 樹冠外輪廓模型不僅能夠估計樹冠任意位置處的樹冠半徑和推導(dǎo)計算整個樹冠體積， 而且可以驅(qū)動林分樹木樹冠三維可視化模擬。 此外， 樹冠形態(tài)與樹冠體積對于間接估算樹木樹冠部分生物量具有重要意義［2］。 為了描述樹冠外輪廓形態(tài)變化規(guī)律和推導(dǎo)計算樹冠體積， 國內(nèi)外學(xué)者常利用以下方法預(yù)測樹冠相關(guān)變量， 包括簡單幾何形狀模擬、 經(jīng)驗?zāi)Ｐ头椒ā?樹冠輪廓模型積分法、 激光掃描方法等。 早期學(xué)者們利用規(guī)則幾何體描述樹冠外輪廓二維形狀， 常采用如下公式定義樹冠形狀：RCR=（1-RDH）k， 其中：RCR表示相對樹冠半徑，RCH為從樹冠基部到任意位置處樹冠長度與最大樹冠長度的比，k的取值決定了其幾何形狀， 當(dāng)k為0、 0.5、 1.0、 1.5 時分別表示圓柱形、 拋物線形、 圓錐形或凹面形等， 該模型被廣泛應(yīng)用于道格拉斯杉Pseudotsuga douglasii［3］、 西部鐵杉Tsuga heterophylla［4］、 火炬松Pinus taeda［5］、 華北落葉松Larix principis-rupprechtii［6］等樹冠形狀模擬和體積計算研究中。 PRETZSCH 等［7］、 SADONO［8］將樹冠分為 “陽冠” “陰冠” 兩部分， 分別采用拋物線體和圓錐體表示。 由于樹冠形狀并不是簡單規(guī)則形式， 不同樹種樹冠形狀并不相同， 樹冠體積大小與冠幅、 冠長、 胸徑、 樹高等林木因子之間存在相關(guān)關(guān)系， 研究者針對不同樹種利用經(jīng)驗?zāi)Ｐ头椒?gòu)建冪函數(shù)、 指數(shù)函數(shù)等多種形式樹冠體積模型， 如楠木Phoebe zhennan［9］、 樟子松Pinus sylvestrisvar.mongolica［10］、杉木Cunninghamia lanceolata［11］、 長白落葉松Larix olgensis［12］、 歐洲赤松Pinus sylvestris和云杉Picea asperata［13］、 黑荊樹Acacia mearnsii［14］， 也有研究者通過構(gòu)建枝條長度、 角度、 半徑等枝條屬性經(jīng)驗方程間接預(yù)測樹冠半徑和模擬樹冠形態(tài)變化規(guī)律， 如李鳳日［15］、 LI 等［16］、 GILMORE 等［17］、 姜立春等［18］， 但是存在樹冠枝條測量困難且誤差較大問題。 三維激光掃描等技術(shù)的發(fā)展， 不僅提供了快速數(shù)據(jù)獲取方式，而且能夠更加準(zhǔn)確獲取樹冠數(shù)據(jù)［19-21］， 但由于數(shù)據(jù)獲取成本高、 處理困難、 數(shù)據(jù)量大等限制了它的使用。 除了以上方法， 為了更加準(zhǔn)確預(yù)測樹冠體積， 一些學(xué)者利用多項式［22］、 冪函 數(shù)［23］、 修正Beta曲線［24］、 連續(xù)分段函數(shù)［25］等曲線形式描述樹冠外輪廓形態(tài)， 并將樹冠曲線圍繞樹干軸旋轉(zhuǎn)得到的旋轉(zhuǎn)體， 通過對旋轉(zhuǎn)體積分推導(dǎo)計算樹冠體積［26-28］。 該方法只需要簡單林木因子就能夠準(zhǔn)確模擬樹冠形態(tài)與樹冠體積。 目前， 樹冠外輪廓模型積分法計算體積研究， 一般將樹冠外輪廓與樹冠體積模型單獨進行研究， 分別擬合不同模型方程得到2 套模型參數(shù)； 由于外輪廓與體積方程之間存在內(nèi)在相關(guān)性， 采用傳統(tǒng)最小二乘方法分別求解模型參數(shù)， 不能保證2 個模型誤差同時最小， 無法滿足模型參數(shù)估計的漸進無偏性、 有效性和一致性。 通過對樹冠輪廓模型積分推導(dǎo)得到樹冠體積方程， 將2 個方程聯(lián)立構(gòu)建一致性方程組， 能夠有效解決以上模型參數(shù)估計存在的問題。 目前一致性方程組的研究主要集中在樹干削度與材積一致性方程組研究方面［29-30］， 關(guān)于利用樹冠外輪廓與樹冠體積一致性模型方面研究還未見報道。 本研究以福建地區(qū)杉木人工林為研究對象， 采用似乎不相關(guān)回歸方法， 構(gòu)建以最大樹冠半徑與相對冠長為自變量的樹冠外輪廓和樹冠體積一致性方程組模型， 研究構(gòu)建的模型能夠用于驅(qū)動樹木樹冠三維可視化模擬， 直觀地反映不同生長條件下樹冠生長活力及林木個體間樹冠重疊程度， 指導(dǎo)林分撫育間伐活動， 同時為通過樹冠體積間接推算樹冠生物量提供新思路。

1 研究區(qū)域與數(shù)據(jù)收集

福建省位于中國東南沿海地區(qū)（23°33′～28°20′N， 115°50′～120°40′E）， 地勢東南低西北高， 地形多丘陵少平原， 不同區(qū)域氣候相差較大， 其中東南沿海屬南亞熱帶氣候， 東北與西北區(qū)域?qū)僦衼啛釒夂颍?年平均氣溫為15.0～22.0 ℃， 年平均降水量為1 400～2 000 mm， 土壤類型主要包括紅土壤、 黃土壤、 山地草甸土壤， 主要用材樹種有杉木、 馬尾松Pinus massoniana、 巨尾桉Eucalyptus grandis×E.urophylla等。 杉木為杉科Taxodiaceae 喬木， 幼樹樹冠尖塔形， 大樹樹冠圓錐形， 為中國長江流域、 秦嶺以南地區(qū)栽培最廣、 生長快、 經(jīng)濟價值高的用材樹種。

研究數(shù)據(jù)來自于福建省順昌縣大歷和嵐下林場布設(shè)的杉木人工純林臨時樣地， 選擇不同齡組、 林分密度和立地條件類型設(shè)置標(biāo)準(zhǔn)樣地， 共設(shè)置98 塊30 m × 20 m 的樣地， 每個標(biāo)準(zhǔn)地內(nèi)選擇3～5 株樹木，共413 株杉木。 測量每株樹木胸徑（D）、 樹高（H）、 冠幅（CW）、 最大冠長（LCL）、 枝下高（HCB）， 樹冠長度（CHi）以及相應(yīng)樹冠半徑（CRi）， 其中：i（i=0.10、 0.25、 0.50、 0.75、 0.90）表示從樹冠基部到樹梢頂端的相對位置。 樹冠因子測量如圖1 所示， 該裝置由透明繪圖板、 三腳架以及照準(zhǔn)裝置組成， 具體使用方法如下： ①根據(jù)相似三角形原理， 在距離所測樹冠一定距離位置處（通常距離1 倍樹高）， 將該繪圖板固定在三腳架之上保持板面垂直于地面， 然后通過透明板觀察樹冠， 往后移動三腳架， 保持透明板與地面垂直， 直到通過透明板能夠觀察到整個樹冠。 ②用筆將樹冠外輪廓繪制在透明板之上， 首先繪制樹冠枝下高位置處的點A、 樹冠頂端的點B以及最大樹冠位置處的點C與點D， 然后從A點開始到B、C點，將樹冠輪廓繪制出來。 并將透明板上的樹冠輪廓草圖復(fù)制到有計算網(wǎng)格方格的透明硫酸紙上。 ③利用測高儀器和皮尺分別測量該樹冠的冠長L（LCL=A′B′）和冠幅（CW=E′F′）， 計算樹冠冠長測量值A(chǔ)′B′與繪圖紙上AB之間比值以及樹冠冠幅測量值E′F′與CD之間比值， 根據(jù)這2 個比值在繪圖紙上分別計算樹冠冠長0.10、 0.25、 0.50、 0.75、 0.90 位置處對應(yīng)的樹冠長度和半徑值。 為了盡可能減少樹冠測量誤差， 應(yīng)從多個方向觀察樹冠， 取測量均值。 通過以上調(diào)查方法收集到杉木樹冠調(diào)查數(shù)據(jù)情況如表1 所示。

圖1 杉木樹冠測量因子示意圖Figure 1 Schematic diagram of crown measurement factors of C. lanceolata

表1 樹冠調(diào)查數(shù)據(jù)的基本概況Table 1 Summary statistics of measurements of tree variables

2 研究方法

2.1 樹冠外輪廓模型與體積模型

本研究收集整理了國內(nèi)外文獻研究中常用于描述樹冠外輪廓形狀的模型方程， 具體形式如模型1、模型2、 模型3 以及模型4 所示。 其中： 模型1 被廣泛用于模擬多個樹種樹冠形態(tài)［3］， 該模型是否可用于杉木樹冠模擬需要進一步驗證， 其他3 個模型常用于模擬杉木樹冠［23］。 本研究將這4 個可積分模型方程作為聯(lián)立方程組中樹冠外輪廓模型的備選模型。 模型因變量為任意位置處樹冠半徑（CR）， 模型自變量為相對樹冠冠長（RCH）、 最大樹冠半徑（LCR）。

采用積分法計算樹冠體積， 從樹冠基部到樹木頂端積分樹冠外輪廓模型， 得到體積方程， 其中：VC為樹冠體積，LCL最大樹冠長度，CH為從樹冠基部到樹冠任意位置處樹冠長度， 推導(dǎo)得出樹冠體積模型方程。

式（5）～（8）中：RCH為相對冠長（RCH=CH/LCL， 樹冠基部為0， 樹冠頂部為1）；a0、a1、a2為模型系數(shù)。

2.2 似乎不相關(guān)非線性模型

假設(shè)定義如下非線性模型聯(lián)立方程組， 存在一組隨機變量Y1， …，Yn與自變量x1， …，xn之間滿足非線性關(guān)系， 如下：

其中：n個模型的同一次觀測模型誤差εi的各分量間是相關(guān)的， 即cov（εi）是非對角矩陣［17］。 這里多個方程之間存在聯(lián)系， 各方程的擾動項之間存在相關(guān)性， 同時估計多個方程能夠提高模型估計效率。

將樹冠外輪廓模型和樹冠體積模型兩兩聯(lián)立為方程組， 即2.1 中式（1）與式（5）、 式（2）與式（6）、 式（3）與式（7）、 式（4）與式（8）共4 個方程組。 每組的2 個方程擁有同一套參數(shù)， 對2 個模型共同進行擬合， 即解決非線性聯(lián)立方程組模型的參數(shù)估計問題。

為了保證參數(shù)估計的一致性和漸進無偏性， 利用SAS 統(tǒng)計軟件proc model 程序提供的似乎不相關(guān)回歸法（SUR）， 選擇SUR 法同時擬合樹冠外輪廓與體積相容性模型。 此外， 在林業(yè)模型擬合過程中， 模型誤差項之間可能存在異方差的問題， 本研究采用模型回歸函數(shù)自身作為權(quán)函數(shù)消除異方差［13］。

2.3 模型評價和檢驗指標(biāo)

模型擬合和檢驗結(jié)果通過以下指標(biāo)評價： 決定系數(shù)（R2）、 均方根誤差（RMSE）和平均絕對誤差（MAE）、平均偏差（MD）。 最優(yōu)模型選擇根據(jù)R2最大，RMSE、MAE、MD絕對值最小的原則進行。 具體公式如下：

其中：yi為第i因變量實際值， y＾i為第i因變量預(yù)測值，yi為因變量實際值平均，n為樣本數(shù)。

3 結(jié)果與分析

3.1 平均斷面積求積法

樹冠體積數(shù)據(jù)可采用分層切割法處理得到。 本研究在對杉木樹冠進行調(diào)查時， 采用平均斷面積求積法， 將樹冠從樹梢到樹冠基部按照相對冠長分6 部分（0～0.10， 0.10～0.25， 0.25～0.50， 0.50～0.75， 0.75～0.90， 0.90～1.00）。 將樹冠最上部分看作圓錐體近似計算體積， 其他各部分采用平均斷面積法。 計算公式如下：

其中：Vc表示整個樹冠體積（m3）；Vi為第i部分樹冠體積（m3）；gi為第i區(qū)分段中央斷面積（m2）；li為第i分段長度（m）；g′為梢頭底端斷面積（m2）；l′為梢頭樹冠長度（m）；n為分段個數(shù)。

3.2 樹冠外輪廓-體積一致性方程系統(tǒng)的總體評價

利用SAS 軟件的proc model 模塊中似乎不相關(guān)回歸過程SUR 方法對樹冠外輪廓模型和樹冠體積預(yù)測模型的一致性方程組同時進行擬合， 表2 給出了不同模型擬合的統(tǒng)計量， 即決定系數(shù)（R2）、 均方根誤差（RMSE）。

表2 樹冠外輪廓-體積模型擬合結(jié)果Table 2 Fitting results of crown profile and crown volume models

由表2 擬合結(jié)果可知： 樹冠外輪廓模型擬合指標(biāo)R2的對比結(jié)果為模型4＞模型3＞模型2＞模型1，輪廓模型擬合指標(biāo)RMSE的對比結(jié)果為模型4＜模型3＜模型2＜模型1， 樹冠體積模型擬合指標(biāo)R2的對比結(jié)果為模型4=模型3＞模型2＞模型1， 體積模型擬合指標(biāo)RMSE結(jié)果為模型4＜模型3＜模型2＜模型1。 根據(jù)R2最大且RMSE相對較小的最優(yōu)模型選擇標(biāo)準(zhǔn)， 不管從樹冠外輪廓還是體積模型結(jié)果來看， 模型4 顯示了較好的擬合結(jié)果， 最終利用模型4 來描述福建地區(qū)杉木樹冠外輪廓和樹冠體積。

3.3 相容模型檢驗

建模數(shù)據(jù)只能反映模型擬合的好壞， 不能反映模型的預(yù)測性能。 模型系統(tǒng)的獨立性檢驗是采用建模時未使用的獨立樣本數(shù)據(jù)， 對各模型系統(tǒng)的預(yù)測性能進行綜合評價。 基于表3 的參數(shù)估計值和檢驗數(shù)據(jù)， 利用SAS 軟件計算各模型系統(tǒng)樹冠外輪廓和體積的絕對誤差、 均方根誤差。 從表3 可以看出： 各樹冠外輪廓模型的檢驗指標(biāo)R2對比結(jié)果為模型4＞模型3＞模型2＞模型1， 外輪廓模型RMSE對比結(jié)果為模型4＜模型3＜模型2＜模型1， 各樹冠體積模型檢驗指標(biāo)R2對比結(jié)果為模型3=模型4＞模型2＞模型1， 體積模型RMSE對比結(jié)果為模型4＜模型3＜模型2＜模型1， 不管是從R2還是RMSE來看， 模型4均優(yōu)于其他模型。 此外， 模型4 的檢驗指標(biāo)MAE和MD絕對值最小， 真實值與預(yù)測值之間誤差最小， 進一步驗證了模型4 作為杉木樹冠外輪廓-體積模型的合理性。

3.4 模型殘差圖評價

模型擬合的總體評價反映了總體樹冠外輪廓和體積的變化， 不能反映各模型是否存在異方差性和無偏性， 評價這2 個指標(biāo)最直觀的方法就是利用殘差分布圖。 為全面評價模型4 效果， 分別繪制模型4 加權(quán)前后的殘差分布圖。 圖2 為未增加權(quán)函數(shù)時相容模型的殘差圖， 左側(cè)為樹冠外輪廓模型殘差分布， 右側(cè)為樹冠體積模型殘差分布， 從圖2 可以看出： 殘差均存在明顯喇叭口形狀， 說明異方差問題顯著。 圖3 為增加權(quán)函數(shù)后相容模型殘差分布圖， 從圖3 可以看出： 加權(quán)后模型的殘差散點圖分布變得均勻， 說明權(quán)函數(shù)明顯消除了異方差。 進一步說明加權(quán)后模型4 顯示了較高的等方差性和無偏性， 效果較好。

表3 不同模型獨立性檢驗Table 3 Validation results of different crown models

圖2 未加權(quán)的SUR 方法擬合的相容模型殘差分布圖Figure 2 Residual distribution of models fitted by unweighted SUR method

圖3 加權(quán)的SUR 方法擬合的相容模型殘差分布圖Figure 3 Residual distribution of models fitted by weighted SUR method

4 結(jié)論

在樹冠模型研究方面， 許多學(xué)者利用規(guī)則幾何體模擬樹冠形狀并計算樹冠體積， 這類方法具有簡單方便特點， 但是不靈活且預(yù)測精度相對較低。 本研究以福建杉木為研究對象， 選擇4 種常用的樹冠外輪廓經(jīng)驗?zāi)Ｐ湍M樹冠形狀曲線， 與簡單幾何體模擬方法相比， 能夠更加準(zhǔn)確合理地描述樹冠形狀的變化， 而且方程形式更加靈活。

本研究通過樹冠外輪廓模型方程推導(dǎo)樹冠體積方程， 構(gòu)建一致性的相容的非線性樹冠外輪廓-體積聯(lián)立方程組模型， 利用SAS 軟件模塊中的似乎不相關(guān)回歸過程（SUR）解決復(fù)雜分段聯(lián)立方程組模型系統(tǒng)的參數(shù)同時估計， 確保了2 個方程參數(shù)估計的一致性， 模型預(yù)測效果較好。 本研究構(gòu)建的樹冠外輪廓-體積一致性模型方程， 可以預(yù)測樹冠外輪廓形態(tài)， 驅(qū)動林分樹木三維可視化并預(yù)估樹冠體積， 實現(xiàn)了樹冠外輪廓與體積模型之間互相推導(dǎo)， 同時也為進一步估測樹木地上部分生物量提供了理論依據(jù)。