王杰芬,夏磊,林露花,胡璐璐,徐懷興,王聚中,徐小軍
(1. 浙江省龍泉市林業(yè)局,浙江 龍泉 323700;2. 浙江農(nóng)林大學(xué) 環(huán)境與資源學(xué)院,浙江 杭州 311300)
全球氣候變化給人類生產(chǎn)和發(fā)展帶來巨大的影響,開展森林碳循環(huán)研究對減緩氣候變化具有重要意義。森林生態(tài)系統(tǒng)是陸地生態(tài)系統(tǒng)中最大的碳庫,其年固碳量占整個陸地生態(tài)系統(tǒng)年固碳量的2/3,具有很強的碳儲存能力[1]。我國森林年固碳量約為4.34 億t,人工林對我國森林碳庫持續(xù)增長起到重要的作用[2],尤其是杉木Cunninghamialanceolata人工林種植面積大、分布范圍廣,其碳儲量約為3.03 億t,位居前列[3]。因此,開展杉木碳儲量估算研究,有利于準(zhǔn)確掌握杉木碳儲量現(xiàn)狀和動態(tài)變化,對提升杉木林碳儲量和固碳潛力具有重要的意義。
目前,森林碳儲量估算的方法主要包括生物量法、模型模擬法和遙感估算法等[4-5]。通過樹高、胸徑等測樹因子建立異速生長方程估算碳儲量,是最常用的方法。Mei 等用胸徑、樹高和木材密度建立了杉木個體生物量模型,對福建將樂國有林場的杉木林分進(jìn)行了生物量估算[6]。莊紅蕾等用異速生長方程法對崇明島水杉Metasequoiaglyptostroboides幼齡林、中齡林和成熟林進(jìn)行了碳儲量估算,取得較高的精度[7]。但通過樣地調(diào)查方法獲取樹高、冠幅和胸徑等測樹因子構(gòu)建異速生長方程,不僅費時費力、效率低,而且難以大面積應(yīng)用[8]。隨著無人機遙感技術(shù)的發(fā)展和普及,為樹高、冠幅和胸徑等參數(shù)的獲取提供了新方法[9]。徐永勝等用改進(jìn)的分水嶺分割算法提取單木冠幅和樹冠投影面積,精度分別為79.84%和76.04%[10]。王越等采用無人機傾斜攝影測量技術(shù)構(gòu)建三維模型,進(jìn)而提取樹高和冠幅[11]。謝巧雅等利用無人機影像構(gòu)建數(shù)字表面模型和數(shù)字高程模型,通過兩者相減,從而提取杉木樹高,精度達(dá)90.8%[12]。但是,對于郁閉度較大的林分,地面信息獲取難度很大,難以獲取精準(zhǔn)的數(shù)字地形模型,造成樹高提取誤差大[13]。另一方面,樹高提取通常需要單木樹冠頂點信息,當(dāng)樹冠頂點提取存在偏差時,會導(dǎo)致樹高低估的現(xiàn)象[12]。因此,通過無人機可見光影像提取山區(qū)高郁閉度林分樹高仍然存在較大的不確定性。冠幅信息在無人機可見光影像中能夠得到直觀的反映,且與樹高和胸徑有較強的相關(guān)性,冠幅越大,胸徑也越大[14-16]。張玉薇等將提取的冠幅和樹冠面積與胸徑建立回歸模型,能夠準(zhǔn)確地估算胸徑,相對誤差低于5%[17]。賈鵬剛等采用無人機影像提取銀杏Gingkobiloba單木樹冠面積、冠幅和樹高三個因子可較高精度地估算出胸徑,精度滿足森林資源調(diào)查要求[18]。Gonzalez-Benecke 等研究得出樹冠面積對提高單木胸徑和蓄積估算精度起到重要作用[15]。因此,以冠幅為自變量,可較準(zhǔn)確地估算單木胸徑,為基于無人機影像估算單木和林分水平杉木碳儲量提供可行性。
浙江省是杉木主要分布省份,麗水市是浙江省杉木主要分布區(qū)域。本研究借助大疆無人機M300 獲取可見光影像,對影像進(jìn)行正射校正后,采用放射線法提取單木冠幅,并對冠幅進(jìn)行精度驗證,然后,建立冠幅與胸徑以及胸徑和樹高之間的回歸模型,估算出胸徑和樹高;最后,將樹高和胸徑代入杉木碳儲量異速生長模型,估算出單木和林分水平碳儲量,并進(jìn)行精度評價。本研究成果將驗證基于無人機影像提取的冠幅信息在估算杉木胸徑、樹高以及單木和林分水平碳儲量上的可行性,為快速準(zhǔn)確估算杉木林碳儲量提供技術(shù)和方法支撐。
研究區(qū)位于浙江省西南部的麗水龍泉市,地理坐標(biāo)為27°42' ~ 28°20' N,118°42' ~ 119°25' E,東西長為70.25 km,南北寬為70.80 km,總面積為 304 486 hm2。屬中亞熱帶季風(fēng)氣候區(qū),四季分明、雨量充沛、溫暖濕潤。全市森林資源豐富,以針葉林、闊葉林、針闊混交林為主,林業(yè)用地面積為 26.6 萬hm2,活立木總蓄積量為1 987 萬m3,森林覆蓋率高達(dá)84.38%,位居浙江省前列[20]。全市杉木林面積為11.15 萬hm2,平均蓄積量為 97 m3·hm-2[19]。
在龍泉市查田鎮(zhèn)和八都鎮(zhèn)分別選擇1 個杉木純林樣地,林地面積分別為774.2 m2和920.4 m2。樣地1 位于查田鎮(zhèn),地理坐標(biāo)為118°57'56.8'' E,27°53'11.4'' N,為18 年生人工中齡林,初植密度為3 300 株·hm-2。樣地2位于八都鎮(zhèn),地理坐標(biāo)為118°57'18.8'' E,27°59'33.7'' N,為27 年生人工成熟林,初植密度為2 700 株·hm-2。2個樣地均開展間伐撫育、割灌除草、除藤等經(jīng)營措施。樣地1 和樣地2 的坡度分別為12°和5°。采用十步抬頭法,沿樣地兩條對角線記錄被枝葉遮擋次數(shù),計算樣地1 和樣地2 杉木林冠層郁閉度分別為0.8 和0.4。樣地1林下以土壤、枯枝為主,有少量蕨類植物;樣地2 林下植被茂密,覆被灌木、蕨類和草本植物。樣地大小參考浙江省森林資源連續(xù)清查標(biāo)準(zhǔn)樣地(28.28 m×28.28 m)[20],并將樣地劃分為9 個小樣方。為了能夠在無人機影像中準(zhǔn)確地找出每個小樣方的位置,在每個小樣方角點上豎立標(biāo)桿,標(biāo)桿上放置易識別的紅布。由于小樣方角點被杉木遮擋,導(dǎo)致在無人機上無法識別,因此將標(biāo)桿移至距離小樣方角點最近的林窗內(nèi),最終,將每個大樣地劃分為9 個不規(guī)則小樣方(圖1)。記錄每個小樣方內(nèi)的每株杉木坐標(biāo)并繪制杉木分布圖。對樣地內(nèi)的杉木進(jìn)行每木檢尺,起測徑階為5 cm,調(diào)查樣地內(nèi)每棵杉木的胸徑和冠幅等。另外,從森林資源二類調(diào)查數(shù)據(jù)中獲取查田鎮(zhèn)和八都鎮(zhèn)杉木小班平均胸徑和平均樹高,查田鎮(zhèn)有1 016 個杉木小班,八都鎮(zhèn)有1 360 個杉木小班,用于構(gòu)建樹高和胸徑關(guān)系式。
圖1 樣地和杉木分布圖Fig. 1 Sample plots and distribution of C. lanceolata
2022 年2 月,選擇晴朗無風(fēng)天氣,使用大疆無人機M300 搭載的禪思H20T 可見光傳感器獲取影像數(shù)據(jù),航拍高度設(shè)置為100 m。樣地1 拍攝時間在11:20—11:40,樣地2 拍攝時間在09:50—10:10。在2 個樣地分別拍攝16 張單幅影像,然后將無人機單幅影像導(dǎo)入Agisoft PhotoScan Pro 軟件做拼接和正射校正處理,生成正射影像,影像的空間分辨率為3.48 cm·pixel-1。
放射線法是單木冠幅估計的一種常用方法,該方法以樹冠中心點為起點,向四周發(fā)射射線,根據(jù)樹冠中心點亮度值大于樹冠邊緣點亮度值的前提假設(shè),即由中心點向四周延伸的射線經(jīng)過的像元亮度值呈遞減的趨勢,并且在樹冠邊緣處遞減的梯度最大,因此可以通過對射線亮度值曲線做多項式擬合平滑處理,然后找出平滑后曲線亮度值變化最大處所對應(yīng)的位置,即為樹冠的邊緣點,計算各邊緣點到樹冠中心點的距離就可得到冠幅[21]。本文基于放射線法的原理,在找尋樹冠邊緣點時,將各個方向射線上亮度值最小且與樹冠中心點距離最近的像元作為樹冠邊緣,并且加入了樹冠邊緣點位置調(diào)整系數(shù)(α)以及各射線冠幅估計值異常值剔除百分?jǐn)?shù)(β),來提高單木冠幅估算精度。具體步驟如下:
(1)以樹冠中心點(通過目視解譯方法獲?。橹行?,以N為窗口大小,構(gòu)建一個N×N大小的窗口。以窗口左上角像元為起點,依次按順時針方向提取經(jīng)過樹冠中心點的射線。對于一個N×N的窗口,共有2×N-2 條射線。窗口越大,射線數(shù)量越多,N通常稍大于影像中樹冠冠幅最大值。
(2)依次提取出各射線的特征曲線(亮度值)并進(jìn)行平滑處理。找出特征曲線的最大值及其對應(yīng)的像元,以該像元為分割點,分別搜索分割點兩邊離分割點最近的最小亮度值。然后,分別計算分割點兩邊亮度值的變幅,即最大值與最小值之差,查找兩邊亮度值與α×變幅+最小值最接近的像元,記為左側(cè)邊緣點和右側(cè)邊緣點,然后計算出左側(cè)邊緣點和右側(cè)邊緣點之間的距離,即為該射線得到的冠幅估計值。依次計算出每條射線的冠幅估計值。
(3)為了剔除由于樹冠內(nèi)大孔隙背景、樹冠周邊相似背景以及樹冠重疊等原因造成的過小或過大的冠幅估計值,計算出所有射線冠幅估計值的百分位數(shù)β和100-β對應(yīng)的冠幅估計值CLβ和CL100-β,取冠幅估計值位于CLβ與CL100-β之間的數(shù)值,并取均值作為最終的單木冠幅估計值。
首先,根據(jù)樣地調(diào)查中實測的378 棵杉木胸徑和冠幅數(shù)據(jù),建立胸徑和冠幅回歸模型;其次,根據(jù)小班杉木平均胸徑和平均樹高數(shù)據(jù),建立胸徑和樹高回歸模型;然后,將基于無人機提取的冠幅值代入胸徑和冠幅回歸模型,得到胸徑估計值,將胸徑估計值代入胸徑和樹高回歸模型,得到樹高估計值;最后,將胸徑和樹高估計值代入杉木全株碳儲量相對生長模型,計算出樣地內(nèi)每棵杉木的碳儲量。杉木全株生物量相對生長模型采用由劉雯雯等構(gòu)建的區(qū)域尺度杉木生物量通用相對生長方程,見式(1)[22]。該方程適用胸徑和樹高范圍廣、通用性強,決定系數(shù)(R2)為0.93。杉木生物量轉(zhuǎn)換為碳儲量的系數(shù)為0.5[23-24]。
式中,W為全株生物量;D和H分別為胸徑和樹高。
根據(jù)每棵杉木的坐標(biāo)及繪制的杉木位置分布圖,在無人機影像中標(biāo)記出每棵杉木,確保每棵杉木實測冠幅與無人機影像提取冠幅一一對應(yīng)。采用實測的冠幅對無人機影像提取的冠幅進(jìn)行精度評價。采用樣地實測胸徑對胸徑估計值進(jìn)行精度評價。將實測胸徑和樹高估計值代入杉木碳儲量相對生長模型計算出碳儲量作為參考值,用于評價基于估算的胸徑和樹高得到的杉木碳儲量估計值的精度。評價指標(biāo)包括R2、均方根誤差(RMSE)、相對均方根誤差(RMSEr)、相對誤差(RE)和偏置(Bias)。
樣地1 和樣地2 的小樣方面積及杉木株數(shù)、胸徑和冠幅情況見表1。由表1 可知,樣地1 和樣地2 的面積分別為774.2 m2和920.4 m2,兩個樣地的小樣方面積分別在64.0 ~ 121.2 m2和77.3 ~ 122.9 m2之間。樣地1 共有257 棵杉木,其平均胸徑和冠幅分別為13.12 cm 和2.16 m,小樣方中的杉木株數(shù)在14 ~ 37 棵之間,冠幅在2.36~ 2.96 m 之間,胸徑在11.97 ~ 15.36 cm 之間。樣地2 共有121 棵杉木,其平均胸徑和冠幅分別為15.98 cm 和2.77 m,小樣方中的杉木株數(shù)在9 ~ 21 棵之間,胸徑在15.01 ~ 17.94 cm 之間,冠幅在2.23 ~ 3.11 m 之間。樣地1 中的株數(shù)要大于樣地2,樣地2 中杉木的胸徑和冠幅大于樣地1 的。
表1 樣地內(nèi)杉木每木檢尺調(diào)查數(shù)據(jù)Tab. 1 Tally of the sample plots
2.2.1 胸徑與樹高建模 利用2016 年森林資源二類調(diào)查小班杉木的平均胸徑與平均樹高,建立以胸徑為自變量,樹高為因變量的回歸模型,兩者散點圖和表達(dá)式見圖2。由圖可知,樣地1 和樣地2 中杉木胸徑和樹高回歸關(guān)系符合線性函數(shù),其擬合決定系數(shù)R2分別為0.91 和0.89。對樣地1 和樣地2 中杉木胸徑和樹高關(guān)系式進(jìn)行交叉驗證得出,樣地1 中杉木樹高實測值與估計值之間的RMSE=0.79 m、RMSEr=11.29%、Bias=0.28 m;樣地2 中杉木樹高實測值與估計值之間RMSE=0.82 m、RMSEr=11.80%、Bias=-0.30 m。樣地1 和樣地2 中杉木胸徑和樹高回歸模型精度分別為88.71%和88.20%,均能準(zhǔn)確地預(yù)估杉木樹高。
圖2 胸徑與樹高模型構(gòu)建與評價Fig. 2 Model establishment and evaluation between DBH and height
2.2.2 冠幅與胸徑建模 利用實測的杉木冠幅和胸徑數(shù)據(jù),建立冠幅與胸徑之間的回歸模型,其散點圖和表達(dá)式見圖3A。由圖可知,樣地1 和樣地2 中杉木的冠幅與胸徑之間都存在較好的線性關(guān)系,擬合決定系數(shù)R2分別為0.55 和0.61。對樣地1 和樣地2 中杉木胸徑和樹高關(guān)系式進(jìn)行交叉驗證得出,樣地1 中杉木實測胸徑和估計胸徑之間的RMSE=3.29 cm,RMSEr=25.41%;樣地2 中杉木實測胸徑和估計胸徑之間的RMSE=3.31 cm,RMSEr=20.90%。樣地1 和樣地2 中杉木采用實測冠幅估算胸徑的精度分別為74.59%和79.10%。
圖3 冠幅與胸徑模型構(gòu)建與評價Fig. 3 Model establishment and evaluation between crown width and DBH
以無人機影像綠光波段亮度值為特征值,采用放射線法提取樣地1 和樣地2 中杉木的冠幅。參考樣地實測杉木冠幅最大值,樣地1 和樣地2 的窗口大小N分別為4.2 m 和4.9 m。樹冠邊緣點位置調(diào)整系數(shù)(α)和異常值剔除百分?jǐn)?shù)(β)分別為0.1 和10%。分別以提取冠幅和實測冠幅為圓直徑繪制杉木樹冠邊界,兩個樣地中杉木的冠幅提取結(jié)果見圖4。從圖4 可見,放射線法提取的杉木樹冠邊界與實測樹冠邊界能夠較好地重合,但存在少部分樹冠邊界被明顯放大和縮小問題。
圖4 提取冠幅與實測冠幅對比圖Fig. 4 Comparison on detected crown width and measured one
杉木的提取冠幅與實測冠幅之間存在顯著(P<0.01)相關(guān)性,決定系數(shù)R2=0.22(圖5A)。但大多數(shù)樣本點聚集在1∶1 線附近,說明大部分樹冠的冠幅能夠被準(zhǔn)確地提取,但是也存在較小的冠幅被高估和較大的冠幅被低估的情況。樣地1 中杉木提取冠幅與實測冠幅之間的RMSE、RMSEr和Bias分別為0.63 m、24.10%和0.09 m;樣地2 中杉木提取冠幅與實測冠幅之間的RMSE、RMSEr和Bias分別為0.79 m、28.45%和-0.22 m。從Bias可得出,樣地1 中杉木的提取冠幅與實測冠幅之間系統(tǒng)誤差很小,說明杉木的提取冠幅均值能夠準(zhǔn)確反映實測冠幅均值,但樣地2 中杉木的提取冠幅存在系統(tǒng)高估。
圖5 杉木測樹因子實測值與估計值之間散點圖Fig. 5 Scatter plot of measured and estimated growth factors
將無人機影像提取的杉木冠幅代入杉木冠幅與胸徑關(guān)系模型估算出胸徑。由估計胸徑與實測胸徑對比結(jié)果表明(圖5B),估計胸徑和實測胸徑之間R2= 0.36(P<0.01),大多數(shù)樣本靠近1∶1 線,表明能夠得到準(zhǔn)確估算。樣地1 中杉木估計胸徑與實測胸徑之間RMSE=2.85 cm、RMSEr=21.76%、Bias=0.27 cm,胸徑估計值系統(tǒng)偏?。粯拥? 中杉木估計胸徑與實測胸徑的RMSE=3.47 cm、RMSEr=21.69%、Bias=-1.05,胸徑估計值系統(tǒng)偏大。將杉木估計胸徑代入胸徑和樹高關(guān)系模型進(jìn)行樹高估算,杉木估計樹高與參考樹高對比結(jié)果表明(圖5A),估計樹高和參考樹高之間的R2=0.32(P<0.01),大多數(shù)樣本能夠得到準(zhǔn)確估算。樣地1 中杉木估計樹高與參考樹高的RMSE=1.79 m、RMSEr=20.16%、Bias=0.17 m,樹高估計值系統(tǒng)偏??;樣地2 中杉木估計樹高與參考樹高的RMSE=2.06 m、RMSEr=20.20%、Bias=-0.63 m,樹高估計值系統(tǒng)偏大。
將基于無人機影像提取的冠幅估計的胸徑和樹高代入杉木碳儲量相對生長模型估計單株碳儲量,并與參考的單株碳儲量進(jìn)行比較,結(jié)果如圖6 所示。由圖可知,樣地1單株杉木碳儲量估計值和參考值分別在3.92 ~ 35.49 kg 之間和4.46 ~ 86.53 kg 之間,兩者平均值分別為19.27 kg 和21.46 kg。樣地2 單株杉木碳儲量估計值和參考值分別在15.07 ~ 53.72 kg 之間和5.55 ~ 78.36 kg 之間,兩者平均值分別為34.40 kg 和31.50 kg。杉木的估計單株碳儲量與參考單株碳儲量的線性擬合R2為0.31(P<0.01),單株杉木碳儲量較大的值存在明顯的低估,但大部分樣本靠近1∶1線,說明大部分樣本能夠被準(zhǔn)確估計。杉木的估計單株碳儲量與參考單株碳儲量精度分析表明,樣地1 中杉木單株碳儲量估計值與單株碳儲量參考值之間的RMSE=11.59 kg、RMSEr=54.00%、Bias=2.20 kg;樣地2 中的RMSE=14.30 kg、RMSEr=45.40%、Bias=-2.90 kg,估計值系統(tǒng)偏大。樣地1 和樣地2 中的杉木單株碳儲量估計精度分別為46.00%和54.60%,精度相對較低。
圖6 單株碳儲量估計值與參考值散點圖Fig. 6 Scatter plot of estimated and reference carbon storage of individual tree
將兩個樣地內(nèi)單株杉木碳儲量累加得到樣地1 碳儲量估計值和參考值分別為63.98 t·hm2和71.27 t·hm2,樣地2 碳儲量估計值和參考值分別為45.23 t·hm2和41.42 t·hm2。由樣地尺度杉木測樹因子和碳儲量估計值與參考值的精度比較得出(表2),樣地1 的9 個小樣方中杉木的平均冠幅、胸徑、樹高和碳儲量估計值的RE變化幅度分別在-8.64% ~ 13.64%、-9.19% ~ 15.22%、-8.45% ~ 14.25%和-16.80% ~ 33.35%之間。樣地2 的 9 個小樣方中杉木的平均冠幅、胸徑、樹高和碳儲量估計值的RE變化幅度分別在-34.90% ~ 5.50%、-35.11% ~ 5.58%、-32.13% ~ 5.23%和-80.52% ~ 13.16%之間。整個樣地的平均冠幅、胸徑、樹高和碳儲量估計值的RE都低于11%。樣地2 的第9 個小樣方中的杉木冠幅估計值明顯高估,導(dǎo)致其他因子估計精度也較低。相對單株尺度精度而言,在小樣方和整個樣地水平上的各測樹因子和碳儲量估計值精度都得到明顯提高。樣地1 和樣地2 水平各測樹因子及碳儲量精度分別在89.77%和90.79%以上。
表2 樣地水平測樹因子和碳儲量估計值相對誤差Tab. 2 Relative errors of the growth factors and carbon storage estimation at the sample level
(1)結(jié)合無人機可見光影像和放射線法提取的杉木林中的杉木冠幅與實測冠幅之間存在較高的誤差,主要原因是實測冠幅采用自下而上觀測方式,受冠層重疊遮擋影響小,而基于無人機影像自上而下觀測樹冠,受冠層重疊遮擋影響較大,往往造成提取冠幅小于實測值[17]。另外,受樹冠內(nèi)孔隙陰影等影響,使得樹冠內(nèi)陰影像元被誤判為樹冠邊緣點,造成提取冠幅小于實測值[25]。受林下植被等影響,樹冠邊緣像元與背景像元非常接近,將背景像元誤判為樹冠邊緣點,使得提取冠幅高于實測值。
(2)單木碳儲量估計精度較低的主要原因是基于無人機提取的單木冠幅、胸徑和樹高精度還不夠準(zhǔn)確,存在高值低估和低值高估的問題。此外,僅采用冠幅一個變量估算胸徑也存在不足,研究表明以冠幅和樹冠面積為自變量能夠提高胸徑估算精度[17]。胸徑和樹高之間的關(guān)系受立地條件、林分條件和競爭等影響,不同區(qū)域或樣地之間具有差異性[26]。構(gòu)建含林分變量或基于混合效應(yīng)的樹高和胸徑關(guān)系式,能提高樹高預(yù)測精度[26]。對于整個樣地,高估和低估的冠幅估計值相互抵消后,樣地的平均冠幅、胸徑和樹高估計值精度明顯提高,進(jìn)而樣地碳儲量精度也較高。本研究樹冠中心點采用目視解譯方法獲取,未考慮株數(shù)提取誤差對樣地碳儲量估計精度影響。
本研究利用無人機可見光影像,采用放射線法實現(xiàn)了杉木單木冠幅自動提取,然后基于無人機影像提取的冠幅數(shù)據(jù),結(jié)合樣地調(diào)查數(shù)據(jù)建立的冠幅-胸徑以及胸徑-樹高回歸模型,估計出單木胸徑和樹高,將兩者代入杉木碳儲量相對生長模型估計單株和樣地碳儲量。由上述方法得到的單木冠幅、胸徑、樹高和碳儲量精度相對偏低,存在高估和低估問題,說明本研究方法在估算單木尺度杉木結(jié)構(gòu)參數(shù)和碳儲量上精度還有待提高。但將尺度上升到與浙江省森林資源連續(xù)清查標(biāo)準(zhǔn)樣地面積時,樣地平均冠幅、胸徑、樹高和碳儲量估計精度在89%以上。本研究方法能夠滿足樣地尺度的冠幅、胸徑、樹高和碳儲量估算精度要求,具有較高的實用性。