耿志偉，薛 偉，王海濱，張華超，，郭德禹

（1.東北林業(yè)大學 工程技術(shù)學院，黑龍江 哈爾濱 150040；2.白河林業(yè)局，吉林省 延邊朝鮮族自治州 133000）

植被聚集指數(shù)是表征植被冠層元素空間分布特征的重要參數(shù)，描述了真實葉片空間分布與理想狀態(tài)下隨機分布的偏離程度[1]。精確估計植被聚集指數(shù)，對研究植被生物量、陸地生態(tài)系統(tǒng)碳循環(huán)等具有重要意義[2]。此外，聚集指數(shù)還可以將有效葉面積指數(shù)轉(zhuǎn)化為真實葉面積指數(shù)，因此也被定義為植被有效葉面積指數(shù)與真實葉面積指數(shù)之比[3-5]。當葉片隨機分布時，聚集指數(shù)為1；當葉片聚集分布時，聚集指數(shù)小于1；當葉片規(guī)則分布時，聚集指數(shù)大于1。自然界中大多數(shù)葉片呈現(xiàn)聚集分布狀態(tài)。

目前，獲取植被聚集指數(shù)的方法主要有儀器實地測量和遙感技術(shù)反演。實地測量一般采用TRAC 儀器（Tracing radiation and architecture of canopies，TRAC）和數(shù)字半球攝影（Digital hemispheric photography，DHP）。但兩者都需在光照條件下工作，TRAC 設(shè)備應(yīng)在陽光直射下進行手動操作[6]，DHP 需設(shè)置最佳曝光時間并在漫射光下獲取影像[7]。傳統(tǒng)光學遙感反演植被聚集指數(shù)的方法多采用中高空間分辨率遙感影像（MODIS、TM 影像）和多角度多光譜數(shù)據(jù)（POLDER、MISR 數(shù)據(jù)）[8-9]。然而由于影像熱點和冷點易受植被各向異性、太陽角度等因素影響，聚集指數(shù)估算具有較大不確定性，因此該技術(shù)常用于全球或區(qū)域尺度的聚集指數(shù)粗略計算[10]。

相比之下，激光雷達（Light detection and ranging,LiDAR）受天氣影響小，可快速精準獲取植被三維結(jié)構(gòu)信息，已廣泛應(yīng)用于植被結(jié)構(gòu)參數(shù)反演[11]。針對星載激光雷達，Yang 等[12]聯(lián)合中分辨率遙感影像提出了植被冠間聚集指數(shù)物理反演方法。針對機載激光雷達，Thomas 等[13]利用地面實測聚集指數(shù)和機載點云數(shù)據(jù)提取的幾何參量建立經(jīng)驗?zāi)Ｐ停籋u 等[14]引入路徑長度分布函數(shù)表征樹冠形狀，建立了基于機載激光雷達的聚集指數(shù)反演模型。然而，星載激光雷達全波形數(shù)據(jù)無法直接提供大光斑內(nèi)地物的精細分布，機載激光雷達聚集指數(shù)提取精度極易受點云密度影響。

與星載、機載激光雷達相比，地基激光雷達（Terrestrial laser scanning,TLS）具有空間分辨率高、光斑尺寸小等特點，能夠精確刻畫植被三維結(jié)構(gòu)，已成為提取植被聚集指數(shù)的重要手段[15-22]?，F(xiàn)有地基激光雷達反演植被聚集指數(shù)的方法可分為三類：一是基于二維圖像的強度法。如Zhao 等[16]將激光雷達強度數(shù)據(jù)通過半球投影生成強度圖像，并應(yīng)用于植被聚集指數(shù)提取。二是基于三維體素的幾何法。通常將激光點云按照一定規(guī)則劃分為規(guī)則體素獲取其間隙率，并假設(shè)體素內(nèi)的點云均勻分布[17-20]。然而，聚集指數(shù)反演精度和效率與規(guī)則體素的劃分方法有關(guān)，如何依據(jù)葉片特征、點云密度等來選取最優(yōu)體素仍是該方法的一大難點。如Zheng 等[19]提出徑向半球面切片法進行體素劃分，較好地實現(xiàn)了冠層間隙率的提取；在此基礎(chǔ)上，Ma 等[20]分別探討了徑向半球面切片法、定向切片法等體素劃分方法在聚集指數(shù)反演中的優(yōu)劣，并結(jié)合間隙大小分布理論實現(xiàn)樣方尺度聚集指數(shù)反演。三是基于點的幾何法，常見于單木聚集指數(shù)估算研究[21-23]。如Yun 等[21]基于地基點云數(shù)據(jù)建立三角網(wǎng)模型，通過識別葉片點來計算冠層真實葉面積指數(shù)；Li 等[22]結(jié)合激光點云的多回波信息與領(lǐng)域連通法提取冠層間隙，利用間隙大小分步法實現(xiàn)了針對單站地基掃描的植被聚集指數(shù)快速提取。

相較而言，基于點的幾何法對點云信息的利用較為充足，減少了資源浪費。然而，現(xiàn)階段基于點的地基激光雷達聚集指數(shù)反演研究相對較少。因此，本研究提出了一種利用多站地TLS 數(shù)據(jù)反演植被聚集指數(shù)的方法——基于三維間隙的聚集指數(shù)提取方法。該方法利用TLS 獲得的點云數(shù)據(jù)構(gòu)建Delaunay 四面體格網(wǎng)模型，通過計算四面體體積統(tǒng)計得到三維間隙累積分布函數(shù)，最后采用間隙大小分布法反演得到植被冠層聚集指數(shù)。該方法相比基于二維間隙的方法更能反映植被的空間聚集特征。同時多站TLS 數(shù)據(jù)保證了植被聚集指數(shù)的高精度提取，為地基激光雷達提取樣方尺度植被聚集指數(shù)提供一種新思路。

1 數(shù)據(jù)獲取

本研究收集了位于內(nèi)蒙古根河林區(qū)（38°29′N～38°59′ N，100°12′ E～100°20′ E）和甘肅張掖林區(qū)（50°01′ N～50°51′ N，121°20′ E～121°38′ E）的20 塊30 m×30 m 植被樣方的地基激光雷達數(shù)據(jù)和野外實測驗證數(shù)據(jù)（圖1）。2 個林區(qū)的地面實驗分別在2016年8月和2012年7月展開。各樣方內(nèi)植被茂密，死樹較少，地勢相對平坦。

1.1 地基激光雷達數(shù)據(jù)

2 個林區(qū)使用的地面激光掃描儀都為RIEGL VZ-1000。最大脈沖頻率和射程分別為300 kHz 和1 400 m，測距精度為2 mm@50 m 處。儀器架設(shè)高度約為1.8 m，掃描視野為360°×100°。為減少枝葉遮擋效應(yīng)、精確獲取冠層完整結(jié)構(gòu)，本研究對每個樣地均進行5 個測站的掃描（圖2）。同時在每個樣方中布設(shè)6 個固定黑白標靶，方便在后續(xù)操作中將各測站數(shù)據(jù)統(tǒng)一在一個坐標系下。樣方中心測站能夠看到所有標靶，其余測站至少掃描3 個標靶。地基激光掃描在清空無風的條件下進行。

1.2 野外實測數(shù)據(jù)

圖1 研究區(qū)域 Fig.1 The position of the study area

地面實測在地基掃描樣方中同步開展，測量方式為TRAC 儀器。測量在晴朗無云條件下開展，且太陽光較為穩(wěn)定。操作人員手持儀器以約0.3 m/s速度勻速前進，行走方向垂直于太陽入射方向，每個樣地進行150～200 m 樣線采集（圖2）。

圖2 樣方掃描基站分布及TRAC 測量樣線（紅色虛線）Fig.2 The distribution of scan positions and TRAC sampling line (red dotted line)

2 植被聚集指數(shù)提取方法

本研究通過地基激光掃描數(shù)據(jù)預(yù)處理、植被點云提取、四面體格網(wǎng)模型建立、間隙累積分布函數(shù)計算、間隙大小分布法等過程實現(xiàn)植被冠層聚集指數(shù)的反演，具體流程如圖3所示。

2.1 點云數(shù)據(jù)預(yù)處理

圖3 植被聚集指數(shù)反演技術(shù)流程Fig.3 Workflow diagram for clumping index retrieval

TLS 點云數(shù)據(jù)預(yù)處理工作主要包括點云配準和點云去噪。點云配準的目的是將多站地基數(shù)據(jù)轉(zhuǎn)換到同一坐標系下，得到目標整體數(shù)據(jù)集。本研究選用基于公共標靶的配準方法，將標靶中心作為兩個測站數(shù)據(jù)的公共點，通過自動添加約束條件實現(xiàn)基于標靶的配準。該方法與手工選取同名特征點的拼接方法相比，精度高且效率快。最終20 個樣地的配準誤差均不超過3 mm，滿足后續(xù)實驗精度要求。然而，受系統(tǒng)誤差、環(huán)境干擾等影響，配準后的點云數(shù)據(jù)中仍存在明顯噪聲。因此，本文采用基于局部空間統(tǒng)計的去噪算法對配準后點云進行處理。該方法將局部點密度與整體點密度相差較大的點標記為噪聲點，實現(xiàn)噪聲去除，在保留目標有效細節(jié)特征的同時，保證了點云數(shù)據(jù)密度相對均勻。

2.2 植被點云提取

由于研究對象為植被點云，因此采用Zhang等[24]提出的布料模擬濾波算法實現(xiàn)地面信號與植被信號的快速精確分離。具體實現(xiàn)過程為：首先將整體點云進行反置，隨后將一塊剛性布料覆蓋在翻轉(zhuǎn)后的表面上；模擬布料在重力作用下的下落過程，確定布料節(jié)點的位置，即地形表面位置；最終對比地形表面和原始點云，分離地面點云與植被點云。

基于濾波獲得的地面點云，使用一次樣條有限元法進行空間內(nèi)插[25]，建立空間分辨率為0.2 m 的數(shù)字地形模型（Digital elevation model,DEM）。之后，將點云濾波后得到的植被點云以X、Y坐標為依據(jù)，根據(jù)生成的DEM 影像分辨率尺寸劃分到不同格網(wǎng)。將落在格網(wǎng)內(nèi)的植被點高程減去對應(yīng)的DEM 格網(wǎng)值，即得到歸一化植被點云。此時植被點云中的點高度全部轉(zhuǎn)換為其相對地面的高程，即植被高度。

2.3 四面體格網(wǎng)模型構(gòu)建

為獲取植被冠層的間隙累計分布函數(shù)，本文提出通過構(gòu)建Delaunay 四面體格網(wǎng)模型來對植被點云進行處理。四面體格網(wǎng)模型能夠保留原始點云數(shù)據(jù)，同時精確表達植被與間隙的空間分布。但常用的構(gòu)網(wǎng)方法普遍存在效率低、速度慢等缺點，極其不適用于密度高、數(shù)據(jù)量大的地基點云數(shù)據(jù)。因此本研究采用分而治之算法進行四面體格網(wǎng)模型構(gòu)建。

首先，將所有的植被點云進行虛擬三維規(guī)則格網(wǎng)劃分，保證每個格網(wǎng)中至少包含一個數(shù)據(jù)點；并設(shè)定格網(wǎng)點云密度閾值判定各格網(wǎng)是否為密集格網(wǎng)。接著，在格網(wǎng)內(nèi)選取種子點，采用局部優(yōu)化算法（Local optimization procedure,LOP）[26-27]，即選擇最小角最大法連接種子點，建立得到初始四面體。在此基礎(chǔ)上依次加入非密集格網(wǎng)內(nèi)的點，直至所有非密集格網(wǎng)內(nèi)的點都連接完成（圖4）。然后處理數(shù)據(jù)密集三維格網(wǎng)，先在密集格網(wǎng)內(nèi)部構(gòu)建Delaunay 四面體，保證密集格網(wǎng)內(nèi)部點均被連接，再將新生成的內(nèi)部Delaunay 四面體與其外部四面體格網(wǎng)進行合并。直至遍歷所有密集格網(wǎng)點云數(shù)據(jù)，形成一個完整四面體網(wǎng)絡(luò)。最終，構(gòu)建的四面體格網(wǎng)模型形成了大量小四面體，所有小四面體一起構(gòu)成植被點云凸包。

圖4 Delaunay 四面體格網(wǎng)生成Fig.4 Delaunay tetrahedral network generation

該四面體格網(wǎng)模型滿足：1）網(wǎng)格模型中的每個面均為平面凸三角形；2）網(wǎng)格中每條邊有且只有2 個三角形共有；3）三角形組成的四面體之間不存在重疊或遺漏，且不包含植被點云以外的體積。對于任意一個四面體，設(shè)點O、A、B、C為其4 個頂點，其三維空間坐標分別為（XO，YO，ZO）、（XA，YA，ZA）、（XB，YB，ZB）、（XC，YC，ZC），則四面體TOABC的體積可計算得到（式（1））。

所有小四面體的體積累加起來即為整個植被點云凸包的體積。將小四面體體積除以凸包體積，統(tǒng)計不同體積的四面體出現(xiàn)的概率，即得到間隙累積分布函數(shù)Fm(λ)。該函數(shù)定義為體積大于λ的所有間隙所占的比例。

2.4 聚集指數(shù)計算

將基于四面體格網(wǎng)模型計算得到的間隙累積分布函數(shù)，輸入到Chen 和Chair 提出的間隙大小分布法（CC 算法）中[6]，計算得到聚集指數(shù)。該方法將葉片隨機分布下的間隙累積分布函數(shù)與實際分布函數(shù)進行對比，經(jīng)過逐步迭代移除大間隙的方式定量衡量植被聚集效應(yīng)。對于葉片隨機分布的冠層，其間隙累積分布函數(shù)為：

式中：Fr(λ)是隨機分布冠層的間隙累積分布函數(shù)。Lp是投影葉面積指數(shù)，其初值為-ln[Fmr(0)]，其中Fm(0)是實際總間隙率。Wp是葉片投影后的平均寬度，其計算方式為：

當實際間隙出現(xiàn)概率超過Fr(λ)時，則被認為是非隨機大間隙，將其移除。每刪除一輪大間隙，重新計算間隙累積分布函數(shù)Fmr(λ)，并通過新的-ln[Fmr(0)]計算新的Lp。經(jīng)多次迭代后，間隙累積分布函數(shù)Fmr(λ)將接近隨機狀態(tài)下的間隙累積分布函數(shù)。迭代終止條件為Lp增長幅度小于規(guī)定閾值，或Fmr(λ)的一部分落到Fr(λ)曲線下面。最終該算法計算得到的聚集指數(shù)表達式[28]為：

式中：Fm(0)是原始冠層總間隙率，F(xiàn)mr(0)是大間隙被移除后的冠層總間隙率。

結(jié)合四面體格網(wǎng)模型和間隙大小分布法的植被聚集指數(shù)反演方法，以三維分布的冠層間隙形式消除了二維間隙樣線內(nèi)部的不均勻性，并通過對比聚集分布和隨機分布下三維間隙累積分布函數(shù)，最終達到高精度反演植被聚集指數(shù)的目的。

3 結(jié)果與分析

3.1 點云數(shù)據(jù)處理

以根河某一樣地數(shù)據(jù)為例，使用布料濾波算法對預(yù)處理后點云數(shù)據(jù)（圖5a）進行地面和植被的分離，結(jié)果如圖5b 所示，紅色表示地面點云，綠色代表植被點云。經(jīng)地面點內(nèi)插與植被點云相對高程計算，最終得到歸一化植被點云，如圖5c所示?？梢钥闯?，布料濾波實現(xiàn)了復(fù)雜地形下地面點和非地面點的精確分離。然而對于植被下層存在的雜亂點，濾波結(jié)果仍存在一些誤差，這與樣方中的低矮植被有關(guān)。低矮植被通常生長濃密且分布均一，導(dǎo)致激光脈沖到達地表概率降低。濾波處理時，低矮植被點云可能被判別為地面點，可選擇人工修正來進一步處理。此外，由于本文主要研究對象為森林冠層，低矮植被在樣方植被中占比較小，且野外TRAC 儀器實測時測量高度一般高于低矮植被高度，因此低矮植被對樣方內(nèi)聚集指數(shù)反演精度影響不大。

圖5 點云數(shù)據(jù)處理結(jié)果Fig.5 TLS point cloud data processing

3.2 植被聚集指數(shù)反演

3.2.1 植被點云四面體格網(wǎng)模型

將歸一化植被點云劃分為三維規(guī)則格網(wǎng)（圖6a），遍歷所有格網(wǎng)內(nèi)植被點云，構(gòu)建Delaunay四面體格網(wǎng)模型（圖6b）。可明顯看出，樹冠內(nèi)的植被點云密集，生成的三角面片極多，四面體體積也相對較小。同時，葉片的聚集分布導(dǎo)致了大間隙的出現(xiàn)，表現(xiàn)為四面體格網(wǎng)模型中較大體積四面體出現(xiàn)的概率增加。利用頂點點云三維坐標計算獲得所有小四面體體積，將其累加即為植被點云凸包體積（圖6c）。統(tǒng)計不同體積的四面體出現(xiàn)的概率，得到該植被冠層的間隙累積分布函數(shù)，如圖7所示。

圖6 植被點云四面體模型Fig.6 Delaunay tetrahedral network model of vegetation point cloud

本研究采用四面體模型來獲取植被三維間隙累積分布函數(shù)。傳統(tǒng)方法的間隙累計分布函數(shù)多基于間隙尺寸，即二維間隙。然而離散激光點云的二維間隙尺寸很難給出明確定義，且計算二維間隙會忽略第三維的聚集特征。因此，本研究得到的間隙累積分布函數(shù)更符合三維空間實際情況，為后續(xù)聚集指數(shù)反演提供了準確的間隙分布參量。

圖7 間隙累積分布函數(shù)Fig.7 Gap cumulative distribution function

3.2.2 聚集指數(shù)反演結(jié)果

通過多次迭代移除大間隙的方式，從三維間隙累積分布函數(shù)中計算得到植被聚集指數(shù)，間隙累積分布函數(shù)迭代前后的結(jié)果如圖8所示。可以看出，與理想隨機分布狀態(tài)相比，實際聚集分布下的植被冠層內(nèi)大間隙頻數(shù)明顯較多，小間隙頻數(shù)較少。經(jīng)多次迭代后，實際間隙累積分布函數(shù)與理想場景中幾乎吻合，但仍存在一定程度的抖動。這是因為理想隨機分布場景下的Fr(λ)是一個光滑曲線，而Fm(λ)是存在鋸齒狀的實測曲線。即使Fmr(λ)滿足迭代終止條件，也不可能與Fr(λ)完全重合。

圖8 間隙累積分布函數(shù)迭代結(jié)果Fig.8 Gap cumulative distribution functions before and after iteration

為驗證該方法提取植被聚集指數(shù)的精度，選取決定系數(shù)R2、相對均方根誤差rRMSE作為指標，對20 個樣方內(nèi)TLS 數(shù)據(jù)提取的聚集指數(shù)結(jié)果進行評價。圖9為實測聚集指數(shù)與TLS 反演的聚集指數(shù)的相關(guān)關(guān)系圖。由圖9可得，該方法反演的聚集指數(shù)與實測聚集指數(shù)表現(xiàn)出較好的正相關(guān)性（R2=0.74，rRMSE=8.94%）。然而，由三維間隙計算得到的聚集指數(shù)與TRAC 儀器測量值相比，存在略微低估。其原因除了受間隙閾值設(shè)置的部分影響外，更主要在于兩者原理不同。TRAC 儀器基于觀測樣線采集的是二維間隙，而多站地基激光雷達通過大量脈沖覆蓋整個三維空間?；谒拿骟w模型的三維間隙比二維間隙更能反映出植被的真實聚集情況。由此可推論，目前植被聚集指數(shù)實地測量方法得到的測量值仍偏低，與Ryu 等[29]的結(jié)論一致。

圖9 TLS 提取的聚集指數(shù)與實測聚集指數(shù)相關(guān)關(guān)系Fig.9 Relationship between the TLS-derived CIs and the TRAC-measured CIs

4 結(jié)論與討論

本研究經(jīng)多站地基激光雷達掃描，獲取了內(nèi)蒙古根河林區(qū)和甘肅張掖林區(qū)多個樣地點云數(shù)據(jù)，然后通過數(shù)據(jù)預(yù)處理、植被點云提取、四面體格網(wǎng)模型構(gòu)建和間隙大小分布法，實現(xiàn)了樣方內(nèi)植被聚集指數(shù)反演。主要結(jié)論如下：

1）首次提出利用四面體格網(wǎng)模型來獲取植被三維間隙，并計算得到植被間隙累積分布函數(shù)。該模型充分利用了激光雷達精確獲取目標三維信息的能力，得到的間隙累積分布函數(shù)也更能體現(xiàn)植被三維結(jié)構(gòu)信息。

2）將四面體格網(wǎng)模型與經(jīng)典的間隙大小分布法相結(jié)合，通過逐步移除大間隙的方式，使得聚集分布的森林冠層逐漸接近隨機分布。

3）經(jīng)TRAC 儀器測量的聚集指數(shù)驗證，該算法能夠較好地提取植被冠層聚集指數(shù)（R2=0.74，rRMSE=8.94%，n=20），但存在略微低估，其原因主要為計算原理不同。二者采用的雖然都是間隙大小分步法，但TRAC 儀器利用樣線采集得到二維間隙，本文方法采用TLS 點云獲取得到三維間隙。

總體而言，本文基于高密度地基激光雷達數(shù)據(jù)，提出了利用三維間隙反演植被聚集指數(shù)的方法，能夠更加真實地反映出植被冠層的聚集情況。然而，該方法對數(shù)據(jù)質(zhì)量要求較高，要求點云密度相對均勻，且盡量避免目標的前后遮擋。若點云密度相差較大，由低密度處的植被點云生成的四面體極大可能會被標記為三維間隙，導(dǎo)致聚集指數(shù)估算不準。因此如何將該方法遷移到存在明顯遮擋的單站地基激光雷達或機載激光雷達數(shù)據(jù)，克服該方法在點云密度不均情況下的局限性，是下一步的研究方向。此外，探討基于二維間隙與三維間隙提取的聚集指數(shù)間的根本差異，也是本研究未來需重點開展的工作。