楊曉慧 吳金卓 劉浩然 鐘 浩 林文樹
(1. 東北林業(yè)大學機電工程學院 哈爾濱 150040;2. 東北林業(yè)大學土木與交通學院 哈爾濱 150040)
作為森林生態(tài)系統(tǒng)的重要組成部分,人工林不僅改善了我國森林資源短缺的問題,而且也為國民經(jīng)濟建設提供了大量木材和林副產(chǎn)品,在發(fā)揮經(jīng)濟效益的同時對維持生態(tài)系統(tǒng)平衡、減緩全球氣候變暖等起著重要作用(銀彬吾等,2019)。 冠層是森林與外界環(huán)境相互作用最直接和最活躍的界面層,是植物進行呼吸作用、蒸騰作用和光合作用等生理過程的主要場所,本身承載了森林生物多樣性的主體部分。郁閉度是森林冠層的主要參數(shù),是反映森林結(jié)構(gòu)和森林環(huán)境的重要因子,在水土流失、水源涵養(yǎng)、林分質(zhì)量評價、森林景觀建設等方面有著廣泛應用。
林分郁閉度的傳統(tǒng)測定方法主要有目測法、樹冠投影法、樣點法、魚眼照片法和光學遙感法等,其中,地面實測方法獲取的數(shù)據(jù)通常只在小范圍內(nèi)具有代表性且測量過程耗時耗力(孫釗等,2020),光學遙感因穿透性低一定程度上影響森林垂直結(jié)構(gòu)參數(shù)估測。激光雷達(light detection and ranging,LiDAR)是一種新興的主動遙感技術(shù),能夠有效穿透森林,在多時空尺度上獲取森林生態(tài)系統(tǒng)高分辨率的三維地形、植被結(jié)構(gòu)參數(shù)、葉面積指數(shù)等參數(shù)。Lefsky 等(2002)提出的地基激光雷達技術(shù)(terrestrial laser scanning,TLS)可用于測量單木三維結(jié)構(gòu)和植被冠層結(jié)構(gòu);機載激光雷達技術(shù)(airborne laser scanning,ALS)的出現(xiàn)使得提取高精度森林結(jié)構(gòu)參數(shù)成為可能(瞿帥等,2018);星載激光雷達覆蓋范圍廣,可用于大區(qū)域林業(yè)資源調(diào)查研究。無人機激光雷達相比上述3 種平臺具有靈活性高、成本低、效率高、受地形影響小、影像重疊率大、分辨率高、姿態(tài)角大、相幅小、不受云霧干擾等特點(王娟等,2020;Salamíet al., 2014),近年來在森林結(jié)構(gòu)參數(shù)提取中得到越來越多的應用(Davieset al., 2014)。
隨著遙感技術(shù)的不斷發(fā)展,國內(nèi)外已有利用激光雷達數(shù)據(jù)進行林分郁閉度估測的研究報道。李丹(2012)基于配準后的多站地基激光雷達數(shù)據(jù),設定不同高度閾值提取喬木樣地郁閉度,估測精度達83.91%。趙勛等(2020)以機載 LiDAR 點云數(shù)據(jù)為基礎,采用首次回波歸一化后 2 m 以上點云密度與所有點云密度之比估測郁閉度,估測值與實測值之間的R2達0.478。Arum?e 等(2018)利用機載激光掃描(ALS)數(shù)據(jù),設定不同高度閾值和回波計算郁閉度,并與基于數(shù)字半球照片(digital hemispherical photographs,DHP)和Crown模型(CCRCrown)的估測值進行了比較。García 等(2012)采用地球科學測高系統(tǒng)(geoscience laser altimetry system,GLAS)數(shù)據(jù)中的冠層基高度(canopy base height,CBH)估算林分郁閉度,精度達89%。吳項乾(2018)結(jié)合單木分割方法,先由UAV-LiDAR 點云生成冠層高度模型(canopy height model,CHM),再將基于點云提取的特征變量逐步引入郁閉度估測模型,R2max達0.78。汪霖等(2019)以地面實測數(shù)據(jù)、無人機高分影像、激光雷達點云數(shù)據(jù)為主要信息源,采用局部最大值法和種子點分割法提取林分郁閉度并對其精度進行了檢驗。Aicardi 等(2017)將TLS 點云與UAV影像點云數(shù)據(jù)結(jié)合,建立了一個具有高分辨率和高精度的森林數(shù)據(jù)集提取高郁閉度森林區(qū)域關(guān)鍵參數(shù)。
綜上可知,目前林分郁閉度估測大多基于地基或機載雷達數(shù)據(jù)進行,關(guān)于東北地區(qū)人工林林分郁閉度的研究較少,特別是林分高度、強度、冠層特征變量對針、闊葉林郁閉度反演過程貢獻率方面的研究鮮見報道。鑒于此,本研究選取我國東北地區(qū)具有代表性的興安落葉松(Larix gmelinii)、 樟子松(Pinus sylvestrisvar.mongolica)、 黑 皮 油 松 (Pinus tabuliformis)、 蒙古櫟(Quercus mongolica)、 白樺(Betula platyphylla)為研究對象,基于UAV-LiDAR 獲取的高密度點云數(shù)據(jù),結(jié)合植物冠層分析儀獲取的40 個樣本郁閉度實測值,采用30 個樣本應用統(tǒng)計模型法進行針葉林和闊葉林林分郁閉度建模,采用10個樣本驗證其精度,比較影響針、闊葉林郁閉度的雷達特征變量并進行郁閉度制圖,以期為快速準確估測人工林林分郁閉度提供基礎數(shù)據(jù)和技術(shù)參考。
東北林業(yè)大學城市林業(yè)示范基地(127°35′—127°39′E,45°42′—45°44′N,海拔 136~140 m)地處黑龍江省哈爾濱市中心城區(qū),占地面積 44 hm2,屬中溫帶大陸性季風氣候,年降雨量400~800 mm,年均氣溫3.5 ℃(劉浩然等,2021)?;貎?nèi)栽種有興安落葉松、樟子松、黑皮油松、胡桃楸(Juglans mandshurica)、蒙古櫟、水曲柳(Fraxinus mandshurica)、白樺等針葉和闊葉樹種,目前已形成小面積純林塊狀混交?;氐匦纹露绕骄彛帜据^為密集,郁閉度較高(范偉偉等,2020)。研究區(qū)地理位置及樣地空間分布見圖1。
圖1 研究區(qū)地理位置及樣地空間分布Fig. 1 Geographical location of the study area and spatial distribution of the sample plots
2.1.1 地面數(shù)據(jù)采集 2021 年5 月,在研究區(qū)選取4塊針葉林樣地,包括2 塊興安落葉松樣地、1 塊黑皮油松樣地、1 塊樟子松樣地;2 塊闊葉林樣地,包括1 塊白樺樣地、1 塊蒙古櫟樣地(表1)。采用WinSCANOPY 2010a 植物冠層分析儀獲取林木冠層半球圖像,得到各樣地實測郁閉度。具體做法是在針葉林和闊葉林樣地內(nèi)分別隨機建立20 個半徑為10 m 的圓形樣方,利用 GPS 記錄各樣方地理位置,坐標參考WGS-84。在光線均勻、避免直射的天空條件下于各樣方內(nèi)隨機選取3 個拍攝點,通過WinSCANOPY 2010a 植物冠層分析儀配有三星NV3 型號相機的180°魚眼鏡頭拍攝林分冠層獲取其半球圖像,魚眼鏡頭豎直向上并與水平面垂直。為消除灌木和雜草的影響,魚眼鏡頭離地高度設為1.5 m。
表1 研究樣地內(nèi)5 種人工林概況Tab. 1 Basic information of five kinds of plantation in the sample plots
2.1.2 激光雷達數(shù)據(jù)采集 2021 年5 月,使用四軸大疆經(jīng)緯M300 RTK 無人機搭載禪思L1 激光掃描儀獲取激光點云數(shù)據(jù)。無人機載激光雷達數(shù)據(jù)獲取參數(shù)如下:飛行高度100 m,飛行速度5 m·s-1,飛行間隔60 m,旁向重疊50%,主航線角度25°,平均點密度約197 pts·m-2,采樣頻率160 kHz,采取實時真彩色點云上色模式。
2.1.3 數(shù)據(jù)預處理 首先,對LiDAR 原始點云數(shù)據(jù)進行去噪、裁剪處理,通過CloudCompare 軟件中的布料模擬濾波(cloth simulation filter, CSF)方法獲得地面文件,并將分類的地面點插值以生成數(shù)字高程模型(digital elevation model,DEM)、數(shù)字表面模型(digital surface model,DSM),利用DEM 進行高度(z)歸一化獲得高度歸一化點云數(shù)據(jù)(艾薩迪拉·玉蘇甫等,2020;尤號田等,2020)。然后,應用Matlab 平臺構(gòu)建冠層體積模型(canopy volume model,CVM)計算冠層體積分布,采用PCM 點云數(shù)據(jù)處理軟件中基于CHM 的方法進行單木分割,分割后的點云導入CloudCompare 軟件中進行樹冠點云提取,以計算冠層密度、冠層回波的偏度和峭度等。對采集的樹冠球形真彩色圖像,將其轉(zhuǎn)換為灰度圖像并繪制灰度圖像直方圖,結(jié)合樹干剔除效果,調(diào)整分割閾值以確定合適閾值點(圖2),取3個拍攝點的平均值作為該樣方的郁閉度實測值。各樣方郁閉度測量值見表2。
表2 各樣方郁閉度測量值Tab. 2 Measured value of canopy closure of each quadrate
圖2 圖像二值化Fig. 2 Image binarization
首先基于點云三維坐標獲得冠層高度模型CHM和樹冠點云,由CHM、回波能量值、樹冠點云計算高度、強度、冠層3 類特征變量,然后采用主成分分析法和逐步回歸法將處理后的特征變量與植物冠層分析儀獲取的郁閉度進行建模,最后評價模型精度和適應性并進行郁閉度制圖(圖3)。
圖3 技術(shù)流程Fig. 3 Technical flow chart
2.2.1 LiDAR 特征變量提取與標準化 離散點云歸一化后,計算樣地尺度的激光雷達特征變量,包括樹木高度和強度特征變量;應用Matlab 平臺構(gòu)建冠層體積模型,計算冠層體積分布;基于單木分割提取樣地冠層信息,計算冠層密度、冠層回波的偏度和峭度等。冠層體積分布(canopy volume distribution, CVD)是描述三維點云形態(tài)的重要指標,通過對每塊樣地的點云集合進行三維分析,獲取基于三維體的冠層體積分布點云變量,步驟如下:按照三維坐標,將點云分為若干個5 m×5 m×0.5 m 的小三維體,若小三維體含有點云,則該小三維體屬于“非空三維體”,反之為“空三維體”;將位于冠層以上的“空三維體”定義為“孔隙體”(Open),反之定義為“封閉體”(Close);將位于上65%區(qū)域的“非空三維體”定義為“透光體”(Euphotic),反之定義為“遮光體”(Oligophotic);匯總所有三維體,形成冠層體積分布圖并匯總計算CVD 內(nèi)4 種類別三維體分別占有的體積比例(以某針葉樣本為例,見圖4)??紤]到地面存在干擾點,濾除0.3 m 以內(nèi)的地面點(Voseet al., 1995)。本研究計算高度、強度、冠層3 類特征變量(表3) 。利用IBM SPSS Statistics 軟件對提取完成的3 類特征變量進行標準化處理。
表3 LiDAR 特征變量及其描述①Tab. 3 LiDAR characteristic variables and descriptions
圖4 冠層體積模型Fig. 4 Canopy volume model
2.2.2 主成分分析 本研究共計算30 組×26 個特征變量,在每組特征變量內(nèi)部,變量之間可能存在相關(guān)性,若對每個指標分別進行分析,則不能完全利用數(shù)據(jù)中的信息,因此采用IBM SPSS Statistics 軟件,在保證主成分分析有效的前提下(KMO 檢驗系數(shù)>0.5,巴特利球體檢驗P<0.05),分組對高度、強度、冠層特征變量進行分析,減少分析指標的同時,盡可能多地保留原始變量信息且使變量彼此不相關(guān),作為新的綜合指標。主成分選取上,每組保留特征值大于1 的主成分。主成分載荷矩陣Ui與因子載荷矩陣Ai和特征值λi之間的關(guān)系為:
將Ui與每組特征變量的標準化值Zxi相乘即得主成分F的表達式:
2.2.3 逐步回歸分析 針對針葉林和闊葉林,利用提取的多組LiDAR 特征變量結(jié)合多元逐步回歸方法構(gòu)建郁閉度估測模型,形式如下(吳項乾等,2020):
式中:y為因變量;q為自變量;b為常數(shù);ε 為誤差。
以地面實測林分郁閉度為因變量,以高度、強度和冠層特征變量為自變量,應用逐步進入法(stepwise)將激光雷達特征變量按照與實測林分郁閉度相關(guān)性由高到低的順序引入方程,在保證每一步引入變量檢驗值F達到顯著水平(P<0.05)的前提下,篩選出高相關(guān)性變量,剔除低相關(guān)性變量,以得到最優(yōu)模型。
2.2.4 回歸模型精度評價 選取精度和真實程度高的變量組合作為反演模型。本研究采用反映自變量與因變量之間相關(guān)性的判定系數(shù)R2、調(diào)整判定系數(shù)AdjR2、均方根誤差RMSE(root mean square error)和相對均方根誤差rRMSE(relative root mean square error)評價回歸模型精度(洪奕豐等,2019)。R2越高,則二者之間相關(guān)性越高;但樣本容量一定時,增加解釋變量必定使自由度減少,AdjR2可在模型的復雜程度和衡量模型的優(yōu)良程度上取一個平衡,使模型趨于簡單;RMSE 反映實測值與預測值之間標準誤差大小,與評價對象本身的數(shù)值關(guān)系很大,RMSE 越小,則模型預測效果越好;rRMSE 為 RMSE 與預測值算數(shù)平均值的比值,反映模型總體預測精度,rRMSE 越小,表明模型預測精度越高(劉浩等,2018)。各評價指標計算公式如下:
式中:n為樣地數(shù)量;X為魚眼鏡頭測得的林分郁閉度;為實測郁閉度的平均值;X?為模型估測的林分郁閉度;k為自變量個數(shù)。
本研究提取高度、強度和冠層特征變量并將3 類特征變量標準化。經(jīng)主成分分析,針葉林樣本的高度特征變量得到3 個主成分F1c、F2c、F3c,強度特征變量得到4 個主成分F4c、F5c、F6c、F7c,冠層特征變量得到2 個主成分F8c、F9c;闊葉林樣本的高度、強度和冠層特征變量分別得到3、3 和2 個主成分,分別記為F1b、F2b…F8b。由式(2)得出各組特征變量的主成分得分及排名,以針葉林的高度特征變量得到的3 個主成分為例,X、Z軸分別代表第一和第二主成分,Y軸代表第三主成分,顏色不同的3 類點在空間中的位置表示興安落葉松、黑皮油松、樟子松3 類針葉樹種的第一、第二和第三主成分得分,結(jié)果見圖5。
圖5 針葉林高度特征變量主成分得分Fig. 5 Principal component score of height characteristic variables for coniferous forests
經(jīng)逐步回歸法分析,得到針葉林(式8)、闊葉林(式9)林分郁閉度估測的回歸方程:
回歸方程的精度如表4 所示。在針葉林郁閉度估測中,回歸方程的AdjR2=0.722,模型擬合優(yōu)度較高,顯著水平P=0.000 6<0.05,模型顯著性較高,變量F8c前的系數(shù)絕對值最大,表明冠層特征變量在針葉林郁閉度預測中貢獻最大;在闊葉林郁閉度估測中,回歸模型的AdjR2=0.725,模型擬合優(yōu)度高,P=0.000 6<0.05,模型顯著,強度特征變量對闊葉林郁閉度預測精度影響最顯著。
表4 回歸方程的精度Tab. 4 Accuracy of regression equation
將提取的特征變量進行主成分分析,代入回歸方程可得郁閉度估測值。針、闊葉林各選擇5 個樣方作為驗證樣本,通過回歸方程預測的林分郁閉度與WinSCANOPY 2010a 植物冠層分析儀獲得的郁閉度的交叉驗證結(jié)果見圖6。針葉林與地面實測郁閉度建立的預測模型精度為R2=0.72,闊葉林與地面實測郁閉度建立的預測模型精度為R2=0.77,篩選10 個檢驗點的郁閉度與實測郁閉度顯示出較高相關(guān)性(r=0.859)。通過反演模型計算得到的林分郁閉度與實測郁閉度較均勻地分布在擬合線兩側(cè),表明模型可用來反演研究區(qū)典型針、闊葉林林分郁閉度,反演結(jié)果精度較高。
圖6 林分郁閉度預測模型交叉驗證結(jié)果Fig. 6 Cross validation results of forest canopy closure prediction model
利用ArcGIS 的反距離權(quán)重插值法將激光雷達估測郁閉度與地面實測郁閉度可視化。由圖7 可知,郁閉度范圍在0.81~0.87 之間,樣地郁閉度呈東北西南高、西北東南低的趨勢,與每塊矩形樣地的走向基本一致。無人機激光雷達數(shù)據(jù)在預測林分郁閉度方面優(yōu)勢明顯,多變量模型能充分發(fā)掘森冠層信息,結(jié)合多組UAV-LiDAR 特征變量估測林分郁閉度能夠取得較好效果。
圖7 實測與估測森林郁閉度Fig. 7 Measurement and estimation of forest canopy closure
本研究基于無人機激光雷達點云數(shù)據(jù)和WinSCANOPY 2010a 植物冠層分析儀樣地實測數(shù)據(jù),采用逐步回歸法構(gòu)建人工林林分郁閉度估測模型,基本滿足郁閉度反演需要,可實現(xiàn)基于無人機激光雷達點云數(shù)據(jù)的林分郁閉度快速估測。谷金英等(2014)提取遙感光譜因子和地形因子,基于多元回歸分析法對林地郁閉度進行反演,精度達81.6%;吳飏等(2012)利用Spot5 數(shù)據(jù)結(jié)合96 個紋理特征估測林分郁閉度,主成分分析與逐步回歸分析結(jié)果顯示,近紅外波段的光譜和紋理特征與郁閉度相關(guān)性最顯著,加入紋理特征后估測精度可達84.32%;李擎等(2019)基于GF-2號遙感影像,采用主成分法進行天山云杉(Picea schrenkiana)林郁閉度估測研究,以紋理特征+光譜信息+地形因子為自變量構(gòu)建的估測模型擬合度為0.823。本研究在點云高度、冠層特征變量基礎上引入回波強度特征變量,針葉林回歸模型R2為0.782,闊葉林回歸模型R2為0.784。這是因為針葉樹冠形為圓錐形,無人機激光雷達在發(fā)射激光束的過程中易錯失冠頂,而闊葉樹冠頂呈橢球形,截獲點云概率較高,相對于針葉樹提取冠頂?shù)母怕矢?,高度特征變量的估測精度更高;此外,激光點云穿透針葉林冠層以及通過冠層間的空隙進入林內(nèi)形成地面回波點的概率更高,在植被點與全部回波點的比值運算中更易造成低估現(xiàn)象。因此,闊葉林郁閉度反演精度高于針葉林。由于本研究將對郁閉度影響較大的CC2m歸入冠層特征變量,故估測方程中冠層特征變量超過高度特征變量,成為對郁閉度預測精度影響更為顯著的因素。
與基于機載雷達點云數(shù)據(jù)進行郁閉度估測的研究相比,本研究數(shù)據(jù)獲取方式更靈活,時效性更強,精度更高,如張瑞英等(2016)結(jié)合ALS 點云數(shù)據(jù)與Landsat ETM+數(shù)據(jù),采用Cubist 模型估測溫帶森林郁閉度,R2為0.722;與基于無人機激光雷達點云提取的特征變量構(gòu)建統(tǒng)計模型進行郁閉度估測的研究相比,本研究估測精度有一定提高,如吳項乾(2018)采用雷達高度特征變量和冠層特征變量建模估測闊葉林郁閉度的R2為0.69(rRMSE=6.5%),而當通過高度特征變量和密度特征變量聯(lián)合進行郁閉度估測時R2為0.78(rRMSE=5.9%);與使用某一高度以上點云所占比例直接代表郁閉度的研究相比,本研究構(gòu)建模型在不同森林類型郁閉度估測中穩(wěn)定性更好,估測精度也有所提高,如穆喜云等(2015)采用高于2 m 的植被點回波與樣地內(nèi)全部點云回波的比值代表郁閉度,全部60 塊樣地的R2為0.721,闊葉林估測精度為56.62%。因此,本研究基于無人機激光雷達數(shù)據(jù)提取高度、冠層和回波強度特征變量構(gòu)建的人工林林分郁閉度估測模型具有較高精度,可滿足森林資源調(diào)查中的精度要求。
傳統(tǒng)人工調(diào)查方法實現(xiàn)森林資源動態(tài)監(jiān)測費時費力,針對樣地或區(qū)域等尺度的林木調(diào)查,利用無人機LiDAR 掃描是一種靈活、高效和低成本的地面數(shù)據(jù)獲取方法,能夠滿足快速的林分郁閉度估測需求。本研究提出的郁閉度估測方法只需獲取目標區(qū)域較高密度和較強穿透力的LiDAR 點云及少量地面實測驗證數(shù)據(jù),提取相關(guān)特征變量就能建立適配研究區(qū)森林特征的反演模型,進而得到精度較高的林分郁閉度估測值。本研究建立的人工針/闊葉林林分郁閉度模型,預估精度能達到森林資源調(diào)查相關(guān)技術(shù)規(guī)定要求,可在實踐中推廣應用。
需要說明的是,盡管無人機激光雷達在森林結(jié)構(gòu)參數(shù)提取和估測中具有諸多優(yōu)勢,但本研究還存在一些不足之處:因飛行數(shù)據(jù)無法同時保證點云的高密度性和強穿透性,本研究僅選取林分結(jié)構(gòu)單一的人工林為研究對象,未對林間情況更復雜的天然林進行郁閉度估測;由于魚眼鏡頭使用條件限制,所選取的樣本點分布不夠均勻;另外,由于UAV-LiDAR 技術(shù)處于發(fā)展階段,歷史數(shù)據(jù)缺乏,現(xiàn)階段對比分析不同時期的森林監(jiān)測數(shù)據(jù)難度較大。在未來的針對天然林郁閉度估測中,有待選取大量具有代表性的天然林樣本,開發(fā)更精確、更簡潔的算法增強模型的魯棒性,如采用多元線性回歸(multiple linear regression,MLR)、偏最小二乘回歸法和隨機森林(random forest,RF)等多種方法構(gòu)建反演模型,討論評價不同方法對估測精度的影響;另外,聯(lián)合GF、Landsat-TM、ETM+數(shù)據(jù),對激光雷達計算的森林參數(shù)進行準確性和可靠性驗證,可為人工林精細化管理和森林資源監(jiān)測提供技術(shù)支持。
1) 本研究構(gòu)建的基于無人機激光雷達數(shù)據(jù)的人工林郁閉度估測模型具有較高精度,且人工闊葉林郁閉度估測模型的判定系數(shù)、調(diào)整判定系數(shù)和均方根誤差等精度指標均優(yōu)于人工針葉林。
2) 引入多類特征變量估測林分郁閉度效果較好,采用逐步回歸法建立統(tǒng)計模型反演的郁閉度與魚眼照片提取的郁閉度之間有良好的關(guān)系。影響針葉林和闊葉林林分郁閉度反演的特征變量不同,冠層特征變量在針葉林郁閉度反演中所占權(quán)重更大。
3) 利用回歸預測模型和空間插值法得到整個樣地的郁閉度,檢驗點郁閉度與實測郁閉度顯示出較高相關(guān)性(r=0.859),表明利用該預測模型和反距離權(quán)重插值法反演制圖預測林分郁閉度是可行的且優(yōu)勢明顯。本研究對東北地區(qū)常見的針葉林興安落葉松、黑皮油松、樟子松,闊葉林白樺、蒙古櫟進行林分郁閉度預測對比研究,取得較好預測效果,表明UAVLiDAR 點云在估測東北地區(qū)針葉林和闊葉林林分郁閉度上具有較好潛力。