基于地面激光雷達點云的樹木冠層變化定量分析

楊水榮，黃洪宇，唐麗玉，陳崇成

(福州大學空間數據挖掘與信息共享教育部重點實驗室，福州大學地理空間信息技術國家地方聯合工程研究中心，福建 福州 350108)

0 引言

樹木空間結構及動態變化信息，對于樹木經營管理、 生態環境建模和植物表型學的研究等具有重要意義. 森林尺度的樹木變化檢測以森林樹木在物候期及人工干擾等情況下生物量及結構變化研究為主[1]. 目前對于單木變化檢測大部分側重于樹木生長生理變化，借助微樹芯、 生長錐[2]等方式進行單木對象的變化檢測對樹木有一定破壞性，且樹木徑向生長測量精度不高. 地面激光掃描儀(terrestrial laser scanner，TLS)具有獲取單木尺度的樹木形態結構信息的能力，能獲取高精度的單木枝干、 冠層結構信息和對樹木無破壞性獲取等優點[3]. 目前，基于TLS點云數據的單木尺度樹木冠層變化檢測主要從兩方面展開. 一方面，通過空間統計分析對不同階段植物冠層點云的百分位高度進行比較，量化植物在特定狀態或時間段冠層移動及節律性變化[4]. 該方式雖然簡單，但對實驗條件要求嚴格，需排除風擾動等外界影響因素，難以在常規室外環境中實施，無法體現樹冠在表型結構上的變化. 另一方面，通過重建不同時期單木的定量結構模型(quantitative structure model，QSM)[5]或基于體素的枝干結構模型[6]，監測冠層結構參數變化情況，反映樹木冠層擴展[7]及枝干體積量變化[8]. 該方式只適用于落葉或無葉樹木，不適用于含有茂盛枝葉的樹木對象. 有葉冠層可采用冠層體積描述變化，樹木冠層體積包括了冠層輪廓以內所占的體積[9]. 目前已有多種方式[10]進行冠層體積估計. 由于樹木冠層結構復雜，規則幾何體無法描述真實輪廓，精準測定冠層體積較為困難[11]. 冠層內部結構復雜程度對于地面激光雷達采集數據質量影響較大，且冠層茂密存在遮擋，導致掃描獲取的冠層點云數據存在一定的缺失. 目前，研究主要從冠層體積定義出發，以逼近冠層輪廓的體積作為樹木冠層體積[12].

近年來，基于TLS點云數據的冠層變化檢測研究較少，本研究利用體素概括性，采用基于體素的冠層分析方法從不同參考角度對冠層進行組分變化分析； 根據冠層體積定義，采用切片后點云逐層構成Alpha shape3d模型進行冠層體積估計，通過冠層外圍輪廓與點云間的距離關系分析冠層變化情況. 以臺風“瑪利亞”過境前后園林樹種芒果樹為例，探究冠層變化定量分析方法的有效性.

1 研究區和數據來源

1.1 研究區概況

研究的案例樹木為位于福州大學旗山校區校園內部道路兩側的芒果樹，如圖1所示.

圖1 研究區及研究對象Fig.1 The research area and the research object

1.2 數據獲取及預處理

使用Riegl VZ-400地面激光掃描儀，獲取樹木空間坐標信息及反射率屬性. 通過同名點拼接多站掃描數據. 根據現場環境，臺風前后掃描儀設站位置如圖2所示，掃描設站數均為5站，掃描角分辨率均為0.04°，垂直掃描角范圍為30°～130°，水平掃描角范圍為0°～360°，多站點云拼接標準誤差控制在0.004～0.005 m. 通過對兩期樹木點云構建KD-Tree數據結構，計算樹木點云近鄰點間間距，拼接后的原始樹木點云間平均間距約為0.003 m. 采用半徑濾波器(radius filter)進行噪聲濾波，過濾后樹木近鄰點云平均間距為0.003 7 m(詳見圖3(a)). 由圖3(b)窗口圖亦可看出臺風對樹木冠層的破壞程度，冠層外圍葉片點云少，果實點云消失(紅點表示臺風后點云； 黑點表示臺風前點云).

圖2 臺風前后測站空間位置示意圖Fig.2 Schematic diagram of scanner location before and after typhoon

圖3 樹木兩期點云配準及冠層受損情況Fig.3 Two-phase tree point cloud registration and canopy damage

跟蹤樹木受臺風影響時間短，成熟期樹木枝干生長變化可忽略不計，分析冠層在臺風期間變化情況，進行樹木點云的枝葉分割. 由于地面激光掃描儀獲取的樹木點云數據僅含坐標及反射率信息，參考點云反射率信息，通過多人對臺風前后樹木對象進行分揀枝干點云與葉片點云并評估分揀結果，臺風前后枝葉分類結果如圖4所示.

圖4 枝葉點云分揀結果Fig.4 Branches and leaves point cloud sorting results

2 研究方法

2.1 冠層組分分布變化分析方法

體素作為三維空間分割最小基本單位，對三維空間的概括性，在一定程度上不受噪聲及遮擋現象的影響.

樹木點云體素化前，首先提取樹木根部地面以上5 cm切片點云，利用最小二乘法擬合圓，求取圓心，并以根節點為空間坐標原點，從不同參考角度分析樹木整體變化情況. 點云體素化公式為：

(1)

其中：res指的是體素分辨率； coord(x, y, z)為點云三維坐標； round()函數為四舍五入函數；V(x, y, z)為體素中心坐標.

針對體素分布情況，從不同參考軸及角度分析，首先計算體素相對距離：

(2)

式中：v為點云數量是n的體素集； (xvi,yvi,zvi)為體素集內體素中點坐標；p(xp,yp,zp)為定義參考點坐標. 體素相對坐標軸距離的計算公式如下式所示：

(3)

這里，a為定義參考坐標軸，當a為x軸，xvi=0； 當a為y或z，yvi或zvi=0. 根據體素向量與坐標軸向量的內積公式，體素相對坐標軸角度的計算公式如下式所示：

(4)

其中：β1為體素向量；β2為選定參考坐標軸向量，本次實驗取|β2|=100，當a為x軸時，β2為x軸向量，當a為y或z軸時，β2為y或z軸向量.

冠層體素投影圖可以反映冠層部分受生長、 環境等因素的影響情況. 兩期的點云處理后，對得到的葉片點云進行體素化，然后進行平面投影，每個投影像素包含點云個數Np和體素個數Nv，分析其體素投影密度分布、 冠幅和投影面積變化情況，以下為冠幅計算公式：

(5)

其中：ΔXp和ΔYp分別為投影體素中橫、 縱坐標的極值差. 投影面積計算公式為：

P=Np×res2

(6)

這里，Np為投影面體素個數.

2.2 冠層體積變化分析方法

各種冠層體積估計方法有各自的優缺點，適用于不同的冠型. 根據冠層體積的定義，目前常用的冠層體積估計方法是基于凸殼或逐層疊加臺體的方法. 圖5(a)為常用冠層體積估計模型凸殼模型; 圖5(b)～(c)分別為α=1的Alpha shape 3D模型及其切片形式. 由Alpha shape 3D構建原理可知，當參數α=∞時，網格模型將為凸殼模型. 取冠層一定厚度切片點云，根據凸殼模型與Alpha shape 3D算法原理[13]形成圖6(a)～(b)兩種方法線框圖以及圖6(c)的兩種方法模型疊加圖，結合樹的形態結構特點的不規則性，Alpha Shape 3D切片模型更加貼近于切片外圍輪廓. 相關研究說明通過切片形式能夠提高冠層體積估計精度[9, 14]，采用切片化Alpha shape 3D模型估算冠層體積，根據冠型特點及點云分散情況將冠層點云按一定厚度水平切片后，逐層構建Alpha shape 3D模型，并累加模型體積作為冠層體積估計值.

圖5 三種冠層體積估計模型Fig.5 Three canopy volume estimation methods

圖6 Alpha shape 3D與凸殼法方法對比Fig.6 The comparison between convex hull method and Alpha shape 3D

進一步分析冠層變化點云，分析的基本流程如圖7所示. 首先，使用R語言圖形庫對臺風后非枝干冠層點云構建Alpha shape 3D網格模型，并輸出該模型. 由于輸出模型三角面片未帶面片法向量信息，使用開源點云處理軟件cloudcompare對網格模型進行表面重采樣，產生模型表面高密度點云，并重新計算點云法向量. 計算臺風前點云與重建網格模型間的距離(cloud to mesh distance，C2M distance)[15]，通過直通濾波器對臺風前冠層點云C2M distance屬性進行濾波，獲取網格外冠層點云，即冠層變化點云. 基于提取的冠層變化點云，分析變化體積量及變化點云分布.

圖7 基于網格模型的冠層變化分析基本流程Fig.7 Basic processing steps of canopy changes analysis based on mesh model

3 結果與分析

根據公式(3)，結合高斯核密度函數統計臺風前后體素相對z軸的分布變化情況(具體詳見圖8). 由圖8可知，臺風前后樹木整體體素密度分布趨勢相近，相對z軸0.3～1.5 m區間內，臺風后體素密度略高于臺風前; 相對于z軸2.1 m外，臺風前體素密度大于臺風后體素密度. 由此可知，臺風過后，遠離主干的冠層外部組分受到臺風的影響較靠近主干的冠層內部組分要大，冠層內部組分變得更加緊致，而外冠層對應的數據點減少. 根據體素相對軸角度公式，設方位角0°為測站的x軸朝向，并將空間坐標原點移至樹木根節點，臺風前后體素變化詳見圖9. 由圖9可知，方位角55°～165°區間內，樹木體素密度增加，該方位角現場區域相鄰為另一株樹木，方位角165°～280°區間內臺風后密度降低. 這表明同一棵樹在不同方位(水平和垂直方向)受到的風力作用的大小有差異，并體現了風場與環境(地形、 建筑、 樹木分布)相互作用的復雜性.

圖8 臺風前后體素密度變化情況(相對z軸距離)Fig.8 Changes of canopy voxel density before and after typhoon (relative z axis distance)

圖9 臺風前后體素密度變化情況(相對x軸角度)Fig.9 Changes of canopy voxel density before and after typhoon (relative x axis angle)

圖10 臺風前后冠層變化情況(x-y平面)Fig.10 Changes of crown projection area and crown width before and after the typhoon(x-y plane)

根據樹木冠層組分的變化情況，利用枝葉分割提取非枝干點云，體素化后投影至x-y平面，以每個投影像素中含有的體素個數作為其像素密度值，將冠層沿垂直高度三等分后，對冠層上、 中和下部冠層分別取臺風前后冠層2.787、 4.186和5.687 m處20 cm厚度的點云切片，通過對臺風前后的切片體素化后進行差分，分析其變化情況, 如圖10所示. 其中x代表體素柵格距x軸距離，y代表體素柵格距y軸距離(單位：m). 將冠層體素投影像素分為五種變化狀況：密度未變化體素，即臺風前后同一位置投影像素密度值未發生變化； 密度減少像素，即臺風后同一位置投影像素密度減少； 密度增大像素，即臺風后同一位置投影像素密度增大，可由臺風擾動造成的冠層組分位移造成； 臺風后形成像素，即由于受臺風影響造成冠層組分位移，這些像素來源可由臺風后消失像素及密度變化像素產生； 臺風后消失像素，即臺風后同一位置投影像素消失. 由圖10可知，切片01損失像素個數最多，且隨高度增加切片層消失體素個數逐漸降低，而臺風后形成像素個數隨高度增加而增加，從臺風后形成像素個數與消失體素個數可以看出，由于冠層下部冠層枝葉更為茂盛，臺風對冠層下部損壞較嚴重，而冠層上部受風擾動造成冠層位移較大. 密度減少像素個數變化趨勢與消失像素個數趨勢一致.

表1 冠幅及冠層投影面積

利用公式(5)～(6)分析臺風前后樹木冠層的投影面積及冠幅變化情況, 詳見表1，從中可以看出臺風過境后樹木冠幅和冠層投影面積均減少.

樹冠結構參數主要包括樹木冠層枝條數量、 冠層體積和葉面積等. 由于樹木對象的冠層稠密度對地面激光掃描儀獲取樹木點云數據質量的影響，冠層枝條數量除現場驗證外，并無更準確的方式. 本研究主要分析臺風前后樹木冠層體積變化情況.

對于冠層樹木的冠層體積估計方法，基于切片形式的Alpha shape 3D冠層模型相較于其他模型更接近冠層外圍輪廓，因此本實驗采用基于切片形式的Alpha shape 3D冠層估計方法對臺風前后樹木冠層進行體積估計，冠層體積分別為40.678 3、 39.010 9 m3，冠層體積減少量為1.667 4 m3.

基于網格模型提取臺風前后變化點云，實驗獲取臺風前后冠層變化點云如圖11所示.

圖11 變化點云提取及分布情況Fig.11 Extraction and distribution of change point clouds

通過計算點云與主干平面位置夾角，分析變化點云分布情況(詳見圖11(b))，冠層變化點云大部分處于距離冠層主干2.0～2.5 m范圍內. 考慮變化點云分散程度，選用基于體素的變化體積量估計方式，體素尺寸選擇0.05 m，計算得變化體積量為0.786 5 m3. 圖11(c)為變化點云體積剖面圖，由圖可知，臺風過后，冠層中下位置體積量變化較大.

4 結語

樹木冠層結構及組分分布的變化影響著林木生長發育，而樹木冠層的復雜性給冠層變化檢測帶來一定的困難. 本研究從樹木冠層點云本身特點入手，提出從點云再到體素兩個尺度上的樹木冠層變化定量分析方案. 基于體素具有空間概括性的特點，利用體素化方法，從體素相對空間位置、 體素相對軸距離、 體素相對軸夾角和體素平面投影的角度，分析冠層組分分布變化情況. 鑒于目前未出現一種普適性、 魯棒性較強的冠層體積估計方案，提出基于切片形式的Alpha shape 3D冠層體積估計方法. 該方法計算原理更接近于冠層實際輪廓. 在構建Alpha shape 3D模型基礎上，結合冠層網格模型與點云間距離關系，提取冠層變化點云，并從方位角和垂直方向分析變化點云分布情況. 以受臺風影響的芒果景觀樹冠層變化分析為例，分析冠層在臺風影響下冠層組分分布、 冠層體積變化情況，并分析受臺風影響變化點云分布情況. 該定量化分析方案可應用于分析植物生長過程的冠層變化情況，有利于探索冠層內各組分生長發育的空間模式，目前僅針對單木對象，未來將進一步擴大研究尺度. 由于不同樹木冠層復雜程度不同，采用的冠層體積估計方法有所差異. 鑒于冠層體積精準測量仍是研究難點，故本研究所提出的切片形式Alpha shape 3D冠層估計方法在不同樹木冠層適用性有待驗證，未來將結合光譜信息等屬性進行精度評價.