王楚虹，劉浩然，鐘 浩，林文樹

（東北林業(yè)大學(xué) 工程技術(shù)學(xué)院，黑龍江 哈爾濱 150040）

激光雷達(dá)作為一種高精度的主動(dòng)遙感技術(shù)，在林業(yè)應(yīng)用領(lǐng)域具有傳統(tǒng)光學(xué)遙感無可替代的優(yōu)勢[1-3]。按照平臺(tái)的不同，激光雷達(dá)主要分為星載激光雷達(dá)、機(jī)載激光雷達(dá)、車載激光雷達(dá)及地基激光雷達(dá)等[4]。單一激光雷達(dá)技術(shù)平臺(tái)測量森林結(jié)構(gòu)參數(shù)時(shí)具有一定局限性，如地面移動(dòng)激光雷達(dá)技術(shù)可快速準(zhǔn)確地獲得樹干點(diǎn)云數(shù)據(jù)信息，但對于樹冠點(diǎn)云信息的獲得存在一定缺失，而無人機(jī)激光雷達(dá)技術(shù)雖然可以有效獲取林分冠層信息，但對于樹干點(diǎn)云信息卻無法做到精確獲取。因此，要利用激光雷達(dá)技術(shù)同時(shí)獲取森林的水平和垂直結(jié)構(gòu)信息，就需要進(jìn)行多源遙感平臺(tái)數(shù)據(jù)的融合來彌補(bǔ)單一激光雷達(dá)點(diǎn)云數(shù)據(jù)在林業(yè)中應(yīng)用的技術(shù)瓶頸。

目前已經(jīng)有相關(guān)學(xué)者將遙感影像作為激光雷達(dá)數(shù)據(jù)的補(bǔ)充應(yīng)用于森林參數(shù)提取及三維模型建立中，其中部分研究是將遙感影像與激光雷達(dá)點(diǎn)云數(shù)據(jù)配準(zhǔn)后進(jìn)行聯(lián)合應(yīng)用。張吳明等[5]基于人為設(shè)置控制點(diǎn)的方法將無人機(jī)影像與地基激光雷達(dá)數(shù)據(jù)進(jìn)行配準(zhǔn)，并轉(zhuǎn)換為在同一坐標(biāo)系下的點(diǎn)云，以此提取樹高。曹明蘭等[6]將無人機(jī)傾斜影像生成密集點(diǎn)云后，使用迭代最近點(diǎn)（ICP,Iterative Closest Point）算法將密集點(diǎn)云與背包式三維激光掃描系統(tǒng)獲取的點(diǎn)云配準(zhǔn)，并構(gòu)建三維模型，使2 種技術(shù)實(shí)現(xiàn)了優(yōu)勢互補(bǔ)。李丹[7]采用聯(lián)合地基激光雷達(dá)與機(jī)載成像激光雷達(dá)的方法對森林參數(shù)進(jìn)行提取，并以地基提取的 DEM 為基準(zhǔn)，對機(jī)載影像解算的 DSM 進(jìn)行地理坐標(biāo)配準(zhǔn)。激光雷達(dá)數(shù)據(jù)聯(lián)合無需配準(zhǔn)的同參考坐標(biāo)系下的平臺(tái)影像數(shù)據(jù)也可以有效提高森林參數(shù)獲取的精度。駱社周等[8]處理星載激光雷達(dá)GLAS 數(shù)據(jù)建立反演森林葉面積指數(shù)模型，并利用神經(jīng)網(wǎng)絡(luò)融合星載激光雷達(dá)GLAS 與TM 光學(xué)遙感數(shù)據(jù)，以此實(shí)行區(qū)域連續(xù)高精度森林葉面積指數(shù)的反演，為生態(tài)環(huán)境研究精準(zhǔn)輸入?yún)?shù)提供了新思路。云增鑫等[9]利用航空激光雷達(dá)數(shù)據(jù)對高光譜數(shù)據(jù)進(jìn)行林下植被信息的剔除，研究發(fā)現(xiàn)處理后的高光譜遙感影像能夠有效改善森林冠層有效葉面積指數(shù)的估算精度。

不同平臺(tái)的激光雷達(dá)點(diǎn)云數(shù)據(jù)也可以相互補(bǔ)充，因而有學(xué)者針對無標(biāo)識(shí)不同平臺(tái)激光雷達(dá)點(diǎn)云數(shù)據(jù)的融合進(jìn)行了研究。2017年，Paris 等[10]利用基于點(diǎn)云數(shù)據(jù)柵格化的圖像配準(zhǔn)方法實(shí)現(xiàn)了地基和機(jī)載點(diǎn)云數(shù)據(jù)融合，但該方法適用于樹木形狀不規(guī)則的且較為開闊的林區(qū)，算法復(fù)雜度高。Polewski 等[11]使用模擬退火算法對二維樹的位置點(diǎn)之間相對距離進(jìn)行比較，完成地面移動(dòng)激光雷達(dá)點(diǎn)云數(shù)據(jù)與空中激光雷達(dá)點(diǎn)云數(shù)據(jù)的配準(zhǔn)，精度較好。Guan 等[12]提出了一種多平臺(tái)激光雷達(dá)數(shù)據(jù)自動(dòng)配準(zhǔn)框架，即建立不規(guī)則三角網(wǎng)并通過三角形相似度投票策略來完成粗配準(zhǔn)，為無人機(jī)激光雷達(dá)與地面移動(dòng)激光雷達(dá)森林點(diǎn)云數(shù)據(jù)的融合提出了有效解決方法。Dai 等[13]通過樹冠密度分析來融合機(jī)載激光掃描和地面激光掃描的森林點(diǎn)云，利用均值偏移法提取基于模型的關(guān)鍵點(diǎn)，并利用最大似然估計(jì)對關(guān)鍵點(diǎn)進(jìn)行排序，該方法避免了對森林點(diǎn)云的幾何描述。

隨著以地面手持式、背包式及無人機(jī)為平臺(tái)的輕小型激光雷達(dá)掃描系統(tǒng)的快速發(fā)展[14]，其在林業(yè)中的應(yīng)用也越來越廣泛。目前，基于地面移動(dòng)與無人機(jī)平臺(tái)的林地點(diǎn)云數(shù)據(jù)配準(zhǔn)研究尚在起步階段，大多數(shù)地空平臺(tái)的激光雷達(dá)點(diǎn)云數(shù)據(jù)融合研究多應(yīng)用于城市建筑方面[15-17]，而在林業(yè)應(yīng)用方面的研究較少。另外，林業(yè)激光雷達(dá)點(diǎn)云數(shù)據(jù)配準(zhǔn)中，部分研究通過利用周圍建筑物或設(shè)立明顯標(biāo)志物完成[18]，選擇的森林樣地中林分也比較稀疏。因此面對復(fù)雜的林分環(huán)境，在不易設(shè)立地面標(biāo)志物的情況下，如何將地面移動(dòng)激光雷達(dá)數(shù)據(jù)與無人機(jī)激光雷達(dá)數(shù)據(jù)的精準(zhǔn)與高效融合仍有待進(jìn)一步研究?；诖?，本研究提出一種基于Delaunay 三角網(wǎng)的手持移動(dòng)激光雷達(dá)數(shù)據(jù)與無人機(jī)激光雷達(dá)數(shù)據(jù)的無標(biāo)識(shí)融合方法，旨在以較高的效率及精準(zhǔn)度完成2 種平臺(tái)的配準(zhǔn)，為提取完整的森林結(jié)構(gòu)參數(shù)奠定基礎(chǔ)，從而推動(dòng)激光雷達(dá)技術(shù)在林業(yè)上的廣泛應(yīng)用。

1 研究樣地與激光雷達(dá)數(shù)據(jù)獲取

1.1 樣地概況

研究樣地位于哈爾濱城市林業(yè)示范基地，地理坐標(biāo)為45°43′10″N，126°37′15″E。在樟子松人工林與蒙古櫟人工林中分別選取2 塊15 m×15 m的區(qū)域作為研究樣地，其中，樟子松樣地樹木分布較松散，林中灌木較多；蒙古櫟樣地樹木分布較密集，林木生長繁茂，郁閉度較高。選取的樣地與基地內(nèi)道路相鄰，便于在道路上設(shè)置用于評價(jià)配準(zhǔn)精度的標(biāo)志物。為了使無人機(jī)激光雷達(dá)能夠清晰地掃描到標(biāo)志物，利用放置于三腳架上的白色標(biāo)靶球作為標(biāo)志物，并將三腳架放置在無樹冠遮擋的位置，確保標(biāo)志物被手持移動(dòng)激光雷達(dá)與無人機(jī)激光雷達(dá)2 個(gè)平臺(tái)完整掃描。

1.2 激光雷達(dá)數(shù)據(jù)獲取

2020年10月，分別使用手持移動(dòng)激光雷達(dá)系統(tǒng)和無人機(jī)激光雷達(dá)系統(tǒng)進(jìn)行數(shù)據(jù)獲取。手持移動(dòng)激光雷達(dá)掃描系統(tǒng)選取GeoSLAM 手持移動(dòng)SLAM 激光雷達(dá)ZEB-HORIZON，其測程為100 m，采集速度為每秒30 萬點(diǎn)，獲取點(diǎn)云的坐標(biāo)系為相對坐標(biāo)系。設(shè)置間隔3 m 左右的蛇形的行走路線，采集過程平穩(wěn)緩慢，林木密集復(fù)雜處以及轉(zhuǎn)彎處進(jìn)行停留采集數(shù)據(jù)，確保樹木數(shù)據(jù)的完整性。數(shù)據(jù)采集完成后，使用GEOSLAM-Hub 對數(shù)據(jù)進(jìn)行SLAM 拼接，生成LAS 格式文件。

無人機(jī)激光雷達(dá)系統(tǒng)使用蜂鳥Genius 微型無人機(jī)LiDAR 系統(tǒng)，最大測距優(yōu)于250 m，最大測量速度64 萬點(diǎn)/s，實(shí)際作業(yè)飛行密度點(diǎn)大于200點(diǎn)/m2。實(shí)際飛行高度100 m，選取東北林業(yè)大學(xué)哈爾濱城市林業(yè)示范基地300 m × 160 m 的區(qū)域進(jìn)行無人機(jī)激光雷達(dá)數(shù)據(jù)獲取，確保研究樣地包含在內(nèi)。通過POS 數(shù)據(jù)解算，完成對無人機(jī)激光雷達(dá)點(diǎn)云的整合，點(diǎn)云坐標(biāo)系為WGS-84 大地坐標(biāo)系。

2 HMLS 和UAV-LiDAR 數(shù)據(jù)融合方法

在融合前，需要對手持移動(dòng)激光雷達(dá)數(shù)據(jù)和無人機(jī)激光雷達(dá)數(shù)據(jù)進(jìn)行預(yù)處理。首先，將2 種平臺(tái)得到的點(diǎn)云數(shù)據(jù)進(jìn)行濾波去噪。然后將無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)根據(jù)分離出的地面點(diǎn)歸一化處理；將手持移動(dòng)激光雷達(dá)點(diǎn)云數(shù)據(jù)進(jìn)行均勻抽稀處理，保證其在不丟失樹木自身特征的前提下與無人機(jī)激光雷達(dá)點(diǎn)云的密度基本一致。最后將處理后的手持移動(dòng)激光雷達(dá)點(diǎn)云數(shù)據(jù)和無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)進(jìn)行單木分割，分別提取樹干的中心點(diǎn)與樹冠的幾何中心點(diǎn)作為樹木的位置點(diǎn)。

融合手持移動(dòng)激光雷達(dá)數(shù)據(jù)與無人機(jī)激光雷達(dá)數(shù)據(jù)的方法分為粗配準(zhǔn)和精配準(zhǔn)2 個(gè)步驟。將預(yù)處理后的2 個(gè)不同平臺(tái)點(diǎn)云數(shù)據(jù)基于樹木位置點(diǎn)構(gòu)建Delaunay 三角網(wǎng)進(jìn)行粗配準(zhǔn)并利用模擬退火算法優(yōu)化配準(zhǔn)過程。選取一定高度區(qū)域，將粗配準(zhǔn)的手持移動(dòng)激光雷達(dá)點(diǎn)云數(shù)據(jù)基于無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)運(yùn)用ICP 算法進(jìn)一步精配準(zhǔn)調(diào)整。選取樟子松樣地和蒙古櫟樣地對融合方法進(jìn)行驗(yàn)證，下文中的點(diǎn)云融合過程以樟子松樣地為例，點(diǎn)云融合平臺(tái)基于Python 編譯。具體流程如圖1所示。

圖1 點(diǎn)云配準(zhǔn)流程Fig.1 Point cloud registration flow chart

2.1 Delaunay 三角網(wǎng)的構(gòu)建

構(gòu)建Delaunay 三角網(wǎng)是粗配準(zhǔn)的第一步，也是提高粗配準(zhǔn)效率的關(guān)鍵步驟。Delaunay 三角網(wǎng)是一系列相連的、不重疊的三角形的集合，該集合內(nèi)任意一個(gè)三角形的外接圓中不包含平面內(nèi)其他點(diǎn)[19]。對于同一個(gè)點(diǎn)集而言，從任意一點(diǎn)構(gòu)建Delaunay 三角網(wǎng)，所得結(jié)果都是一致的。本研究分別將預(yù)處理提取的手持移動(dòng)激光雷達(dá)和無人機(jī)激光雷達(dá)數(shù)據(jù)樹木位置的三維坐標(biāo)投影于平面，形成2組二維平面位置點(diǎn)。由于樹木位置是固定的，又因相同離散點(diǎn)構(gòu)建的Delaunay 三角網(wǎng)具有唯一性，故使用2 組樹木的二維位置點(diǎn)集合分別構(gòu)建的2 個(gè)Delaunay 三角網(wǎng)形狀應(yīng)基本相同。

本研究選擇效率高、使用最為廣泛的逐點(diǎn)插入算法構(gòu)建Delaunay 三角網(wǎng)[20]。使用逐點(diǎn)插入算法構(gòu)建基于樹木位置點(diǎn)的Delaunay 三角網(wǎng)的基本步驟如下：1）根據(jù)樹木位置點(diǎn)的最大分布來求得隨機(jī)一個(gè)超級三角形，并將其加入臨時(shí)三角形集合中；2）依次插入樹木位置點(diǎn)，找出臨時(shí)三角形集合中外接圓包含插入點(diǎn)的三角形T，刪除三角形T 的公共邊，并將插入的樹木位置點(diǎn)與T 的頂點(diǎn)連接，完成該點(diǎn)的插入，將形成的三角形放入臨時(shí)三角形集合；3）重復(fù)2）步驟，直到所有樹木位置點(diǎn)插入完成。

2.2 對應(yīng)點(diǎn)的確定

理想情況下，2 組樹木位置點(diǎn)構(gòu)建的Delaunay三角網(wǎng)形狀應(yīng)完全相同，但是樹木郁閉度高等因素會(huì)影響無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)單木分割的精度，導(dǎo)致例如樹木位置點(diǎn)未被記錄或記錄錯(cuò)誤等偏差出現(xiàn)。另外，手持移動(dòng)激光雷達(dá)點(diǎn)云單木位置是樹干中心，而無人機(jī)激光雷達(dá)點(diǎn)云單木位置是樹冠中心，在樹干傾斜的情況下2 點(diǎn)在同一平面的投影點(diǎn)不能重合。因此在2 個(gè)由樹木位置點(diǎn)構(gòu)造的三角網(wǎng)仍可能會(huì)出現(xiàn)較大差別的情況下，需通過設(shè)定限制確保2 組樹木位置點(diǎn)之間建立準(zhǔn)確的對應(yīng)關(guān)系。

以Delaunay 三角網(wǎng)為基礎(chǔ)，通過以下限制搜索2 種平臺(tái)所得到的樹木位置點(diǎn)。1）搜索所有線段，設(shè)定合理閾值，記錄差值在閾值內(nèi)的線段及其對應(yīng)線段，兩兩構(gòu)成線段組。2）將線段組整合，記錄手持移動(dòng)激光雷達(dá)點(diǎn)云中同一樹木位置點(diǎn)出發(fā)的線段所對應(yīng)的無人機(jī)激光雷達(dá)點(diǎn)云線段。3）尋找對應(yīng)線段相交的公共點(diǎn)，該點(diǎn)即為手持移動(dòng)激光雷達(dá)點(diǎn)云中樹木位置點(diǎn)的對應(yīng)點(diǎn)。在搜索對應(yīng)點(diǎn)的過程中，設(shè)置線段差閾值為可調(diào)節(jié)的參數(shù)。當(dāng)線段長度差的閾值設(shè)置極小時(shí)，搜索過程中該方法對于線段與其對應(yīng)線段的長度偏差的容忍度極低。由于手持移動(dòng)激光雷達(dá)點(diǎn)云樹木位置點(diǎn)與無人機(jī)激光雷達(dá)點(diǎn)云樹木位置點(diǎn)并非全部垂直對應(yīng)，當(dāng)線段差閾值設(shè)置極小時(shí)，容易出現(xiàn)無法搜索到足夠?qū)?yīng)點(diǎn)的情況。反之，將線段長度差的閾值設(shè)置極大，表明搜索過程中容忍度極高，這會(huì)大大增加搜索進(jìn)程時(shí)間，使得搜索效率降低，同時(shí)有可能出現(xiàn)樹木位置點(diǎn)對應(yīng)錯(cuò)誤的問題，導(dǎo)致配準(zhǔn)失敗。

為了提高粗配準(zhǔn)的效率，本研究利用模擬退火算法優(yōu)化搜索無人機(jī)激光雷達(dá)樹木位置點(diǎn)對應(yīng)的手持移動(dòng)激光雷達(dá)樹木位置點(diǎn)這一過程。模擬退火算法是一個(gè)搜索全局最優(yōu)解的算法，是尋找溫度最低時(shí)概率分布中具有最大概率的狀態(tài)，并在搜索過程中設(shè)置了一個(gè)變化且最終趨于零的概率突跳性，避免搜索陷入局部最優(yōu)解[21]?；谀M退火優(yōu)化算法的具體配準(zhǔn)過程如下：1）設(shè)定較高的初始溫度T，隨機(jī)選取可能的無人機(jī)激光雷達(dá)點(diǎn)云樹木位置點(diǎn)與其手持移動(dòng)激光雷達(dá)點(diǎn)云對應(yīng)點(diǎn)作為初始模型A0，計(jì)算其相應(yīng)的目標(biāo)函數(shù)值S(A0)，目標(biāo)函數(shù)見式（1）；2）對當(dāng)前模型進(jìn)行擾動(dòng)產(chǎn)生一個(gè)新的可能的無人機(jī)激光雷達(dá)點(diǎn)云樹木位置點(diǎn)與其手持移動(dòng)激光雷達(dá)點(diǎn)云對應(yīng)點(diǎn)作為模型A，計(jì)算相應(yīng)的目標(biāo)函數(shù)值S(A)，得到ΔS=S(A)-S(A0)；3）若ΔS＜0，則新模型A被接受；若ΔS＞0，則新模型A按概率P=exp(-ΔS/T)進(jìn)行接受。當(dāng)模型被接受時(shí)，A0=A；4）溫度T下，多次進(jìn)行擾動(dòng)和接受；5）逐漸降低溫度T；6）重復(fù)步驟，至結(jié)果收斂或者到達(dá)迭代次數(shù)，迭代次數(shù)這一參數(shù)可在配準(zhǔn)過程中根據(jù)點(diǎn)云數(shù)據(jù)的數(shù)據(jù)量進(jìn)行相應(yīng)調(diào)節(jié)。模擬退火算法的目標(biāo)函數(shù)設(shè)為下式。

式（1）中，(x,y)表示無人機(jī)激光雷達(dá)點(diǎn)云樹木位置目標(biāo)點(diǎn)在平面的投影點(diǎn)，(xi,yi)為Delaunay 三角網(wǎng)中(x,y)的相鄰點(diǎn)，(xi,yi)與(x,y)兩點(diǎn)構(gòu)成目標(biāo)線段。(X,Y)為疑似手持移動(dòng)激光雷達(dá)點(diǎn)云樹木位置對應(yīng)點(diǎn)在平面的投影點(diǎn)，(X,Y)與其相鄰點(diǎn)(Xi,Yi)的距離為與目標(biāo)線段長度相似的對應(yīng)線段。通過計(jì)算，確定(X,Y)是否為(x,y)真正的對應(yīng)點(diǎn)。尋找所有目標(biāo)點(diǎn)的對應(yīng)點(diǎn)，使用模擬退火算法找到最優(yōu)對應(yīng)點(diǎn)集合，舍棄一些無法確定的點(diǎn)對。

2.3 旋轉(zhuǎn)平移矩陣

利用模擬退火算法搜索得到手持移動(dòng)激光雷達(dá)樹木位置點(diǎn)與無人機(jī)激光雷達(dá)樹木位置點(diǎn)對應(yīng)的點(diǎn)集合，利用奇異值分解法根據(jù)對應(yīng)點(diǎn)對的集合求得旋轉(zhuǎn)矩陣R與平移矩陣t[22]。

2.4 迭代最近點(diǎn)算法（ICP）

使用ICP 算法對粗配準(zhǔn)后的手持移動(dòng)激光雷達(dá)點(diǎn)云數(shù)據(jù)與歸一化的無人機(jī)激光雷達(dá)數(shù)據(jù)進(jìn)行精配準(zhǔn)。ICP 算法是基于最小二乘法的最優(yōu)匹配，通過將點(diǎn)與點(diǎn)不斷進(jìn)行匹配后，對點(diǎn)云旋轉(zhuǎn)平移，直到找到最優(yōu)變換矩陣，使目標(biāo)點(diǎn)云與待配準(zhǔn)點(diǎn)云重合。ICP 算法需要2 組點(diǎn)云形成一致或者包含的關(guān)系，才能確保配準(zhǔn)正常運(yùn)行[23]。手持移動(dòng)激光雷達(dá)點(diǎn)云數(shù)據(jù)包括樹干以及部分樹冠信息，而無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)包括樹冠、部分樹干及少量地面點(diǎn)，由此可知2 組點(diǎn)云數(shù)據(jù)中，重合的點(diǎn)云大部分屬于樹干上部及樹冠中下部區(qū)域。為保證ICP 算法能夠獲取理想配準(zhǔn)結(jié)果，選取一定高度區(qū)域的手持移動(dòng)激光雷達(dá)點(diǎn)云與無人機(jī)激光雷達(dá)點(diǎn)云進(jìn)行精配準(zhǔn)。通過ICP 算法，得到最終精配準(zhǔn)結(jié)果，使手持移動(dòng)激光雷達(dá)點(diǎn)云轉(zhuǎn)換到無人機(jī)激光雷達(dá)點(diǎn)云的同一坐標(biāo)系下。

2.5 融合精度評價(jià)

本研究從2 方面對融合結(jié)果進(jìn)行精度評價(jià)：

1）根據(jù)配準(zhǔn)后2 個(gè)平臺(tái)的樹木位置投影點(diǎn)坐標(biāo)之間的偏移距離評價(jià)。提取手持移動(dòng)激光雷達(dá)點(diǎn)云樹干中心在平面上的投影點(diǎn)坐標(biāo)(XH,YH)及無人機(jī)激光雷達(dá)點(diǎn)云樹冠中心在同一平面上的投影點(diǎn)坐標(biāo)(XU,YU)，計(jì)算無人機(jī)激光雷達(dá)點(diǎn)云樹冠中心點(diǎn)坐標(biāo)相較于手持移動(dòng)激光雷達(dá)點(diǎn)云樹干中心點(diǎn)坐標(biāo)的偏移距離D，計(jì)算公式見式2 以檢驗(yàn)2個(gè)樣地的融合效果。

2）根據(jù)配準(zhǔn)后2 個(gè)平臺(tái)的標(biāo)志物點(diǎn)云數(shù)據(jù)中心點(diǎn)坐標(biāo)之間的均方誤差（mean-square error,MSE）評價(jià)。設(shè)定標(biāo)志物的無人機(jī)激光雷達(dá)數(shù)據(jù)為真值，手持移動(dòng)激光雷達(dá)數(shù)據(jù)為估計(jì)值，由此來計(jì)算均方誤差。均方誤差的計(jì)算公式見式3，其中(xi,yi,zi)為標(biāo)志物的手持移動(dòng)激光雷達(dá)點(diǎn)云中心點(diǎn)坐標(biāo)，(x,y,z)為標(biāo)志物的無人機(jī)激光雷達(dá)點(diǎn)云中心點(diǎn)坐標(biāo)。

3 結(jié)果與分析

3.1 樹木位置提取

經(jīng)過預(yù)處理后，手持移動(dòng)激光雷達(dá)點(diǎn)云與無人機(jī)激光雷達(dá)點(diǎn)云的樹木位置點(diǎn)如圖2所示,與樣地中樹木實(shí)際位置相比較可得分割結(jié)果較為精準(zhǔn)。將圖2a 順時(shí)針旋轉(zhuǎn)約90°后，圖2a 與圖2b 形狀相似，即圖2a 中樹木位置點(diǎn)的分布與圖2b 中樹木位置點(diǎn)的分布大致相同。

圖2 樹木位置點(diǎn)Fig.2 Tree locations

3.2 點(diǎn)云數(shù)據(jù)粗配準(zhǔn)

如圖3所示，使用由2 個(gè)平臺(tái)激光雷達(dá)點(diǎn)云提取的樹木位置點(diǎn)分別構(gòu)造2個(gè)Delaunay三角網(wǎng)。圖3a 三角網(wǎng)中大部分三角形可以在圖3b 三角網(wǎng)中找到對應(yīng)的三角形。Polewski 等通過比較2 種點(diǎn)云數(shù)據(jù)中的單個(gè)位置點(diǎn)與其他位置點(diǎn)間的所有距離來實(shí)現(xiàn)對應(yīng)點(diǎn)的確定[11]，而本研究提出的建立Delaunay 三角網(wǎng)后再比較位置點(diǎn)距離的方法來搜索對應(yīng)點(diǎn)的效率更高。

圖3 樹木位置點(diǎn)構(gòu)造Delaunay 三角網(wǎng)Fig.3 Delaunay triangulation of tree locations

經(jīng)過反復(fù)試驗(yàn)可得，粗配準(zhǔn)中線段長度差閾值為0.8 m 時(shí)，樹木匹配數(shù)量最多且結(jié)果準(zhǔn)確。模擬退火算法中迭代次數(shù)的設(shè)置對于配準(zhǔn)過程的影響較小，如果迭代次數(shù)太小可能會(huì)使得結(jié)果出現(xiàn)偏差，同樣經(jīng)過試驗(yàn)得到設(shè)置迭代次數(shù)為15 時(shí)配準(zhǔn)結(jié)果較準(zhǔn)確且運(yùn)行效率較高。

經(jīng)過粗配準(zhǔn)后，不同平臺(tái)數(shù)據(jù)得到的樹木位置點(diǎn)的對應(yīng)關(guān)系如圖4所示，圖4中同種顏色圈出的點(diǎn)為一組對應(yīng)點(diǎn)，由于樹木傾斜及無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)樹木位置點(diǎn)分割不準(zhǔn)確等原因，部分樹木位置點(diǎn)不能很好地找到其對應(yīng)點(diǎn)，如圖4中用“×”符號(hào)標(biāo)出的樹木位置點(diǎn)。只要求得3組以上對應(yīng)點(diǎn)即可計(jì)算旋轉(zhuǎn)平移矩陣，因此存在未搜索到對應(yīng)點(diǎn)的目標(biāo)點(diǎn)不會(huì)影響粗配準(zhǔn)的進(jìn)行。

圖4 樹木位置點(diǎn)的對應(yīng)關(guān)系Fig.4 The corresponding relationship between tree locations

手持移動(dòng)激光雷達(dá)點(diǎn)云去除高度值后的二維坐標(biāo)利用旋轉(zhuǎn)平移矩陣進(jìn)行二維旋轉(zhuǎn)平移。將旋轉(zhuǎn)平移后得到的粗配準(zhǔn)結(jié)果點(diǎn)云數(shù)據(jù)匹配點(diǎn)云數(shù)據(jù)原有的高度值，從而獲取粗配準(zhǔn)后的三維點(diǎn)。選取樟子松樣地中心位置的任意一棵樹木作為示例，其粗配準(zhǔn)效果如圖5所示，圖中黃色代表手持移動(dòng)激光雷達(dá)點(diǎn)云，藍(lán)色代表無人機(jī)激光雷達(dá)點(diǎn)云。由圖5可得，樹干與樹冠部分均有所偏移。

圖5 粗配準(zhǔn)結(jié)果Fig.5 Rough registration results

3.3 點(diǎn)云數(shù)據(jù)精配準(zhǔn)

針對有偏移、誤差較大的粗配準(zhǔn)結(jié)果，對點(diǎn)云進(jìn)一步進(jìn)行精配準(zhǔn)。選取一定高度區(qū)域的點(diǎn)云數(shù)據(jù)，即手持移動(dòng)激光雷達(dá)點(diǎn)云與無人機(jī)激光雷達(dá)點(diǎn)云均較密集的區(qū)域，運(yùn)用ICP 算法配準(zhǔn)，實(shí)現(xiàn)2 種平臺(tái)點(diǎn)云數(shù)據(jù)的精準(zhǔn)融合。ICP 算法選取的點(diǎn)云區(qū)域與獲取激光雷達(dá)點(diǎn)云數(shù)據(jù)的儀器性能、樹木平均高度以及林區(qū)郁閉度均有很大相關(guān)性。由于手持移動(dòng)激光雷達(dá)點(diǎn)云數(shù)據(jù)與無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)重合度高的區(qū)域?yàn)闃涔趨^(qū)域，因此本研究通過對2 塊樣地中不同樹木點(diǎn)云數(shù)據(jù)高度區(qū)域范圍進(jìn)行選擇，當(dāng)處于平均樹高4～6 m 的范圍時(shí)可以獲取較好的精配準(zhǔn)結(jié)果。

針對樟子松樣地的點(diǎn)云數(shù)據(jù)，選取14～15 m高度區(qū)域的點(diǎn)云數(shù)據(jù)進(jìn)行精配準(zhǔn)。精配準(zhǔn)計(jì)算得到旋轉(zhuǎn)矩陣：平移矩陣：

將粗配準(zhǔn)后的手持移動(dòng)激光雷達(dá)點(diǎn)云利用旋轉(zhuǎn)平移矩陣進(jìn)行坐標(biāo)轉(zhuǎn)換，從而實(shí)現(xiàn)點(diǎn)云數(shù)據(jù)的精配準(zhǔn)。精配準(zhǔn)后，與圖5粗配準(zhǔn)結(jié)果相比，圖6中顯示的同一棵樹的樹冠和樹干部分融合效果得到了改進(jìn)，樹干與樹冠部分基本沒有偏移。

圖6 精配準(zhǔn)結(jié)果Fig.6 Precise registration results

3.4 融合效果與精度

3.4.1 融合效果

經(jīng)過2 次配準(zhǔn)后，基于手持移動(dòng)激光雷達(dá)與無人機(jī)激光雷達(dá)掃描的2 塊樣地的點(diǎn)云數(shù)據(jù)融合效果較好，結(jié)果如圖7所示，圖中藍(lán)色點(diǎn)云為無人機(jī)激光雷達(dá)點(diǎn)云。由圖7可以看出，無人機(jī)激光雷達(dá)點(diǎn)云與手持移動(dòng)激光雷達(dá)點(diǎn)云的樹冠部分基本重合，2 個(gè)平臺(tái)點(diǎn)云相互補(bǔ)充。

圖7顯示了樣地整體的配準(zhǔn)效果，然而為了進(jìn)一步檢驗(yàn)樣地內(nèi)部樹木的融合效果，隨機(jī)選取樣地內(nèi)部區(qū)域的樹木進(jìn)行檢驗(yàn)，選取的區(qū)域見圖8中的方框區(qū)域。

圖7 樣地配準(zhǔn)結(jié)果Fig.7 Registration results of plots

圖8 隨機(jī)區(qū)域選取Fig.8 Selection of random regions

圖9和圖10分別是樟子松樣地與蒙古櫟樣地選取區(qū)域內(nèi)部樹木的激光雷達(dá)點(diǎn)云圖，其中圖9a與圖10a 為無人機(jī)激光雷達(dá)點(diǎn)云圖。由方框區(qū)域可見，無人機(jī)激光雷達(dá)的穿透性使其能夠獲取樣地中樹干的部分點(diǎn)云數(shù)據(jù)。由于樟子松樣地的郁閉度更低且樹木分布更稀疏，與圖10所示蒙古櫟樣地?zé)o人機(jī)激光點(diǎn)云樹干處點(diǎn)云相比，圖9所示樟子松樣地的樹干處點(diǎn)云更加清晰可見。由無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)可以大致確定樹木的生長方向，該方向與圖9b、圖10b 中圖所示手持移動(dòng)激光雷達(dá)點(diǎn)云樹干生長方向基本一致。圖9c、圖10c中黃色點(diǎn)云為手持移動(dòng)激光雷達(dá)點(diǎn)云，藍(lán)色點(diǎn)云為無人機(jī)激光雷達(dá)點(diǎn)云，可以看到隨機(jī)選擇的樣地內(nèi)部區(qū)域中，2 個(gè)平臺(tái)獲取的數(shù)據(jù)在樹干和樹冠區(qū)域均得到有效配準(zhǔn)。

圖9 樟子松樣地Fig.9 Pinus sylvestris plot

圖10 蒙古櫟樣地Fig.10 Quercus mongolica plot

3.4.2 融合精度

1）根據(jù)2 個(gè)平臺(tái)的部分樹木位置坐標(biāo)計(jì)算偏移距離。針對隨機(jī)選擇的區(qū)域，經(jīng)過提取區(qū)域內(nèi)樹木的手持移動(dòng)激光雷達(dá)點(diǎn)云樹干中心及相應(yīng)的無人機(jī)激光雷達(dá)點(diǎn)云樹冠中心，計(jì)算兩個(gè)中心投影點(diǎn)的偏移量可得表1～2，其中XH與YH分別表示手持移動(dòng)激光雷達(dá)點(diǎn)云樹干中心點(diǎn)投影點(diǎn)的平面X 坐標(biāo)與Y 坐標(biāo)，XU與YU分別表示無人機(jī)激光雷達(dá)點(diǎn)云樹冠中心點(diǎn)投影點(diǎn)的平面X 坐標(biāo)與Y坐標(biāo)。通過計(jì)算可得，樟子松的投影點(diǎn)坐標(biāo)偏移距離D 最大和最小值分別為0.30 和0.09 m，而蒙古櫟的投影點(diǎn)坐標(biāo)偏移距離最大和最小值分別為0.59 和0.06 m。

表1 樟子松樣地樹木位置點(diǎn)坐標(biāo)偏移量Table 1 Offset of tree locations in Pinus sylvestris plot

表2 蒙古櫟樣地樹木位置點(diǎn)坐標(biāo)偏移量Table 2 Offset of tree locations in Quercus mongolica plot

2）根據(jù)同一標(biāo)志物在2 個(gè)平臺(tái)中的坐標(biāo)計(jì)算均方誤差。通過計(jì)算樣地上放置的標(biāo)志物在2 種平臺(tái)中的坐標(biāo)之間的偏差來得到點(diǎn)云配準(zhǔn)后的均方誤差，結(jié)果表明，樟子松和蒙古櫟樣地融合后均方誤差分別為0.051 2 和0.080 2，與蒙古櫟樣地相比，樟子松樣地樹木融合精度更高。這可能是因?yàn)檎磷铀蓸拥貫獒樔~林，無人機(jī)激光雷達(dá)點(diǎn)云中樹冠更易被準(zhǔn)確分割，因而由此提取的樹木位置點(diǎn)更加精準(zhǔn)，同時(shí)樟子松樣地的郁閉度較低、樹木分布較離散，所以樟子松樣地的融合精度更高。

4 討 論

與普通測繪相比，林業(yè)測繪具有地形復(fù)雜，樹木樹冠不規(guī)則、非剛體，森林點(diǎn)云特征較普通測繪中建筑物、道路點(diǎn)云特征邊線更難提取等特點(diǎn)，因此特征的選取是森林點(diǎn)云配準(zhǔn)的關(guān)鍵。本研究通過樹木位置構(gòu)成的Delaunay 三角網(wǎng)作為配準(zhǔn)特征，有效解決了森林特征難提取的問題，達(dá)成手持移動(dòng)激光雷達(dá)點(diǎn)云與無人機(jī)激光雷達(dá)點(diǎn)云精確配準(zhǔn)的目標(biāo)。

手持移動(dòng)激光雷達(dá)和無人機(jī)激光雷達(dá)掃描所得數(shù)據(jù)尺度均為樹木實(shí)際尺度，但由于不同平臺(tái)設(shè)備采集數(shù)據(jù)的角度不同，所以不同平臺(tái)點(diǎn)云數(shù)據(jù)的特點(diǎn)也不同，點(diǎn)云的精準(zhǔn)融合需要考慮不同平臺(tái)點(diǎn)云數(shù)據(jù)的特點(diǎn)。手持移動(dòng)激光雷達(dá)自下而上采集數(shù)據(jù)，點(diǎn)云密度自下而上遞減。無人機(jī)激光雷達(dá)自上而下采集數(shù)據(jù)，樹冠處點(diǎn)云較密集，而森林下層數(shù)據(jù)采集數(shù)量較少。2 種設(shè)備的不同特點(diǎn)導(dǎo)致采集到的點(diǎn)云數(shù)據(jù)整體重疊度較低，但是針對樹木中上層點(diǎn)云數(shù)據(jù)，手持移動(dòng)激光雷達(dá)數(shù)據(jù)與無人機(jī)激光雷達(dá)數(shù)據(jù)重疊度較高，較高的重疊度是使用ICP 算法進(jìn)行精配準(zhǔn)的前提。

對于不同的樹種，融合精度有一定差異。針對闊葉林，如本研究中的蒙古櫟，由于郁閉度較高，以及闊葉樹樹冠形狀無規(guī)則的特征，無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)的單木位置提取效果較差，導(dǎo)致后期配準(zhǔn)結(jié)果容易出現(xiàn)誤差。而針對冠幅較小、樹干生長垂直于地面的針葉林，如本研究中的樟子松，無人機(jī)激光雷達(dá)點(diǎn)云數(shù)據(jù)單木分割效果較好，故配準(zhǔn)結(jié)果也更加精準(zhǔn)。

本研究結(jié)果為多源遙感數(shù)據(jù)在林業(yè)中的廣泛應(yīng)用提供了可能，有利于克服單一平臺(tái)無法獲取樹木完整點(diǎn)云的局限性，對于森林結(jié)構(gòu)參數(shù)的精確提取和樹木三維模型的精準(zhǔn)構(gòu)建有重要意義，但仍存在以下不足：1）本研究區(qū)域?yàn)槿斯ち郑匦纹骄?，樹木種類較單一，針對樹種及地形更加復(fù)雜的天然林，本算法的適用性仍有待驗(yàn)證。2）雖然模擬退火算法使得融合方法運(yùn)行效率有所提升，但當(dāng)本方法應(yīng)用于大區(qū)域尺度時(shí)，仍存在花費(fèi)時(shí)間較長的問題。同時(shí)隨著區(qū)域的增大，樹木位置對應(yīng)點(diǎn)的搜索結(jié)果也更易出現(xiàn)誤差。3）本研究雖然將手持移動(dòng)激光雷達(dá)數(shù)據(jù)與無人機(jī)激光雷達(dá)數(shù)據(jù)精準(zhǔn)融合，但是并未對融合后的數(shù)據(jù)進(jìn)一步精簡處理，數(shù)據(jù)仍存在冗余。

在未來研究中將嘗試選取地形起伏較大、林分結(jié)構(gòu)更復(fù)雜的天然林作為研究對象，進(jìn)一步驗(yàn)證算法的適用性。另外，對算法進(jìn)行優(yōu)化以提高算法的自動(dòng)化程度及運(yùn)行效率，減少時(shí)間成本。針對融合后數(shù)據(jù)冗余問題，基于森林結(jié)構(gòu)特點(diǎn)提出點(diǎn)云壓縮算法，為后續(xù)提取高精度森林結(jié)構(gòu)參數(shù)、構(gòu)建三維可視化模型奠定基礎(chǔ)。

5 結(jié) 論

本研究利用樹木位置點(diǎn)，提出了一種基于Delaunay 三角網(wǎng)和ICP 算法的不同平臺(tái)點(diǎn)云數(shù)據(jù)的融合方法，該方法屬于無標(biāo)識(shí)配準(zhǔn)方法，并通過2 塊復(fù)雜程度不同的樣地獲取得到的點(diǎn)云數(shù)據(jù)進(jìn)行了驗(yàn)證。結(jié)果表明，在不同的參數(shù)設(shè)置下，樹木匹配的精度可能不同。粗配準(zhǔn)中設(shè)置的線段差閾值和精配準(zhǔn)中設(shè)置的高度區(qū)域范圍這2 個(gè)參數(shù)直接影響配準(zhǔn)結(jié)果的精準(zhǔn)度。而粗配準(zhǔn)中設(shè)置的迭代次數(shù)的變化對精準(zhǔn)度結(jié)果的影響較小，但是迭代次數(shù)的合理選擇可以提高方法運(yùn)行效率。

本研究方法可以較為準(zhǔn)確地將手持移動(dòng)激光雷達(dá)點(diǎn)云和無人機(jī)激光雷達(dá)點(diǎn)云進(jìn)行融合，融合效果較好。融合后2 塊樣地隨機(jī)選取區(qū)域中的2種平臺(tái)數(shù)據(jù)樹木位置投影點(diǎn)的偏移距離均值分別為0.19 和0.25 m。根據(jù)配準(zhǔn)后標(biāo)志物在2 種平臺(tái)中的位置坐標(biāo)，計(jì)算可得融合后2 塊樣地的均方誤差均小于0.1，融合精度較高。本研究結(jié)果可為單木結(jié)構(gòu)參數(shù)精確提取、樹木地上生物量估算以及三維模型精準(zhǔn)構(gòu)建提供基礎(chǔ)，為多源數(shù)據(jù)在林業(yè)中的應(yīng)用提供參考。