林怡， 曹宇杰，

（1.同濟(jì)大學(xué) 測(cè)繪與地理信息學(xué)院，上海 200092；2.同濟(jì)大學(xué) 遙感技術(shù)應(yīng)用研究中心，上海 200092）

激光雷達(dá)技術(shù)（light detection and ranging，LiDAR）作為一種主動(dòng)遙感技術(shù)已廣泛應(yīng)用于各類林業(yè)資源調(diào)查任務(wù)中［1］。其在信息更新頻率、采集成本、數(shù)據(jù)精度、量化分析等方面比傳統(tǒng)人工測(cè)量方式更具備優(yōu)勢(shì)［2］。LiDAR遙感通過激光脈沖信號(hào)，以點(diǎn)云數(shù)據(jù)形式感知森林三維場(chǎng)景，通過分析處理激光點(diǎn)云數(shù)據(jù)，可計(jì)算得到樹高［3］、胸徑［4］、冠幅［5］等森林結(jié)構(gòu)參數(shù)，進(jìn)而為蓄木量［6］、固碳量［7］、生物量［8］等宏觀森林生態(tài)參數(shù)信息反演提供數(shù)據(jù)基礎(chǔ)。激光雷達(dá)遙感成果在森林動(dòng)態(tài)監(jiān)測(cè)、林業(yè)資源精細(xì)化管理等任務(wù)中日趨重要［9］。

地基激光雷達(dá)（terrestrial laser scanning， TLS）遙感和機(jī)載激光雷達(dá)（airborne laser scanning， ALS）遙感是目前該領(lǐng)域任務(wù)中主流的森林三維結(jié)構(gòu)感知方式［10］。TLS以地面基站為平臺(tái)，常用于單木、樣地尺度測(cè)量［11］。TLS基站通過水平和豎直面內(nèi)的二軸旋轉(zhuǎn)來獲取關(guān)于目標(biāo)對(duì)象千萬級(jí)以上點(diǎn)云數(shù)據(jù)，能夠準(zhǔn)確捕捉到樹干等林下結(jié)構(gòu)［12］。ALS以有人駕駛飛機(jī)或無人機(jī)為載體，掃描覆蓋面廣，在區(qū)域尺度森林資源調(diào)查等任務(wù)中較為常見［13］。ALS激光脈沖穿透樹冠間隙，并收集多次回波信息，以此探測(cè)不同冠層和森林垂直結(jié)構(gòu)［14-15］。但是，由于數(shù)據(jù)采集過程中激光傳感器視角限制、激光脈沖能量衰減、樹冠遮擋等因素，這2種采集方式均存在單一源數(shù)據(jù)缺失、數(shù)據(jù)質(zhì)量不佳等問題［16］。因此，如何高效、準(zhǔn)確配準(zhǔn)機(jī)載-地基點(diǎn)云，對(duì)于多源點(diǎn)云融合、提高數(shù)據(jù)完整性以及觀測(cè)精度具有重要意義。

針對(duì)機(jī)載和地基激光點(diǎn)云匹配問題，最直接的方法是通過手工方式提取實(shí)地勘測(cè)過程中所布設(shè)的地面控制點(diǎn)［17］或者從樹頂［18］和地形［19］中提取特征，實(shí)現(xiàn)匹配。但匹配特征選取易受主觀影響，其結(jié)果具有一定偏向性，推廣性有限。近年來的研究更多集中于自動(dòng)化的無控匹配，主要應(yīng)用在針葉林等較為簡(jiǎn)單森林結(jié)構(gòu)場(chǎng)景中。Liu等［20］以冠層高度模型為基礎(chǔ)，運(yùn)用滑動(dòng)窗口搜索方式提取具備相似灰度值的點(diǎn)對(duì)作為特征對(duì)，以此求解旋轉(zhuǎn)矩陣，并運(yùn)用最近鄰點(diǎn)迭代（ICP）方法精化配準(zhǔn)結(jié)果，配準(zhǔn)結(jié)果距離差異為0.10～0.43m。但該方法在生成冠層高度模型中造成的信息損失會(huì)影響特征點(diǎn)對(duì)提取與最終匹配結(jié)果。Fekry等［21］和Polewski等［22］以ALS和TLS分割后的單木點(diǎn)云為基礎(chǔ)，提取單木平面質(zhì)心、樹頂?shù)忍卣鲗?duì)，進(jìn)而求解旋轉(zhuǎn)矩陣，配準(zhǔn)精度為0.3～0.7m。這類方法的問題在于需預(yù)先進(jìn)行點(diǎn)云分割，無疑大幅增加了計(jì)算成本。此外，Dai等［23］通過提取樹冠點(diǎn)云分布模型，運(yùn)用相關(guān)點(diǎn)漂移方法提取特征，實(shí)現(xiàn)粗配準(zhǔn)，再利用迭代最近點(diǎn)（iterative closest point，ICP）算法精化配準(zhǔn)結(jié)果，精度約為0.3m。

針對(duì)上述問題，提出一種精度更高且可應(yīng)用于多種森林復(fù)雜度結(jié)構(gòu)的自動(dòng)化空地點(diǎn)云融合匹配方法。該方法首先利用ALS和TLS軸對(duì)齊包圍盒中心，與主方向向量進(jìn)行初步定向，初始化相對(duì)位置，在此基礎(chǔ)上分別提取不同源樹冠層點(diǎn)云，并分別進(jìn)行體素化，將數(shù)據(jù)空間轉(zhuǎn)換為更加有序的數(shù)據(jù)形式，然后以TLS點(diǎn)云為模板窗口進(jìn)行滑動(dòng)，在ALS空間中尋找具備最大化匹配得分的位置，并返回匹配特征對(duì)，最后將問題轉(zhuǎn)換為最小化ALS和TLS點(diǎn)云空間差異，通過奇異值分解解算旋轉(zhuǎn)矩陣，并計(jì)算平移矩陣，最終實(shí)現(xiàn)機(jī)載與地基激光點(diǎn)云配準(zhǔn)。

1 機(jī)載-地基點(diǎn)云無控配準(zhǔn)方法

1.1 無控配準(zhǔn)流程

如圖1，所提出的機(jī)載-地基激光點(diǎn)云無控配準(zhǔn)方法主要包括：點(diǎn)云包圍盒定向、樹冠點(diǎn)云濾波、樹冠點(diǎn)云體素化、基于體素模板的特征對(duì)提取以及旋轉(zhuǎn)平移矩陣計(jì)算。

圖1 機(jī)載-地基無控匹配流程Fig.1 Flowchart to show mark-free ALS and TLS point clouds registration

1.2 基于點(diǎn)云包圍盒的位置初始化

該步驟目的是消除較大的平移旋轉(zhuǎn)差異，先對(duì)ALS和TLS點(diǎn)云構(gòu)建軸對(duì)齊包圍盒，分別確定各自包圍盒x、y、z主方向（圖2）。在此基礎(chǔ)上，對(duì)齊ALS和TLS點(diǎn)云包圍盒中心，減小平移差異，按式（1）計(jì)算ALS點(diǎn)云和TLS點(diǎn)云各主方向向量夾角。

圖2 包圍盒定向示意Fig.2 ALS and TLS position initialization with bounding box orientation

式中：A為單一分量的旋轉(zhuǎn)角度，其取值范圍為γ、β、α，分別為相對(duì)于x、y、z軸的旋轉(zhuǎn)角度；Valsi和Vtlsi分別為ALS點(diǎn)云和TLS點(diǎn)云單一主方向分量向量。

當(dāng)角度計(jì)算完成后，分別按照式（2）、（3）、（4）分別計(jì)算相對(duì)于x、y、z軸的旋轉(zhuǎn)矩陣Rx(γ)、Ry(β)、Rz(α)，最后依次將各主方向旋轉(zhuǎn)矩陣應(yīng)用于原始TLS點(diǎn)云上，減小其與目標(biāo)ALS點(diǎn)云的差異，實(shí)現(xiàn)位置初始化。

1.3 基于高度與密度分布的樹冠點(diǎn)云濾波

森林包括不同層級(jí)樹冠。在不同高度，點(diǎn)云頻數(shù)也呈現(xiàn)出一定規(guī)律性，同時(shí)高大樹冠比低矮樹冠形態(tài)更為清晰和穩(wěn)定，因此從中提取的特征更為清晰與穩(wěn)定。由圖3a樹高與點(diǎn)云頻數(shù)分布折線圖可知，示例森林樣區(qū)（圖3b）主要包括2級(jí)樹冠層，分別對(duì)應(yīng)2個(gè)波峰。通過區(qū)間定位波谷位置（局部極小值點(diǎn)），得到濾波閾值點(diǎn)，進(jìn)而提取樹冠上層點(diǎn)云，作為后續(xù)特征提取的數(shù)據(jù)基礎(chǔ)。

圖3 基于樹高與頻數(shù)分布的樹冠點(diǎn)云濾波Fig.3 Canopy point clouds filtering with height and density

1.4 樹冠點(diǎn)云體素化

分別對(duì)ALS和TLS中所提取的冠層點(diǎn)云進(jìn)行體素化，劃分體素格網(wǎng)。通過格網(wǎng)表達(dá)，將無序點(diǎn)云轉(zhuǎn)換為規(guī)則化的數(shù)據(jù)結(jié)構(gòu)（圖4），為后續(xù)模板匹配提供數(shù)據(jù)基礎(chǔ)，同時(shí)也有利于降低數(shù)據(jù)量，提升算法整體效率。

1.5 基于體素模板的特征對(duì)提取

針對(duì)森林點(diǎn)云數(shù)據(jù)難以有效提取不變特征的問題，提出一種區(qū)域模板匹配的方法。以最大化匹配得分（Stotal）為目標(biāo)函數(shù)，尋找最佳搜索位置，從而獲取特征對(duì)。該模塊主要包括體素過濾、模板搜索、最大化得分以及特征對(duì)提取四部分（圖5）。

圖5 模板匹配流程Fig.5 Workflow for template matching process

在ALS和TLS體素質(zhì)心集中，相同平面位置可能存在多個(gè)不同z值冗余體素中心，因此進(jìn)一步通過限定最大高度值提取樹冠上表層體素集。

在模板匹配過程中，以整體面積更小的TLS體素空間為模板，并將其坐標(biāo)通過式（5）與ALS體素空間對(duì)齊，并通過x和y方向上的增量實(shí)現(xiàn)模板滑動(dòng)。

式中：Px和Py分別表示體素空間對(duì)齊后的x和y方向的坐標(biāo)值。

在每個(gè)搜索位置，對(duì)ALS和TLS體素中心集整體建立KDTree［24］，對(duì)KDTree中每一個(gè)TLS體素質(zhì)心點(diǎn)檢索其最近鄰的ALS體素質(zhì)心點(diǎn)，若存在最近鄰點(diǎn)，且兩者間Mahalanobis距離［25］小于一個(gè)體素長(zhǎng)度（ε），則認(rèn)定其為潛在匹配點(diǎn)對(duì)，增加一個(gè)匹配得分，見式（6）：

式中：Si為單一質(zhì)心點(diǎn)的匹配得分；Ptlsi、Palsi分別表示TLS和ALS質(zhì)心集坐標(biāo)分量值。因此，單一搜索位置匹配總得分則為潛在點(diǎn)對(duì)的匹配得分總和，見式（7）：

式中：Stotal為當(dāng)前位置匹配總得分。

完成所有位置搜索后返回具有最大匹配得分位置的TLS和ALS特征點(diǎn)對(duì)，作為旋轉(zhuǎn)平移矩陣求解步驟中的輸入。

1.6 旋轉(zhuǎn)平移矩陣求解

在森林場(chǎng)景下，ALS和TLS點(diǎn)云匹配為剛體變換，不考慮形變和縮放問題，因此，該問題轉(zhuǎn)換為最小化TLS和ALS點(diǎn)云差異，見式（8）：

式中：R、t分別為旋轉(zhuǎn)矩陣和平移矩陣。

分別對(duì)ALS和TLS點(diǎn)云去中心化，并構(gòu)建協(xié)方差矩陣，見式（10）：

式中：μals和μtls分別為ALS與TLS點(diǎn)云質(zhì)心坐標(biāo)向量，在此基礎(chǔ)上通過SVD奇異值分解方法獲取U、S、V矩陣（式（11）），通過式（12）解算旋轉(zhuǎn)矩陣R，最后將R矩陣回代式（8），利用式（13）計(jì)算得到平移矩陣t。

最終將計(jì)算得到的旋轉(zhuǎn)矩陣R和平移矩陣t應(yīng)用于TLS點(diǎn)云，實(shí)現(xiàn)機(jī)載-地基點(diǎn)云配準(zhǔn)。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

采用4種不同森林結(jié)構(gòu)復(fù)雜度的空地點(diǎn)云數(shù)據(jù)，對(duì)所提出的無控配準(zhǔn)方法進(jìn)行驗(yàn)證與比較（表1）。

表1 4種不同森林結(jié)構(gòu)的機(jī)載激光雷達(dá)（ALS）與地基激光雷達(dá)（TLS）點(diǎn)云數(shù)據(jù)Tab.1 Four sets of ALS and TLS point clouds with different forest structures

樣地1為馬來西亞婆羅洲某自然保護(hù)區(qū)熱帶雨林，機(jī)載點(diǎn)云區(qū)域大小約125m×125m，地基點(diǎn)云范圍約為70m×70m，區(qū)域中主要包括石南樹、龍果樹等樹種，最大樹高約達(dá)158m。樣地2為德國(guó)卡爾斯魯厄地區(qū)闊葉、針葉混合林，點(diǎn)云數(shù)據(jù)來源于文獻(xiàn)［19］，機(jī)載點(diǎn)云區(qū)域大小約為145m×145m，地基點(diǎn)云范圍約為200m×200m，其中主要包括樟子松、橡木等樹種，最大樹高約為62m。樣地3和樣地4為分別為荷蘭Veluwe地區(qū)的闊葉林和針葉林，點(diǎn)云數(shù)據(jù)均來源于文獻(xiàn)［17］，機(jī)載點(diǎn)云區(qū)域大小分別約為95m×95m和60m×60m，地基點(diǎn)云范圍分別約為102m×102m和51m×51m，主要包括樺木、云杉等樹種，最大樹高分別約45m和40m。4組數(shù)據(jù)按森林結(jié)構(gòu)由復(fù)雜到簡(jiǎn)單排序?yàn)椋簶拥?、樣地2、樣地3、樣地4。4組樣地ALS點(diǎn)云量均為百萬至千萬級(jí)，其密度為139～3 865個(gè)·m-2，TLS點(diǎn)云數(shù)據(jù)量為千萬至億級(jí)，其密度為2 375～20 381個(gè)·m-2。其他參數(shù)見表2。

表2 機(jī)載激光雷達(dá)（ALS）點(diǎn)云與地基激光雷達(dá)（TLS）點(diǎn)云數(shù)據(jù)指標(biāo)詳情Tab.2 Details for ALS and TLS point clouds

2.2 機(jī)載-地基點(diǎn)云配準(zhǔn)精度評(píng)價(jià)指標(biāo)

為避免倒樹等可能造成無法匹配最近點(diǎn)的情況，通過人工采樣ALS和TLS點(diǎn)云中樹冠與樹干位置的平均距離差異來描述機(jī)載-地基點(diǎn)云匹配成果質(zhì)量。即構(gòu)建KDtree建立空間索引，以轉(zhuǎn)換后點(diǎn)云為基準(zhǔn)，按照式（14）計(jì)算得到其到目標(biāo)點(diǎn)云中最近點(diǎn)的平均距離差異。

式中：D為轉(zhuǎn)換后的TLS點(diǎn)云與目標(biāo)ALS點(diǎn)云的平均距離差異；和分別代表轉(zhuǎn)換后的TLS點(diǎn)云坐標(biāo)向量和目標(biāo)ALS點(diǎn)云坐標(biāo)向量。

2.3 實(shí)驗(yàn)參數(shù)選擇

基于Python3.7.12軟件實(shí)現(xiàn)所提出的方法，分別運(yùn)用Open3D構(gòu)建KDTree與檢索最近鄰點(diǎn)，以Numpy庫實(shí)現(xiàn)矩陣運(yùn)算與SVD分解，以multiprocess庫并行加速。實(shí)驗(yàn)硬件環(huán)境為64G內(nèi)存、英特爾2.6GHz CPU。

在基于樹高-頻數(shù)的樹冠點(diǎn)云濾波環(huán)節(jié)中，通過區(qū)間對(duì)比可知，檢測(cè)到的4塊林區(qū)上層樹冠點(diǎn)云高度分別約為120m、45m、32m、34m。此外，所選擇的體素格網(wǎng)尺寸均定義為1。

對(duì)比實(shí)驗(yàn)選擇了人工配準(zhǔn)、最近鄰迭代法［26］（iterative closest point， ICP）以及上述二者組合的方法。人工配準(zhǔn)方法是通過手動(dòng)選取測(cè)區(qū)周邊和中心區(qū)域9組幾何特征顯著點(diǎn)對(duì)作為匹配特征，進(jìn)而求解旋轉(zhuǎn)平移矩陣，實(shí)現(xiàn)點(diǎn)云匹配。在樣地1中，特征點(diǎn)為樹主干與枝干交叉點(diǎn)。在樣地2、3和4中，樹頂點(diǎn)作為特征點(diǎn)。單獨(dú)使用ICP方法時(shí)，由于其對(duì)初始位置要求較高，因此手動(dòng)將TLS點(diǎn)云移動(dòng)至與對(duì)應(yīng)ALS點(diǎn)云有重疊區(qū)域的位置，并且迭代參數(shù)選擇為1 500，該參數(shù)在組合方法中選擇為1 000。

2.4 機(jī)載-地基點(diǎn)云配準(zhǔn)結(jié)果與分析

對(duì)4組方法在4組樣地中機(jī)載、地基點(diǎn)云匹配結(jié)果按式（15）計(jì)算平均距離差異。

表3表明，本文所提出的方法在4組不同森林復(fù)雜度場(chǎng)景中均表現(xiàn)出較低的平均距離差異，分別為0.245m、0.238m、0.184m和0.020m。同時(shí)，人工方法也具有較高的精度，最大距離差異低于0.37m。相比下，基于位置初始化的單一ICP方法在這4組結(jié)果中均出現(xiàn)較為顯著的不一致性（圖6），其最小數(shù)值超過0.48m，最大差異達(dá)18m。并且在組合方法中，ICP算法在多數(shù)情況下對(duì)最終結(jié)果造成負(fù)面影響。除了樣地2數(shù)據(jù)外，組合方法匹配在平均距離差異上相較于單一人工配準(zhǔn)方法分別上升了0.378m、0.082m和0.325m。

表3 機(jī)載點(diǎn)云與地基點(diǎn)云匹配結(jié)果平均距離差異Tab.3 Average spatial difference of ALS-TLS registration result （單位：m）

圖6 機(jī)載-地基點(diǎn)云配準(zhǔn)結(jié)果細(xì)節(jié)Fig.6 Details of ALS and TLS point clouds registration results

此外，本文所提出方法在應(yīng)對(duì)不同森林場(chǎng)景時(shí)也表現(xiàn)出優(yōu)良的穩(wěn)定性。其中，在結(jié)構(gòu)最簡(jiǎn)單的針葉林中效果最好，其平均距離差異僅為0.020m。在其他更為復(fù)雜的混合林（樣地3）、闊葉林（樣地2）以及熱帶雨林（樣地1）中，精度浮動(dòng)在0.061m內(nèi)。而人工配準(zhǔn)與組合方法精度出現(xiàn)較大波動(dòng)，上述3個(gè)樣地中浮動(dòng)分別在0.142m和0.308m。

對(duì)于ICP算法而言，盡管在執(zhí)行算法前通過人工方式將ALS和TLS點(diǎn)云初始化為相似位置，但森林環(huán)境數(shù)據(jù)更加雜亂，特征對(duì)應(yīng)關(guān)系更加模糊，若無穩(wěn)定不變性的輸入特征，很可能導(dǎo)致算法失效。這也是單一ICP算法失效以及其組合方法在大多數(shù)樣地實(shí)驗(yàn)中出現(xiàn)精度下降的主要原因。

人工配準(zhǔn)是人為選取樹頂和樹干交叉作為特征點(diǎn)，具備較好的穩(wěn)定性，因而配準(zhǔn)精度更高。然而樹頂特征是不穩(wěn)定的，通過圖6可看到，樣地3和4的匹配結(jié)果在冠層以下部分依然出現(xiàn)明顯偏差，而樹干交叉特征相對(duì)穩(wěn)定，其配準(zhǔn)結(jié)果與本文所提出方法基本一致。

相比而言，本文方法更優(yōu)的原因主要有兩方面：第一，樹冠點(diǎn)云體素化解決了原數(shù)據(jù)空間的點(diǎn)云無序性與散亂性問題，體素分布更加規(guī)律，有助于后續(xù)特征提取。第二，在特征對(duì)選擇中通過定義體素模板并通過網(wǎng)格滑動(dòng)方式尋找在最優(yōu)匹配位置下的特征點(diǎn)對(duì)，該過程本質(zhì)是基于局部采樣點(diǎn)集提取具備最佳樹冠層體素拓?fù)浣Y(jié)構(gòu)相似度的配準(zhǔn)特征對(duì)，而非簡(jiǎn)單地基于點(diǎn)特征或基于面特征的配準(zhǔn)，因此算法精度和穩(wěn)定性更優(yōu)。

3 結(jié)論

針對(duì)森林場(chǎng)景，提出了一種基于樹冠體素特征的機(jī)載-地基激光點(diǎn)云無控配準(zhǔn)方法。首先初始化ALS和TLS點(diǎn)云包圍盒方位來消除顯著方位差異，然后以樹冠層點(diǎn)云體素質(zhì)心為基礎(chǔ)，通過體素模板搜索過程計(jì)算最大匹配得分位置提取匹配特征對(duì)，并以此計(jì)算旋轉(zhuǎn)平移矩陣，最終實(shí)現(xiàn)機(jī)載-地基點(diǎn)云配準(zhǔn)?；卺樔~林、混合林、闊葉林和熱帶雨林4組不同結(jié)構(gòu)復(fù)雜度場(chǎng)景的機(jī)載-地基點(diǎn)云數(shù)據(jù)，與ICP、人工配準(zhǔn)、人工配準(zhǔn)與ICP結(jié)合這3種方法進(jìn)行對(duì)比實(shí)驗(yàn)，表明所本文方法在精確度與魯棒性方面均最優(yōu)，平均距離差異分別為0.245m、0.238m、0.184m、0.020m。配準(zhǔn)結(jié)果極大提高了森林場(chǎng)景點(diǎn)云數(shù)據(jù)完整性，有助于多源點(diǎn)云融合、森林結(jié)構(gòu)生態(tài)參數(shù)反演等工作。

作者貢獻(xiàn)聲明：

林 怡：研究方向確定、算法指導(dǎo)、論文修改。

曹宇杰：算法實(shí)現(xiàn)、數(shù)據(jù)采集、數(shù)據(jù)分析、論文撰寫與修改。