基于無人機激光雷達的銀杏人工林有效葉面積指數(shù)估測*

吳項乾 曹 林 申 鑫 汪貴斌 曹福亮

(南京林業(yè)大學 南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心 南京 210037)

人工林是森林生態(tài)系統(tǒng)不可或缺的一部分，不僅可為社會經(jīng)濟發(fā)展提供優(yōu)質(zhì)豐富的木材資源，而且在維持生態(tài)系統(tǒng)平衡以及減緩全球變暖等環(huán)境保護方面也發(fā)揮著重要作用(Brockerhoffetal.， 2008)。我國是世界上人工林面積最大的國家，約69.3萬km2，超過世界人工林總面積的1/4(FAO， 2010)。銀杏(Ginkgobiloba)在我國種植面積廣、時間久，40年及以上銀杏株數(shù)達28萬株，僅果實年產(chǎn)量就達1.6萬t(曹福亮， 2007)。近年來，隨著銀杏人工林種植面積的迅速增長，不僅保障了國家和社會對林木產(chǎn)品的市場需求，而且在森林可持續(xù)發(fā)展和生態(tài)環(huán)境重建等方面也具有重要價值。

森林冠層是植物進行蒸騰、呼吸和光合作用等生理過程的主要場所，葉面積指數(shù)(leaf area index，LAI)是描述冠層結(jié)構(gòu)及葉片數(shù)量和分布情況的有效定量指標(Bequetetal.， 2011； Ryuetal.， 2012； Gonsamoetal.， 2014)，實時、快速、精確估測人工林葉面積指數(shù)并掌握其連續(xù)空間分布，是實現(xiàn)人工林可持續(xù)經(jīng)營的前提，同時也是人工林生態(tài)循環(huán)的研究基礎(chǔ)(Wulderetal.， 2010； Thomasetal.， 2011)。精確估測銀杏葉面積指數(shù)并充分掌握其冠層結(jié)構(gòu)特點，是預(yù)測銀杏人工林地上生物量的重要參考依據(jù)(Runningetal.， 1988)，同時也有助于掌握銀杏人工林的競爭與生長、評價其立地質(zhì)量及理解人工林生態(tài)系統(tǒng)的功能和生產(chǎn)力(Turneretal.， 2000； 張佳華等， 2002)。

傳統(tǒng)以地面實測為主的葉面積指數(shù)測量主要通過直接法和儀器測量，可以獲得較為精確的葉面積指數(shù)； 但地面實測存在效率低、費時費力等問題，且往往只能獲得樣地尺度的葉面積指數(shù)，較難獲得林場尺度及更大范圍的葉面積指數(shù)連續(xù)分布信息； 此外，某些直接法(如鮮重打孔法、描形稱重法)可能會在一定程度上破壞森林冠層結(jié)構(gòu)(程武學等， 2010)。激光雷達(light detection and ranging，LiDAR)作為一種主動遙感技術(shù)，其原理是通過發(fā)射激光脈沖到物體表面，然后接收并分析其返回信號，利用LiDAR可以獲得高精度的地表及地物的垂直結(jié)構(gòu)信息，如地形和植被等的理論垂直精度可達到15～30 cm(李增元等， 2016)。相較于傳統(tǒng)的光學遙感技術(shù)，LiDAR數(shù)據(jù)獲取不易受天氣及光照條件影響，在獲取森林三維結(jié)構(gòu)方面具有獨特優(yōu)勢(李德仁等， 2012； 曹林等， 2014； Wulderetal.， 2010； Thomasetal.， 2011)。近年來，隨著無人機(unmanned aerial vehicle，UAV)遙感技術(shù)的迅速發(fā)展，無人機激光雷達在估測林分結(jié)構(gòu)參數(shù)方面的研究也逐漸增多(Linetal.， 2011； Wallaceetal.， 2012)。Nagai等(2009)將低頻激光掃描儀(SICK LMS-291)和單反相機(Canon EOS 10D)搭載在無人直升機上用于提取冠層數(shù)字表面模型(DSM)和紋理特征，制圖平均誤差為3～10 cm； Jaakkola等(2010)將激光掃描儀(Sick LMS151)、CCD相機(AVT Pike F-421C)、光譜儀(Specim V10H)和熱紅外相機等搭載在輕巧型多旋翼無人機激光雷達系統(tǒng)上(FGI Sensei)用于估測芬蘭針葉林樹高，其標準差為13 cm，R2=0.92； Wallace等(2014； 2016)研制出一款低成本激光雷達系統(tǒng)(TerraLuma)，該系統(tǒng)合成了低頻掃描儀(Ibeo LUX)和高清錄像機，用于澳大利亞東南部桉樹(Eucalyptus)單木探測和樹木高度估測，單木探測準確率達98%，樹木高度估測的R2=0.69～0.84，RMSE=0.92～1.3 m。

目前，基于LiDAR數(shù)據(jù)估測葉面積指數(shù)的方法可歸納為2類。一類是“統(tǒng)計模型法”，即利用LiDAR數(shù)據(jù)提取的冠層結(jié)構(gòu)特征變量與地面實測LAI構(gòu)建統(tǒng)計模型，再采用該統(tǒng)計模型進行估測。如Lim等(2003)利用LiDAR點云計算的特征變量與49塊地面樣地的有效葉面積指數(shù)【effective leaf area index，eLAI，即未考慮到葉片的聚集度(clumping index)】實測數(shù)據(jù)進行回歸分析，估測了加拿大北部闊葉林的eLAI； Roberts等(2005)使用LiDAR數(shù)據(jù)結(jié)合多元非線性模型方法，預(yù)測了火炬松(Pinustaeda)人工林單木尺度的eLAI(R2=0.61，RMSE=0.23 m2·m-2)。另一類是“孔隙度模型法”，即通過計算激光點云的穿透率，根據(jù)比爾-朗博定律(Beer-Lambert law)計算葉面積指數(shù)。如Morsdorf等(2006)使用LiDAR數(shù)據(jù)結(jié)合Beer-Lambert定律估測瑞士石松(Lycopodiumjaponicum)的eLAI，通過機載小光斑離散點云數(shù)據(jù)冠層回波數(shù)與地面回波數(shù)的比值來遙測eLAI，同時對樣地尺度圓形(樣地半徑2～25 m)進行敏感性測試，結(jié)果發(fā)現(xiàn)在樣地半徑為15 m時，LiDAR估測eLAI的效果最佳，其中R2達到0.69，RMSE僅為0.01 m2·m-2； Sasaki等(2008)也利用機載LiDAR數(shù)據(jù)點云和近紅外高光譜數(shù)據(jù)估測了位于日本大阪府秀塔市溫帶常綠闊葉林的有效葉面積指數(shù)和冠層孔隙度，通過點云的冠層孔隙度計算得到eLAI，結(jié)果發(fā)現(xiàn)LiDAR估測eLAI的精度(R2=0.76～0.83，RMSE=0.23～0.28 m2·m-2)優(yōu)于光學數(shù)據(jù)估測值； Zhao等(2009)對比LiDAR點云數(shù)據(jù)、光學影像數(shù)據(jù)以及2種數(shù)據(jù)融合3種方式來估測eLAI，得出LiDAR估測效果更好，更能滿足樣地到區(qū)域不同尺度eLAI制圖的要求。

綜觀以往研究，大多為單獨利用孔隙度模型法或統(tǒng)計模型法估測有效葉面積指數(shù)，且多針對歐美北寒帶、溫帶的天然林，而關(guān)于人工林方面類似的研究較少。鑒于此，本研究以我國具有典型性和代表性的成片實生銀杏林為對象，利用無人機激光雷達獲取的高密度LiDAR點云數(shù)據(jù)，結(jié)合45塊地面實測樣地數(shù)據(jù)，使用孔隙度模型法和統(tǒng)計模型法對研究區(qū)內(nèi)樣地尺度的有效葉面積指數(shù)進行估測并驗證精度。具體內(nèi)容包括： 1) 分別使用孔隙度模型法和統(tǒng)計模型法估測eLAI并評價精度； 2) 在通過LiDAR高度特征變量預(yù)測eLAI的基礎(chǔ)上，逐步引入冠層密度變量、冠層容積比特征和強度特征變量，并評估幾組無人機LiDAR高密度點云特征變量對eLAI預(yù)測的效果。

1 研究區(qū)概況與研究方法

選取江蘇省邳州市鐵富鎮(zhèn)銀杏人工林為研究區(qū)(34°33′49″—34°34′23″N，118°05′1″—118°06′06″E)，該區(qū)氣候?qū)侔霛駶櫆貛Ъ撅L氣候，年降雨量約903 mm，總面積約5 400 hm2(本研究從中選取2 197 hm2銀杏人工林)，林下有部分大豆(Glycinemax)等林農(nóng)復(fù)合經(jīng)營的低矮植被作物。基于八旋翼無人機平臺(UAV)獲取的高密度激光雷達點云數(shù)據(jù)和樣地實測有效葉面積指數(shù)(eLAI)，結(jié)合孔隙度模型法、統(tǒng)計模型法估測樣地尺度eLAI。步驟如下： 首先，對LiDAR原始點云數(shù)據(jù)進行去噪、分類處理，通過插值方法得到數(shù)字高程模型(digital elevation model，DEM)，利用DEM進行高度歸一化以獲得歸一化點云LiDAR數(shù)據(jù)； 然后，將樣地實測eLAI與提取的多組LiDAR特征變量結(jié)合多元逐步回歸方法構(gòu)建eLAI估測模型(統(tǒng)計模型法)，同時，根據(jù)比爾-朗博定律(Beer-Lambert law)結(jié)合LiDAR點云提取的孔隙度(P)估測eLAI(孔隙度模型法)； 最后，利用樣地實測有效葉面積指數(shù)對模型的精度和實用性進行評價。研究技術(shù)路線如圖1所示。

圖1 研究技術(shù)路線Fig.1 Technological route of this research

1.1 樣地選取

基于森林資源二類調(diào)查數(shù)據(jù)、林場種植歷史記錄數(shù)據(jù)和高分辨率影像數(shù)據(jù)，在研究區(qū)銀杏核心分布范圍選取5塊1 km×1 km的正方形無人機采樣區(qū)，按照典型選樣方法在5塊樣地內(nèi)布設(shè)圓形小樣地各9塊(半徑為15 m)，樣地中心坐標由Trimble GeoXH6000 Handhelds手持GPS獲取，并通過JSCROS廣域差分系統(tǒng)進行位置校正以獲取亞米級精度的坐標信息。研究區(qū)及樣地分布如圖2所示。

1.2 地面數(shù)據(jù)采集

地面調(diào)查時間為2016年10月，在研究區(qū)設(shè)置45塊圓形樣地(半徑為15 m)。樣地調(diào)查過程中，對胸徑超過5 cm的單木逐一測定其胸徑、樹高、枝下高和冠幅，并將測量所得單木信息統(tǒng)計得到樣地尺度林分統(tǒng)計特征，即平均胸徑(cm)、Lorey’s平均樹高(m)、樣地株數(shù)、樣地密度(hm-2)和有效葉面積指數(shù)(m2·hm-2)。根據(jù)提取的激光雷達Cover特征量(即首次回波中高于1 m的激光返回點占所有返回點的比例)大小將樣地劃分為3組(第1組Cover： 0.08～0.20； 第2組Cover： 0.21～0.31； 第3組Cover： 0.32～0.83)。通過對比h95(森林冠層首次回波高度垂直分布的95分位數(shù))、hcv(歸一化點云高度的變異系數(shù))和Cover這3個激光雷達特征變量的分組效果發(fā)現(xiàn)，Cover分組對葉面積指數(shù)的估測精度最高，可以用來評價樣地的樹冠覆蓋度以表示樣地內(nèi)單木生長程度。比較分組估測模型與總體估測模型的精度發(fā)現(xiàn)，總體模型估測精度高于分組模型。研究區(qū)樣地分組及參數(shù)匯總?cè)绫?所示。

圖2 研究區(qū)正射影像及樣地空間分布Fig.2 The orthophoto of study area and the spatial distribution of the plots

表1 研究區(qū)樣地林分參數(shù)調(diào)查信息匯總Tab.1 The summary of plot-level stand characteristics

有效葉面積指數(shù)(eLAI)通過LAI-2200測量獲得，考慮天頂角A的最外圈數(shù)據(jù)可能存在地物遮擋以及藍光多次散射效應(yīng)的影響，使用1～4內(nèi)圈數(shù)據(jù)計算有效葉面積指數(shù)。具體操作方法如下： 在樣地內(nèi)選取2條平行的長度30 m測線(方向垂直于太陽光照)，2條測線中點分別距樣地中心7.5 m。將LAI-2200置于距離地面1 m高的位置進行測量，以防止低矮植被的影響。測量前在空曠的林窗下匹配2支探桿的時間和測量參數(shù)，并保證視角蓋均為90°且方向一致； 匹配后，其中一支探測桿置于林窗下每隔10 s自動采樣1次(A)，另一支在樣地內(nèi)沿著測線每隔4 m采樣1次(B)； 測量結(jié)束后，A和B通過時間匹配，聯(lián)合計算eLAI。每個樣區(qū)內(nèi)16個測點的eLAI平均值作為該樣區(qū)的eLAI。有效葉面積指數(shù)地面實測方法如圖3所示。

1.3 LiDAR點云獲取及預(yù)處理

LiDAR數(shù)據(jù)由八旋翼無人機平臺搭載的Velodyne激光掃描儀(VLP-16)于2016年10月13—17日獲取。無人機平臺飛行高度60 m，飛行速度4.8 m·s-1，條帶之間間隔50 m，LiDAR掃描儀激光波長903 nm，掃描角設(shè)定±30°，發(fā)散度3 mrad，地面光斑直徑18 cm，脈沖發(fā)射頻率21.7 kHz，掃描頻率16 s-1，最終獲取的LiDAR點云平均點密度約159 pts·m-2。對于獲取的原始數(shù)據(jù)， 首先根據(jù)閾值法對噪聲點進行識別并刪除，同時對點云地面點進行識別分類； 然后通利用反距離加權(quán)法(inverse distance weighted，IDW)結(jié)合分類的地面點插值以生成DEM； 最后利用DEM對點云進行高度(z)歸一化以獲得高度歸一化點云數(shù)據(jù)。

1.4 LiDAR特征變量提取

使用DEM將離散點云歸一化后，計算樣地尺度的激光雷達特征變量。本研究共計算6類(4組)特征變量，特征變量及描述如表2所示。

圖3 有效葉面積指數(shù)地面實測方法示意Fig.3 The sketch of eLAI measured on ground

表2 LiDAR特征變量及其描述Tab.2 The description of LiDAR metrics

枝葉剖面(foliage profile，F(xiàn)P)是樹木枝葉在垂直空間結(jié)構(gòu)上的分布狀況，描述了冠層內(nèi)葉面積指數(shù)隨著高度增長的變化情況(Lovelletal.， 2003)，F(xiàn)P值為樹木冠層指定高度上單位體積內(nèi)的枝葉面積之和(Riaoetal.， 2003)。FP與LAI 在垂直高度上的分布密切相關(guān)(Radtkeetal.， 2011)，其表達式為：

(1)

式中：L(z)為在高度z處的累積有效葉面積指數(shù)；z1、z2為冠層高度。

采用Weibull函數(shù)擬合FP，本研究Weibull參數(shù)不是根據(jù)胸徑算出的表示樣地徑階的參數(shù)(模型因變量)，而是利用LiDAR垂直枝葉剖面的Weibull擬合計算的特征變量(模型自變量)，公式如下：

(2)

式中：Hmax為最大樹高；α為尺度參數(shù)；β為形狀參數(shù)(Coopsetal.， 2007)。

Weibull函數(shù)中，α和β通常通過最大似然估計法(maximum likelihood estimation，MLE)解算求得。先前研究表明，MLE在預(yù)測森林Weibull分布函數(shù)時擬合效果較好(Gorgosoetal.， 2007)，其計算公式如下：

(3)

(4)

式中：n為樣本個數(shù)；xi為每株單木的冠層容積。

1.5 孔隙度模型法估測有效葉面積指數(shù)

根據(jù)光傳輸特性，當光通過某種介質(zhì)后，能量會有所減弱。如果將森林冠層當成一種均勻介質(zhì)，那么影響光穿透后衰減的重要原因是冠層“孔隙度”，可通過Beer-Lambert定律算出。該定律假設(shè)光進入冠層前的初始入射能量為I0，通過距離為m的均勻介質(zhì)后，衰減后的能量為I，公式如下：

I=I0exp(-m)。

(5)

植物冠層消光系數(shù)直接決定其對輻射的截獲量，因為穿過的介質(zhì)為均勻介質(zhì)。若設(shè)消光系數(shù)為a，介質(zhì)濃度為C，那么有：

(6)

將式(4)、(5)合并，得：

I=I0exp(-aCm)；

(7)

(8)

式中：I0/I表示光的穿透率，設(shè)為P。

消光系數(shù)G(θ)是與葉傾角密切相關(guān)的系數(shù)，P可以用下式表示：

P(θ)=e-G(θ)L/cosθ。

(9)

根據(jù)激光雷達計算葉面積指數(shù)：

(10)

式中：θ為激光雷達掃描角。

本研究區(qū)內(nèi)為銀杏純林樣地，利用19塊單木定位樣地(使用Postex林木定位儀確定每株銀杏單木位置)，通過最近鄰體分析法(nearest neighbor analysis，NNA)計算研究區(qū)內(nèi)樹木水平分布情況，結(jié)果表明銀杏單木在水平方向上均呈隨機分布(胡艷波等， 2014； 趙中華等， 2016)。因銀杏為喜光闊葉樹種，為了最大化提高光能利用率進行光合作用，冠層中的葉片會向各方向生長，整體葉傾角在空間上近似為球形分布。假設(shè)激光垂直入射冠層，因此可將消光系數(shù)設(shè)為0.5(Weissetal.， 2004； Morsdorfetal.， 2006)。P(θ)可通過下式算出：

(11)

式中：n為LiDAR發(fā)射的總光點數(shù)；nground為LiDAR穿透森林冠層后到達地面的光點數(shù)。

基于以上公式，可通過孔隙度模型法算出有效葉面積指數(shù)，使用地面實測有效葉面積指數(shù)對孔隙度模型法估測的有效葉面積指數(shù)進行精度評價。

1.6 統(tǒng)計模型法估測有效葉面積指數(shù)

通過逐步回歸法構(gòu)建的線性方程預(yù)測有效葉面積指數(shù)，模型形式如下：

y=c0+c1q1+c2q2+…+cnqn+ε。

(12)

式中：y為因變量；q為自變量；c為常數(shù)；ε為誤差。

本研究中，地面實測有效葉面積指數(shù)為因變量，LiDAR提取的多個特征變量組為自變量，將自變量逐個引入模型，每引入1個自變量都進行F檢驗(P<0.05)，當原先引入的自變量由于后面自變量的引入變得不再顯著時，則將其刪除，直到將所有自變量篩選完為止，最終篩選出與有效葉面積指數(shù)相關(guān)性最高的3個自變量(自變量的數(shù)量并未增加)進行建模(Nssetetal.， 2008； Breidenbachetal.， 2010)。采用決定系數(shù)(R2)評價模型擬合效果，采用均方根誤差(RMSE)、相對均方根誤差(rRMSE)和赤池信息量準則(AIC)評價回歸模型精度(Akaike，1974； Getisetal.， 2004)。計算公式如下：

(13)

(14)

(15)

AIC=-lgL+2k。

(16)

在模型參數(shù)已確定的情況下，對45塊地面實測樣地全部進行交叉驗證(cross-validation)，首先隨機選取任意1塊樣地的eLAI作為驗證數(shù)據(jù)，將其余44塊樣地的eLAI進行建模，然后利用擬合的模型與隨機選取1塊樣地的eLAI進行檢驗，如此重復(fù)進行直到所有樣地全部驗證完畢。

2 結(jié)果與分析

表3為選取不同LiDAR特征變量組合參與建模估測eLAI的模型形式及精度評價。由表可見，僅使用高度特征變量建模預(yù)測eLAI時精度最低(R2=0.44，rRMSE=48.3%，AIC=202.5)，逐步引入其他特征變量后，預(yù)測精度明顯提升。使用所有LiDAR特征變量建立模型估測eLAI時精度最高(R2=0.78，rRMSE=22.1%，AIC=158.3)，高于使用高度和密度變量建立模型估測eLAI的精度(R2=0.67，rRMSE=43.2%，AIC=172.1)，再其次是使用高度和冠層容積變量建立模型變量預(yù)測eLAI的精度(R2=0.65，rRMSE=43.8%，AIC=175.2)。

表4為不同Cover分組的有效葉面積指數(shù)估測模型比較及交叉驗證精度評價。由表可見，通過Cover將45塊樣地分為3組后，每組模型交叉驗證精度相比于總體樣地模型都有所提升。

表3 不同特征變量組合作為預(yù)測變量的有效葉面積指數(shù)估測模型比較及精度評價Tab.3 eLAI model by different metrics and their accuracy assessments

表4 不同Cover分組的有效葉面積指數(shù)估測模型比較及交叉驗證精度評價Tab.4 eLAI models by different Cover and cross-validation results

圖4 地面實測eLAI與統(tǒng)計模型法預(yù)測eLAI交叉驗證結(jié)果Fig.4 The cross-validation results of the filed-measured and statistical model predicted eLAIa為基于高度特征變量建模； b為基于高度+密度特征變量建模； c為基于高度+冠層容積特征變量建模； d為基于所有提取的特征變量(含強度特征變量)建模。a is modelled by height percentile； b is modelled by height percentile and density metrics； c is modelled by height percentile and canopy volume metrics； d is modelled by all metrics(contain intensity metrics).

統(tǒng)計模型法預(yù)測有效葉面積指數(shù)模型交叉驗證結(jié)果如圖4所示。僅通過LiDAR高度特征變量參與建模的精度最低，R2為0.38(rRMSE=54%)(圖4a)； 使用高度和密度特征變量組合參與建模的精度為R2=0.64(rRMSE=26%)(圖4b)； 通過高度和冠層容積特征變量組合與地面實測eLAI建立模型的精度為R2=0.61(rRMSE=28%)(圖4c)； 當使用所有LiDAR特征變量(包括強度特征變量)建模時精度最高R2為0.74(rRMSE=23%)(圖4d)。加入其他類型特征變量后，比僅通過高度特征變量建立模型的預(yù)測效果都有不同程度提高，其中，基于所有提取的特征變量建模預(yù)測效果最佳，通過高度和密度特征變量組合估測eLAI效果次之，最后是通過高度和冠層容積特征變量組合估測eLAI的精度。

根據(jù)Cover分組后，不同組統(tǒng)計模型法預(yù)測有效葉面積指數(shù)估測模型交叉驗證結(jié)果如圖5所示。分組后模型估測精度皆優(yōu)于不分組估測精度(R2=0.74，rRMSE=23.0%)(圖5d)； 在分組估測模型中，第2組估測有效葉面積指數(shù)精度最高(R2=0.88，rRMSE=12.6%)(圖5b)，其次是第1組(R2=0.86，rRMSE=15.6%)(圖5a)，第3組估測有效葉面積指數(shù)精度最低(R2=0.79，rRMSE=21.9%)(圖5c)。

孔隙度模型法預(yù)測eLAI精度如圖6所示。相較于4種點云特征變量統(tǒng)計模型法估測精度，孔隙度模型法估測精度(R2=0.71，rRMSE=32.0%)高于僅通過高度百分數(shù)建模估測、高度和密度特征變量組合建模估測以及高度和冠層容積特征變量組合估測eLAI的精度，但比通過所有特征變量(包括強度特征變量)建模預(yù)測精度低。

根據(jù)激光雷達算出的銀杏人工林Cover特征量，將45塊樣地分為3組，選取3組中具有典型代表性樣地如圖7所示。對比3組樣地正射影像，第1組樣地密度最高，第2組樣地密度最低，第3組樣地密度較第2組有所增高(a1—c1)； 對比3組樣地球面鏡影像，第1組樣地天空所占比例最大，第2組樣地次之，第3組樣地天空所占比例最小(a2—c2)； 觀察3組樣地激光雷達點云及點云剖面圖，第1組樣地平均樹高最低(約8.14 m)但激光雷達點密度最高，第2組樣地平均樹高比第1組樣地更高(約11.52 m)但激光雷達點密度降低，第3組樣地平均樹高最高(約12.76 m)但激光雷達點密度比第2組高(a3—c3，a4—c4)。

圖5 不同Cover組地面實測eLAI與統(tǒng)計模型法預(yù)測eLAI交叉驗證結(jié)果Fig.5 The cross-validation results of the filed-measured and statistical model predicted eLAI indifferent Covera為第1組結(jié)果； b為第2組結(jié)果； c為第3組結(jié)果； d為所有樣地結(jié)果。a is the result of Group 1； b is the result of Group 2； c is the result of Group 3； d is the result of all plots.

圖6 孔隙度模型法估測eLAI與地面實測eLAI擬合精度Fig.6 eLAI estimated by gap-fraction model

3組樣地冠層垂直結(jié)構(gòu)分布剖面【計算方法詳見Coops等(2007)】如圖8所示。由3組典型樣地冠層高度分布(canopy height distribution，CHD)可知，第1組樣地激光雷達點密度最高處在3～5 m之間，第2組樣地激光雷達點密度最高處在7.2～10.5 m之間，第3組樣地激光雷達點密度最高處在3～11 m之間，3組樣地平均樹高依次增高，冠幅增大(a1—c1)； 由3組樣地枝葉剖面(FP)隨冠層高度變化情況可知，第1組樣地枝葉剖面峰值最大在約3.5 m樹高處，第2組樣地枝葉剖面峰值最大在約8.5 m樹高處，第3組樣地枝葉剖面峰值最大在約9.0 m樹高處，3組樣地最大峰值在冠層中的高度分布依次增高，表明冠幅最大的枝葉的高度變高(a2—c2)，且CHD和FP分別與Weibull函數(shù)曲線擬合形式較好； 由3組樣地的有效葉面積指數(shù)累積曲線可知，最高處的eLAI為該樣地的eLAI，第1組樣地eLAI為0.48 m2·m-2，第2組樣地eLAI為1.8 m2·m-2，第3組樣地eLAI為2.3 m2·m-2。

3 討論

圖7 3組典型樣地正射影像(a1—c1)、球面鏡影像(a2—c2)、激光雷達點云(a3—c3)及激光雷達點云剖面(含50、75、95高度分位數(shù)位置)(a4—c4)Fig.7 The orthophotos(a1—c1), spherical HemiVeiw photos(a2—c2), point clouds (a3—c3)and profile(contain 50th, 75th, 95th height percentiles)(a4—c4)of three typical plots

本研究基于無人機激光雷達數(shù)據(jù)并結(jié)合45塊地面銀杏人工林樣地實測數(shù)據(jù)進行eLAI估測，結(jié)果發(fā)現(xiàn)，使用統(tǒng)計模型法構(gòu)建eLAI估測模型時，引入多種變量(密度、冠層容積和強度特征變量)后，估測精度高于僅使用高度特征變量建立模型的精度。Nsset(2002)利用機載激光雷達反演挪威不同林齡森林結(jié)構(gòu)參數(shù)，通過提取激光雷達點云冠層高度和冠層密度變量，得出高度和密度變量組合的估測精度優(yōu)于僅通過高度或密度變量估測精度； Cao等(2016)通過激光雷達估測亞熱帶森林地面生物量時，利用激光雷達點云特征變量和強度相關(guān)信息聯(lián)合建模精度高于僅通過點云特征變量建模的精度； Morsdorf等(2006)通過建立離散小光斑激光雷達不同尺度(圓形半徑2～25 m)的冠層回波數(shù)和地面回波數(shù)的比值與基于半球影像估測的LAI進行回歸分析，建立LAI估測經(jīng)驗方程，發(fā)現(xiàn)光斑半徑為15 m時效果最好，決定系數(shù)為0.69，其孔隙度模型法估測eLAI精度與本研究相近； Peduzzi等(2012)利用多回波LiDAR數(shù)據(jù)計算樣地尺度的LAI以及不同高度處的覆蓋度等變量與LAI地面實測值進行回歸分析，并用于卡羅萊納州和弗尼吉亞州北部松樹林LAI估測，決定系數(shù)最高達0.77，RMSE為0.42，與本研究統(tǒng)計模型法估測eLAI精度相近。利用統(tǒng)計模型法估測eLAI時，引入多種特征變量后的估測精度高于僅通過高度特征變量估測精度，這是因為在引入其他點云特征變量后，新的特征變量組合對冠層結(jié)構(gòu)刻畫更為全面，從而提升了有效葉面積指數(shù)估測精度。本研究選取強度特征變量估測有效葉面積指數(shù)，取得了較高精度(R2=0.74，rRMSE=23%)，與前人研究結(jié)果相似。Mariano等(2010)利用機載LiDAR的高度和強度特征變量估測地中海地區(qū)的森林生物量，結(jié)果發(fā)現(xiàn)使用高度和強度變量組合建模的精度明顯優(yōu)于僅通過高度特征變量建模的精度，同時，預(yù)測模型也選取了較高的強度特征變量作為預(yù)測因子(I95)。本研究中由于無人機飛行高度恒定，且研究區(qū)處于平原地區(qū)； 同時，樣地大多處于無人機激光雷達掃描的下方，避免了能量因為傳輸距離增大(掃描角增大)的損耗，因此使得未校正的強度變量適用于本研究有效葉面積指數(shù)的估算(Kiddetal.， 2016； Brunaetal.， 2019)。

圖8 3個典型樣地冠層高度分布(CHD)(a1—c1)、枝葉剖面(FP)(a2—c2)與Weibull擬合曲線及有效葉面積指數(shù)(eLAI)累積曲線(a3—c3)Fig.8 Canopy height distribution(a1—c1)and foliage profile(a2—c2)fitted with Weibull function, and eLAI accumulate(a3—c3)

對比圖7中的 3個組別，第1組樹高較低，但種植密度較大，激光雷達點云多集中在冠層中上部，結(jié)合銀杏人工林的生長情況，第1組為幼齡林； 相對于第1組，第2組密度有所降低，但是樹木高度增大，此組類型為隨著樹木生長發(fā)育進行一系列撫育間伐，林分覆蓋度有所增加，因此激光雷達點云也多集中在冠層中上部，穿透到地面的點云數(shù)量較少； 第3組中，樹高較前2組有所增大，密度與第2組相比較為相近，但覆蓋度增加，根據(jù)銀杏生長規(guī)律，銀杏屬于強喜光樹種，隨著樹齡增加，對光照的要求也愈加迫切(曹福亮， 2007)，其中上部冠層會因為趨光性而逐漸向有陽光照射區(qū)域生長，冠幅變大，只有極少部分激光雷達點云可以穿透樹冠到達地面，此組銀杏為成熟林，相較于前2個生長階段，激光雷達點云分布更偏于上層。綜合3個生長階段激光點云分布情況，與模型選取的75%、95%高度分位數(shù)來估測eLAI的情況一致。

從冠層枝葉以及非光合作用組織在垂直高度上的點云分布情況(圖8)可知，3組樣地樹木依次增加，所對應(yīng)的生長時期也不同。第1組樣地內(nèi)平均樹高較低，栽植密度大，中下部點云密度較多，由于是幼齡林，銀杏的粗壯枝葉集中在中層。第2組樣地內(nèi)平均樹高較第1組增加，密度降低，隨著銀杏生長，進行撫育間伐，主枝變大并朝上生長，因此冠層中上部點云分布較多。第3組樣地平均樹高比第2組樣地增加，但冠層中高密度點云向下偏移，冠層中部點云密度最大，這是由于銀杏林逐漸成長為成熟林，主枝不再隨著樹高增長而升高，因銀杏趨光性而向有光照區(qū)域生長，主枝變粗大，冠層中部枝葉覆蓋度增加，所以在樹高6～10 m范圍內(nèi)，激光雷達點云數(shù)分布較多。圖8b1—b3為枝葉剖面(FP)，描述了在垂直高度上冠層枝葉面積的分布情況(Juppetal.， 2009)。隨著樹齡增長，銀杏林分不斷成熟，其枝葉剖面分布更加均勻且峰值向上偏移； CHD和FP與Weibull函數(shù)擬合較好(圖8a1—a3、b1—b3)，基本可以得到分布的形狀和范圍。圖8c1—c3反映了有效葉面積指數(shù)隨高度的變化情況，3組樣地中，冠層中部eLAI增長速度較快，底部和頂部增長速度較慢，這是由于銀杏中部冠層生長較為茂盛，集中了較多的粗大枝。本研究利用無人機激光雷達數(shù)據(jù)對研究區(qū)內(nèi)樣地尺度的有效葉面積指數(shù)進行估測和方法比較及精度驗證，未能進行尺度變換實現(xiàn)整個無人機數(shù)據(jù)獲取區(qū)域(包含5個“無人機采樣區(qū)”)的葉面積指數(shù)估測，在今后的研究中，可以嘗試進行有效葉面積指數(shù)空間分布外推，實現(xiàn)整個研究區(qū)的有效葉面積指數(shù)空間分布估測。

4 結(jié)論

本研究利用多旋翼無人機激光雷達系統(tǒng)獲取的點云數(shù)據(jù)，結(jié)合45塊地面實測樣地數(shù)據(jù)，使用孔隙度模型法(通過計算點云的冠層穿透率，根據(jù)Beer-Lambert 定律計算有效葉面積指數(shù))和統(tǒng)計模型法(首先通過地面實測的有效葉面積指數(shù)和所提取的LiDAR特征變量建模，然后借助擬合的模型估測有效葉面積指數(shù))對我國典型平原銀杏人工林進行樣地尺度的有效葉面積指數(shù)估測。結(jié)果表明： 統(tǒng)計模型法估測有效葉面積指數(shù)時，僅通過高度特征變量估測精度為R2=0.38(rRMSE=54%)，分別引入多種特征變量，即冠層密度變量、冠層容積比特征和強度特征后，估測葉面積指數(shù)精度有所提高，精度分別為R2=0.64(rRMSE=26%)、R2=0.61(rRMSE=28%)、R2=0.74(rRMSE=23%)； 利用Cover分組后建模結(jié)果表明，分組建模的精度優(yōu)于不分組建模的精度； 孔隙度模型法估測有效葉面積指數(shù)的精度為R2=0.71(rRMSE=32.0%)。