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

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

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

0 引言

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

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

1 研究區(qū)和數(shù)據(jù)來源

1.1 研究區(qū)概況

研究的案例樹木為位于福州大學(xué)旗山校區(qū)校園內(nèi)部道路兩側(cè)的芒果樹，如圖1所示.

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

1.2 數(shù)據(jù)獲取及預(yù)處理

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

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

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

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

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

2 研究方法

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

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

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

(1)

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

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

(2)

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

(3)

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

(4)

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

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

(5)

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

P=Np×res2

(6)

這里，Np為投影面體素個數(shù).

2.2 冠層體積變化分析方法

各種冠層體積估計方法有各自的優(yōu)缺點，適用于不同的冠型. 根據(jù)冠層體積的定義，目前常用的冠層體積估計方法是基于凸殼或逐層疊加臺體的方法. 圖5(a)為常用冠層體積估計模型凸殼模型; 圖5(b)～(c)分別為α=1的Alpha shape 3D模型及其切片形式. 由Alpha shape 3D構(gòu)建原理可知，當參數(shù)α=∞時，網(wǎng)格模型將為凸殼模型. 取冠層一定厚度切片點云，根據(jù)凸殼模型與Alpha shape 3D算法原理[13]形成圖6(a)～(b)兩種方法線框圖以及圖6(c)的兩種方法模型疊加圖，結(jié)合樹的形態(tài)結(jié)構(gòu)特點的不規(guī)則性，Alpha Shape 3D切片模型更加貼近于切片外圍輪廓. 相關(guān)研究說明通過切片形式能夠提高冠層體積估計精度[9, 14]，采用切片化Alpha shape 3D模型估算冠層體積，根據(jù)冠型特點及點云分散情況將冠層點云按一定厚度水平切片后，逐層構(gòu)建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語言圖形庫對臺風(fēng)后非枝干冠層點云構(gòu)建Alpha shape 3D網(wǎng)格模型，并輸出該模型. 由于輸出模型三角面片未帶面片法向量信息，使用開源點云處理軟件cloudcompare對網(wǎng)格模型進行表面重采樣，產(chǎn)生模型表面高密度點云，并重新計算點云法向量. 計算臺風(fēng)前點云與重建網(wǎng)格模型間的距離(cloud to mesh distance，C2M distance)[15]，通過直通濾波器對臺風(fēng)前冠層點云C2M distance屬性進行濾波，獲取網(wǎng)格外冠層點云，即冠層變化點云. 基于提取的冠層變化點云，分析變化體積量及變化點云分布.

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

3 結(jié)果與分析

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

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

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

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

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

表1 冠幅及冠層投影面積

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

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

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

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

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

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

4 結(jié)語

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