基于無(wú)人機(jī)激光雷達(dá)的銀杏人工林有效葉面積指數(shù)估測(cè)*

吳項(xiàng)乾 曹 林 申 鑫 汪貴斌 曹福亮

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

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

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

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

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

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

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

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

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

1.1 樣地選取

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

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

地面調(diào)查時(shí)間為2016年10月，在研究區(qū)設(shè)置45塊圓形樣地(半徑為15 m)。樣地調(diào)查過(guò)程中，對(duì)胸徑超過(guò)5 cm的單木逐一測(cè)定其胸徑、樹(shù)高、枝下高和冠幅，并將測(cè)量所得單木信息統(tǒng)計(jì)得到樣地尺度林分統(tǒng)計(jì)特征，即平均胸徑(cm)、Lorey’s平均樹(shù)高(m)、樣地株數(shù)、樣地密度(hm-2)和有效葉面積指數(shù)(m2·hm-2)。根據(jù)提取的激光雷達(dá)Cover特征量(即首次回波中高于1 m的激光返回點(diǎn)占所有返回點(diǎn)的比例)大小將樣地劃分為3組(第1組Cover： 0.08～0.20； 第2組Cover： 0.21～0.31； 第3組Cover： 0.32～0.83)。通過(guò)對(duì)比h95(森林冠層首次回波高度垂直分布的95分位數(shù))、hcv(歸一化點(diǎn)云高度的變異系數(shù))和Cover這3個(gè)激光雷達(dá)特征變量的分組效果發(fā)現(xiàn)，Cover分組對(duì)葉面積指數(shù)的估測(cè)精度最高，可以用來(lái)評(píng)價(jià)樣地的樹(shù)冠覆蓋度以表示樣地內(nèi)單木生長(zhǎng)程度。比較分組估測(cè)模型與總體估測(cè)模型的精度發(fā)現(xiàn)，總體模型估測(cè)精度高于分組模型。研究區(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)通過(guò)LAI-2200測(cè)量獲得，考慮天頂角A的最外圈數(shù)據(jù)可能存在地物遮擋以及藍(lán)光多次散射效應(yīng)的影響，使用1～4內(nèi)圈數(shù)據(jù)計(jì)算有效葉面積指數(shù)。具體操作方法如下： 在樣地內(nèi)選取2條平行的長(zhǎng)度30 m測(cè)線(方向垂直于太陽(yáng)光照)，2條測(cè)線中點(diǎn)分別距樣地中心7.5 m。將LAI-2200置于距離地面1 m高的位置進(jìn)行測(cè)量，以防止低矮植被的影響。測(cè)量前在空曠的林窗下匹配2支探桿的時(shí)間和測(cè)量參數(shù)，并保證視角蓋均為90°且方向一致； 匹配后，其中一支探測(cè)桿置于林窗下每隔10 s自動(dòng)采樣1次(A)，另一支在樣地內(nèi)沿著測(cè)線每隔4 m采樣1次(B)； 測(cè)量結(jié)束后，A和B通過(guò)時(shí)間匹配，聯(lián)合計(jì)算eLAI。每個(gè)樣區(qū)內(nèi)16個(gè)測(cè)點(diǎn)的eLAI平均值作為該樣區(qū)的eLAI。有效葉面積指數(shù)地面實(shí)測(cè)方法如圖3所示。

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

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

1.4 LiDAR特征變量提取

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

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

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

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

(1)

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

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

(2)

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

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

(3)

(4)

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

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

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

I=I0exp(-m)。

(5)

植物冠層消光系數(shù)直接決定其對(duì)輻射的截獲量，因?yàn)榇┻^(guò)的介質(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ù)激光雷達(dá)計(jì)算葉面積指數(shù)：

(10)

式中：θ為激光雷達(dá)掃描角。

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

(11)

式中：n為L(zhǎng)iDAR發(fā)射的總光點(diǎn)數(shù)；nground為L(zhǎng)iDAR穿透森林冠層后到達(dá)地面的光點(diǎn)數(shù)。

基于以上公式，可通過(guò)孔隙度模型法算出有效葉面積指數(shù)，使用地面實(shí)測(cè)有效葉面積指數(shù)對(duì)孔隙度模型法估測(cè)的有效葉面積指數(shù)進(jìn)行精度評(píng)價(jià)。

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

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

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

(12)

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

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

(13)

(14)

(15)

AIC=-lgL+2k。

(16)

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

2 結(jié)果與分析

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

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

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

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

圖4 地面實(shí)測(cè)eLAI與統(tǒng)計(jì)模型法預(yù)測(cè)eLAI交叉驗(yàn)證結(jié)果Fig.4 The cross-validation results of the filed-measured and statistical model predicted eLAIa為基于高度特征變量建模； b為基于高度+密度特征變量建模； c為基于高度+冠層容積特征變量建模； d為基于所有提取的特征變量(含強(qiáng)度特征變量)建模。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)計(jì)模型法預(yù)測(cè)有效葉面積指數(shù)模型交叉驗(yàn)證結(jié)果如圖4所示。僅通過(guò)LiDAR高度特征變量參與建模的精度最低，R2為0.38(rRMSE=54%)(圖4a)； 使用高度和密度特征變量組合參與建模的精度為R2=0.64(rRMSE=26%)(圖4b)； 通過(guò)高度和冠層容積特征變量組合與地面實(shí)測(cè)eLAI建立模型的精度為R2=0.61(rRMSE=28%)(圖4c)； 當(dāng)使用所有LiDAR特征變量(包括強(qiáng)度特征變量)建模時(shí)精度最高R2為0.74(rRMSE=23%)(圖4d)。加入其他類(lèi)型特征變量后，比僅通過(guò)高度特征變量建立模型的預(yù)測(cè)效果都有不同程度提高，其中，基于所有提取的特征變量建模預(yù)測(cè)效果最佳，通過(guò)高度和密度特征變量組合估測(cè)eLAI效果次之，最后是通過(guò)高度和冠層容積特征變量組合估測(cè)eLAI的精度。

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

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

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

圖5 不同Cover組地面實(shí)測(cè)eLAI與統(tǒng)計(jì)模型法預(yù)測(cè)eLAI交叉驗(yàn)證結(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 孔隙度模型法估測(cè)eLAI與地面實(shí)測(cè)eLAI擬合精度Fig.6 eLAI estimated by gap-fraction model

3組樣地冠層垂直結(jié)構(gòu)分布剖面【計(jì)算方法詳見(jiàn)Coops等(2007)】如圖8所示。由3組典型樣地冠層高度分布(canopy height distribution，CHD)可知，第1組樣地激光雷達(dá)點(diǎn)密度最高處在3～5 m之間，第2組樣地激光雷達(dá)點(diǎn)密度最高處在7.2～10.5 m之間，第3組樣地激光雷達(dá)點(diǎn)密度最高處在3～11 m之間，3組樣地平均樹(shù)高依次增高，冠幅增大(a1—c1)； 由3組樣地枝葉剖面(FP)隨冠層高度變化情況可知，第1組樣地枝葉剖面峰值最大在約3.5 m樹(shù)高處，第2組樣地枝葉剖面峰值最大在約8.5 m樹(shù)高處，第3組樣地枝葉剖面峰值最大在約9.0 m樹(shù)高處，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)、激光雷達(dá)點(diǎn)云(a3—c3)及激光雷達(dá)點(diǎn)云剖面(含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

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

圖8 3個(gè)典型樣地冠層高度分布(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)

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

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

4 結(jié)論

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