應用空載全波形光達資料於波形分析與地物分類

國 立 交 通 大 學

土 木 工 程 學 系

碩 士 論 文

應用空載全波形光達資料於波形分析

與地物分類

Waveform Analysis and Landcover Classification

Using Airborne Full-Waveform Lidar Data

研 究 生：林郁珊

指導教授：張智安

應用空載全波形光達資料於波形分析與地物分類

Waveform Analysis and Landcover Classification Using

Airborne Full-Waveform Lidar Data

研 究 生：林郁珊

Student：Yu-Shan Lin

指導教授：張智安

Advisor：Tee-Ann Teo

國 立 交 通 大 學

土 木 工 程 學 系

碩 士 論 文

A Thesis

Submitted to Department of Civil Engineering College of Engineering

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master in

Civil Engineering July 2012

Hsinchu, Taiwan, Republic of China

i

應用空載全波形光達資料於波形分析與地物分類

學生：林郁珊 指導教授：張智安

國立交通大學土木工程學系

中文摘要

空載光達是一種主動式對地進行測量的技術，可快速取得地表物之三維坐標，並提 供高密度及高精度之三維坐標。而全波形光達技術中重要的發展，其可以記錄回波的連 續波形，可藉由分析波形得到更多的地表反射物理特性與地表細節與變化，讓使用者擁 有較豐富及完整的資訊，有助於地形重建及地物判識。 在處理全波形光達資料的過程中，波形萃取與分析是首先且必要的一環。本研究的 第一階段，先進行資料整理與檔案格式讀取，接著利用對稱與不對稱兩種函數(高斯函 數及韋伯函數)對波形進行擬合(fitting)，並進行原始資料與擬合成果兩者間的精度評估， 討論與驗證不同擬合函數對於全波形光達訊號處理的適切性。經過精度評估後，萃取光 達 波 形 參 數 包 括 波 寬 (echo width) 、 振 幅 (amplitude) 、 背 向 散 射 參 數 (backscatter cross-section coefficient)，配合光達幾何參數包括高程(Z)、高程差(dZ)、回波數(echo number)、多重回波百分比(echo ratio)，提供為接下來分析應用的特徵。 取得波形的參數後，接著將這些萃取得的波形參數作為分類之用。所以本研究的第 二階段，就萃取出的波形參數配合人工判識資料，使用支持式向量機(SVM)與隨機森林 (Random Forest)兩種方法做地表物的特徵分類，並就分類成果進行精度評估，藉此探討 使用全波形光達資料與一般空載光達進行分類與兩種分類法的優缺點。 本研究結果顯示，波形擬合分析的部分，使用不對稱函數的擬合成果較好，擬合誤 差較小；但是在於波形峰值的位置與對稱函數的成果差異不大，因此高斯函數為一個簡

ii

單省時的擬合方式。而在地物分類的部分，全波形光達所提供的參數對於樹木的分析的 結果顯著；相較於支持式向量機，隨機森林分類法的成果較為佳。

iii

Waveform Analysis and Landcover Classification Using Airborne

Full-Waveform Lidar Data

Student：Yu-Shan Lin Advisor：Tee-Ann Teo

Department of Civil Engineering

National Chiao Tung University

Abstract

Airborne Lidar is an active remote sensing system. It can obtain the three dimensional coordinates effectively, and provide high density and high precision 3-D point cloud. Full-waveform (FWF) lidar is a new generation of airborne laser scanner which receives one dimensional continuous signal. It offers useful information about the structure of the target. Therefore, the analysis of received signal of FWF lidar and obtaining the implicit information is helpful for landcover classification.

In the processing of full waveform Lidar data, the waveform parameter extraction and analysis are the important steps. The major objective of this study is to analyze the received waveform and extract its parameters. We select Gaussian distribution as a symmetric function and Weibull distribution as an asymmetric function in waveform decomposition. Then, we calculate several accuracy assessment indicators between raw waveform data and fitting function for quality assessment. We use echo width, amplitude, backscatter cross-section coefficient, elevation, elevation difference, echo number, and echo ratio as waveform parameter of classification.

After waveform parameter extraction, we select Support Vector Machine (SVM) and Random Forests (RF) as classifier for landcover classification. This study employs echo width,

iv

amplitude, backscatter cross-section coefficients and other features for classification. Error matrix is used to compare the performance of the classifiers.

The experimental results indicate that the accuracy of asymmetric function is slightly better than symmetric function. However, the extracted peak positions from the Gaussian and Weibull are very close. Moreover, Gaussian distribution is relatively simple and easy to implement in the waveform analysis. The result of landcover classification shows that waveform parameters are helpful for classification and Random Forests classifier is better than SVM in our study cases.

Keywords : Airborne Lidar, Full-Waveform, waveform fitting, landcover classification, Support Vector Machine (SVM) , Random Forests (RF)

v

致謝

短暫的研究所生涯匆匆的過去了，在這交大最後的兩年日子裡，每天都過得相當的 充實與開心。很幸運可以在大四的時候認識張智安老師，成為老師在交大第一屆的專題 生，接著也成為老師第二屆的研究生。很感謝與老師相處的這些日子，老師總是不厭其 煩的教授我們做研究的態度與方法，不論再忙也都會關心與照顧我們，沒有老師的指導 與幫助絕對無法完成此碩士論文。另外也感謝口試委員林玉菁老師、林士淵老師，為我 的論文提供了相當多寶貴的建議，使學生的論文能夠更加的完善。此外還有感謝平時很 關心我的史天元老師，與史老師的團隊一起meeting總是能激發思考力。 感謝悶悶學長，總是照顧我幫我解決很多的問題，祝你事事順心。感謝信瑜學長兩 年的照顧，幫我度過大四與碩一的程式海。感謝暐尊學長對我研究主題上的幫忙。感謝 三大正妹蒨蒨肥肥昭昭為我的研究生涯帶來了許多歡笑與淚水，能夠在這段日子認識妳 們有妳們相伴，讓我並不孤單!!感謝羊咩大師濃濃翟翟同儕們各方面的幫忙。學弟以諾 boy認真向上的研究精神激發了我的動力、宛宜學妹天天都有的笑料感謝妳幫助我校稿、 麵包碩碩大任讓 413 門庭若市。另外還有很多很多的好朋友、大學同學、高中同學、學 長姊學弟妹，你們的相伴支持著我，讓我有源源不絕的動力完成研究與論文。 感謝爸爸媽媽支持我唸研究所，總是聽久久才回家一次的女兒抱怨東抱怨西，買好 吃的東西給女兒吃，讓女兒帶滿滿的食物與愛心離家繼續做研究。感謝阿姊的照顧，每 次放假最期待的就是去高雄找妳住，謝謝妳總是能夠包容我的壞習慣與壞口氣，能跟妳 同時拿到碩士學位真的很開心。最後感謝死神=吸血鬼=阿吸總是一直陪在我身邊，不論 我生氣開心或是生病，總是照顧著我，讓我們一起順利的完成研究所的最後階段吧!! 謝謝每個曾經讓我感受到溫暖與幫助的人 郁珊 2012 年夏天

vi

目錄

中文摘要 ... i Abstract ... iii 致謝 ... v 目錄 ... vi 圖目錄 ... ix 表目錄 ... xii 第 1 章 前言 ... 1 1.1 研究背景動機 ... 1 1.2 研究目的 ... 2 1.3 研究方法與流程 ... 3 1.4 論文架構 ... 3 第 2 章 全波形光達相關研究 ... 4 2.1 文獻回顧 ... 4 2.1.1 全波形光達原理與優點 ... 4 2.1.2 擬合函數模型選擇與求解 ... 7 2.1.3 全波形光達應用於森林與都市地物分類 ... 11 2.1.4 支持向量機(SVM)、隨機森林(Random Forest)分類器應用 ... 15 2.2 全波形光達儀器介紹 ... 16 2.2.1 Leica ALS系列 ... 17 2.2.2 Riegl LMS系列 ... 17 2.2.3 Optech ALTM3100 ... 18 2.2.4 Trimble Harrier56 ... 19

vii 2.3 光達資料格式介紹 ... 20 2.3.1 LAS1.3 格式 ... 20 2.3.2 SDW格式 ... 22 第 3 章 研究方法 ... 23 3.1 波形分析 ... 23 3.1.1 資料前處理 ... 24 3.1.2 計算初始值 ... 25 3.1.3 擬合函數計算 ... 27 3.1.3.1 波形分解、光達方程式與背向散射參數 ... 27 3.1.3.2 對稱與非對稱函數 ... 31

3.1.3.3 Trust Region algorithm ... 33

3.1.4 精度分析 ... 34

3.2 地物分類 ... 35

3.2.1 分類特徵參數 ... 35

3.2.2 分類方法 ... 37

3.2.2.1 隨機森林分類(Random Forests) ... 37

3.2.2.1 Support Vector Machine (SVM)分類 ... 39

3.2.3 精度評估 ... 42 第 4 章 研究資料 ... 43 4.1 測區一 ... 43 4.2 測區二 ... 43 4.3 測區三 ... 46 第 5 章 成果分析 ... 49 5.1 波形模擬 ... 49

viii 5.1.1 模擬變數設定 ... 49 5.1.2 模擬成果分析與討論 ... 51 5.2 波形分析 ... 55 5.2.1 波形與地物 ... 55 5.2.2 擬合成果 ... 57 5.2.3 波形參數圖 ... 60 5.2.4 背向散射參數改正 ... 64 5.3 分類成果 ... 66 第 6 章 結論與建議 ... 77 6.1 結論 ... 77 6.2 建議 ... 78 參考文獻 ... 79

ix

圖目錄

圖 2-1 脈衝式多重回波光達原理(Mallet and Bretar, 2009) ... 4

圖 2-2 全波形光達回波與取樣 (Doneus et al., 2008) ... 5

圖 2-3 多重回波光達與全波形光達點雲剖面(Chauve et al., 2007) ... 5

圖 2-4 地物辨識與足跡波長的關係(Jutzi and Stilla, 2005) ... 6

圖 2-5 地物垂直分布分辨與足跡大小(Mallet and Bretar, 2009) ... 6

圖 2-6 不同地形地物的回波形狀(Jutzi and Stilla, 2006) ... 7

圖 2-7 左：高斯分佈與對數常態分佈、右：廣義高斯分佈形狀參數α的影響 .. 8

圖 2-9 使用不同Level的小波擬合全波形光達--藍：波形，紅：擬合成果 ... 9

圖 2-11 全波形與多重回波光達CHM的差異(Chauve et al., 2007a) ... 11

圖 2-12 地表粗糙度暈渲圖與實際地物(Hollaus and Höfle, 2010) ... 12

圖 2-13 單顆樹樹幹偵測與重建成果(Reitberger et al., 2007)... 12 圖 2-14 分類參數與分類成果的關係(Mallet et al., 2008) ... 13 圖 2-15 波形參數對於不同地物類別的重要性(Guo et al., 2011)... 14 圖 2-16 全波形與脈衝式光達對於地物分類之精度評估(Mallet et al., 2011) ... 15 圖 2-17 Leica ALS60 ... 17 圖 2-18 旋轉多邊形掃描原理 ... 18 圖 2-19 LMS-Q560(左)、Q680i(右) ... 18 圖 2-21 Trimble的Harrier56 ... 19 圖 2-22 LAS 1.0 與LAS 1.3 ... 21 圖 2-23 LAS公開檔頭資訊範例 ... 21 圖 2-24 波形數據紀錄waveform package範例... 21 圖 2-25 SDW檔案內容範例 ... 22

x 圖 3-1 波形分析研究流程 ... 23 圖 3-2 高斯平滑前後成果 ... 24 圖 3-3 濾除背景雜訊前後成果 ... 25 圖 3-4 有錯誤的初始峰值 ... 26 圖 3-5 使用門檻濾除錯誤的初始峰值與擬合成果 ... 26 圖 3-6 波形分解圖例 ... 27 圖 3-8 背向散射示意圖(Wagner et al., 2006) ... 28 圖 3-9 高斯函數 (上) 與韋伯函數 (下) ... 32 圖 3-10 對稱與不對稱函數之差異 ... 33 圖 3-11 地物分類流程圖 ... 35 圖 3-12 地形粗糙度的參數dZ ... 36 圖 3-13 考慮地形影響的dZ ... 36 圖 3-14 決策樹分類(Tan et al., 2006) ... 37 圖 3-15 隨機森林分類步驟示意圖(Guo et al., 2011)... 38 圖 3-16 SVM原理示意圖(線性分割) ... 40 圖 3-17 非線性屬性轉換示意圖 ... 40 圖 4-1 測區一點雲影像及四角坐標 ... 44 圖 4-2 測區二點雲影像及四角坐標 ... 45 圖 4-3 測區二之參考地物類別資料 ... 46 圖 4-4 測區三點雲影像及四角坐標 ... 47 圖 4-5 測區三之參考地物類別資料 ... 47 圖 5-1 模擬波形示意圖 ... 50 圖 5-2 模擬波形分辨力示意圖(左：可以分辨、右：無法分辨) ... 50 圖 5-3 noise deviation=1 ... 53

xi 圖 5-4 noise deviation=2 ... 53 圖 5-5 noise deviation=3 ... 53 圖 5-6 不同波寬與noise deviation相互影響關係圖... 54 圖 5-7 裸露地與其回波波形 ... 55 圖 5-8 波形與地物的關係(測區一) ... 56 圖 5-9 測區一RMSE分布直方圖 ... 58 圖 5-10 萃取更多回波的成果(測區三) ... 60 圖 5-11 回波太接近的成果(測區三) ... 60 圖 5-12 波形參數圖(測區一) ... 61 圖 5-14 波形參數圖(測區三) ... 63 圖 5-16 測區三全區背向散射參數改正成果(左：改正前，右：改正後) ... 65 圖 5-17 改正參數Ccal值分布(左：測區二，右：測區三) ... 65 圖 5-18 測區背向散射參數圖(左：測區二，右：測區三) ... 65 圖 5-19 測區二分類成果 ... 73 圖 5-20 測區三分類成果 ... 75 圖 5-21 測區二隨機森林分類特徵重要性 ... 76 圖 5-22 測區三隨機森林分類特徵重要性 ... 76

xii

表目錄

表 2-1Generalize Gaussian、Weibull、Nakagami、Burr之差異 (Mallet et al.,

2009) ... 9 表 2-2 全波形光達儀器參數 ... 16 表 4-1 研究資料與實驗內容 ... 48 表 5-1 模擬成果之分辨力大小(單位：與波寬之比值) ... 52 表 5-2 模擬成果之分辨力大小(單位：距離(m)) ... 52 表 5-3 測區一擬合成果的XYZ與原始資料的差異 (單位:公尺) ... 57 表 5-4 測區一精度評估成果 ... 58 表 5-5 測區一，三種不同地物的精度評估 ... 59 表 5-6 測區三高斯函數擬合成果 ... 59 表 5-7 測區二精度成果(使用振幅未使用背向散射參數) ... 67 表 5-8 測區二精度成果(使用背向散射參數未使用振幅) ... 68 表 5-9 測區三精度成果(使用振幅未使用背向散射參數) ... 69 表 5-10 測區三精度成果(使用背向散射參數未使用振幅) ... 70 表 5-11 測區三分類成果誤差矩陣 ... 71

1

第1章 前言

1.1 研究背景動機

空載光達是一種主動式的對地量測系統，一般脈衝式空載光達系統的原理是以飛機 做為載具，搭載全球衛星定位系統(Global Positioning System, GPS)與慣性量測裝置 (Inertial measurement unit, IMU)，獲取定位資訊與姿態參數。利用脈衝雷射進行掃描， 記錄反射訊號的回訊脈衝及時間再轉換為距離，進行直接地理定位。其高程精度可達 0.15 公尺，可以將大量的點位坐標(點雲)的地理幾何資訊與反射強度資料，運用在很多 不同的領域，如建立數值高程模型(Digital Elevation Model, DEM、Digital Terrain Model, DTM)、三維空間模型、電力線偵測等等應用。 隨著科技技術的進步與硬體儲存設備上面的演進，空載陸域全波形光達從 2004 年 開始商品化及普及化，目前多數空載光達系統均具有獲取全波形資料的能力。全波形光 達可以接收與紀錄連續的波形，使用者可以自行對波形進行後處理，獲得所需要的資訊。 全波形光達與傳統脈衝式多重回波光達系統不同地方在於，多重回波光達將接收到的訊 號以門檻值過濾，輸出成為多個回波，只記錄該點的強度值與三維離散坐標，其優點是 提升存取的效能並節省儲存空間；而全波形光達藉由回波訊號的密集的取樣，取樣間距 可達 1ns記錄 1 筆回波訊號，以描述接收到連續的回波訊號，能夠提供使用者比傳統多 重回波光達更多的地表資訊與更高密度的點雲，有助於地形重建及地物判識。 全波形光達依其使用目的主要分為以下三種類別：(1)實驗光達(Experimental Lidar Systems) 、 (2) 測 深 光 達 (Bathymetric Lidar Systems) 、 (3) 商 用 光 達 (Commercial Lidar Systems)。實驗光達設計目的為科學任務，主要載具是衛星，雷射足跡(Footprint)直徑約 為在 10 至 100 公尺之間，例如NASA的SLICER (Scanning Lidar Imager of Canopies by Echo Recovery)、SLA (Shuttle Laser Altimeter)、LVIS (Laser Vegetation Imaging Sensor)、 MOLA (Mars Orbiter Laser Altimeter)與ICESat (Ice, Cloud,and land Elevation Satellite)任

2

務，對全球進行地形、植被等環境資料收集。測深光達主要目的是進行水深的量測，利 用紅外光(1064nm)及綠光(532nm)分別量測水表及水底，例如Fugro LADS、Optech SHOALS、AHAB HawkEye等，主要對淺海水質清澈地區量測水深及海底地形。商用光 達大多是進行陸域掃描，利用掃描頻率高與掃描足跡小的特點，提供豐富且細緻的地形 資訊，例如Riegl Q680i、Trimble (Toposys)、Optech ALTM Gemini、Leica ALS系列等。

本研究的研究動機是藉由波形訊號提升光達資料處理品質。在波形數據中，不同的 地物所反射的回波形狀不同，若是可以利用「分析波形」或是「參數萃取」等方法來獲 得全波形光達數據中所隱含的資訊，相信能獲得更詳細清楚的地形地貌，這也是目前全 波形光達相關研究中十分重要的議題。所以本研究希望可以藉由數學函式擬合全波形光 達數據，找出其波形的代表參數；同時比較不同的函數對於波形擬合的成果是否有顯著 影響。並且最後在地物分類的過程中加入波形的參數，藉以提升地物分類之精度。

1.2 研究目的

本研究之目的包含波形擬合及地物分類，波形擬合中比較對稱與不對稱函數的擬合 成果，進行分析；接著整合這些萃取出來的波形參數與光達幾何資訊，進行不同地物的 分類。 波形擬合的目的是為了找出適合擬合的函數，以萃取幾何及波形特徵，本研究將使 用不同的函數進行回波的分解，並且比較方法的優缺點。採用的函數包含：(1)對稱高 斯函數、(2)非對稱Weibull函數。 地物的分類研究重點在於分析全波形特徵對分類的幫助，比較使用全波形資料加入 分類與未使用全波形資料分類的精度評估成果，利用不同的分類法分析全波形資料的必 要性與其對於不同地表物分類的幫助。

3

1.3 研究方法與流程

本研究的方法主要分為兩部分：第一部分為全波形光達的波形分析，包含資料前處 理、計算初始值、函數擬合、參數萃取與精度分析工作；第二部分為應用全波形光達於 地物分類，包含建立分類圖層、隨機森林(Random Forest)分類、支持式向量機(Support Vector Machine, SVM)分類與精度評估工作。 全波形光達波形分析使用三個不同測區、不同地表物、兩種儀器所獲得的資料。從 資料前處理的部分著手進行，依據不同的使用資料給定不同的高斯平滑參數和濾除背景 雜訊。接著對平滑處理過的波形微分進行區域最大值的偵測，其目的是為了當作接下來 擬合解算的初始值。在擬合的函數方面，本研究使用了高斯函數與韋伯函數，主要的因 素是為了探討對稱與不對稱函數對於波形擬合成果的影響。 地物分類的部分，分類的特徵包含全波形光達之三維點雲坐標、高程差、波寬、振 幅、回波數目等等分類特徵。將這些特徵整合不同圖層的網格式資料，同一網格位置包 含不同圖層的特徵。本研究選用的分類器為隨機森林與SVM兩種，最後進行分類的精度 評估與波形特徵之重要性探討。

1.4 論文架構

第一章為本研究之背景、動機與目的，研究方法與流程之簡述。 第二章為相關研究的部分：(1)文獻回顧包含全波形光達原理、擬合函數模型的選擇、全 波形光達的應用和隨機森林與SVM分類器之應用；(2)全波形光達儀器則是介紹目前商 用陸域全波形光達；(3)全波形光達數據資料的LAS標準格式介紹。 第三章為研究之方法，說明波形分析與地物分類方法。 第四章為研究資料的部分，使用三個測區分別代表不同的地表覆蓋物。 第五章為成果分析，分波形模擬參數設定、函數擬合波形分析與地物分類三大項目。 第六章為結論，分別就不同的地表覆蓋物所做的分析討論做總結。

4

第2章 全波形光達相關研究

2.1 文獻回顧

本節是研究相關的文獻回顧，包含全波形光達原理與優點、擬合函數模型選擇、全 波形光達在森林與都市區域分類的應用和隨機森林、SVM分類器的優點與應用。

2.1.1 全波形光達原理與優點

傳統的脈衝式多重回波光達其原理是將接收到的訊號以門檻值過濾，門檻值以上的 才會被記錄下來，輸出成為多個離散的回波，如圖 2-1。多重回波光達數據只記錄點雲 強度值與三維坐標，且其記錄的回波數量有限制，因此使用者沒有辦法取得最原始的完 整資訊。 然而，在硬體儲存設備的提升下，光達一次的飛行計畫中可以儲存更多的資料，得 以讓波形資料被記錄下來。因此紀錄回波波形的方式是經由密集的取樣，紀錄連續的 DN(digital number)值，這樣就可以得到回波波形如圖 2-2，通常取樣間隔為 1~3ns。而 不同的地物所反射的回波形狀也不同，使用者可以自行對波形進行後處理，獲得所需要 的幾何與物理資訊。

5

圖 2-2 全波形光達回波與取樣 (Doneus et al., 2008)

由於全波形光達系統儲存的資料是回波波形資料，與多重回波光達相較之下，其回 波資料較為完整，因此可將波形的峰值位置萃取出來，就能夠得到更為密集的點雲，也 可以利用於地面點的分類。研究成果顯示，波形資料可以偵測出比傳統多重回波光達多 30%到 130%的點，當植物密度愈高的地區其增加的點數愈多(Chauve et al., 2007a) 而這 些多偵測出的點通常位於茂密的樹冠與地面點，見圖 2-3。因此全波形資料具有優勢的 點雲密度，對於地面點的偵測是有很大的幫助。

圖 2-3 多重回波光達與全波形光達點雲剖面(Chauve et al., 2007)

全波形光達的足跡(footprint)與光達的光束發散角和飛行高度有關。不同大小的足跡 會影響地表物的解析度。Hollaus and Höfle(2010)說明地形粗糙度與足跡大小和波長有關， 分為三個等級：(1)Macro structure、(2)Meso structure、(3)Micro structure，分別可以判釋

6

出不同的屋頂類型、地面高程不同的小結構物與表面材質粗糙度，見圖 2-4。

圖 2-4 地物辨識與足跡波長的關係(Jutzi and Stilla, 2005)

足跡同時也會影響地物垂直分布的分辨力，以小足跡(Small footprint, 0.2-3m)光達系 統為例，可能會錯過樹冠層與樹梢，造成光達紀錄的高度與實際高度差異，也會較難分 辨地面點(Dubayah and Drake, 2000)，因此小足跡的系統必須提升掃描頻率，以提升密度。 大足跡(Large-footprint,10-70 m) 光達系統可以增加雷射掃描到樹頂跟地面的機率，獲得 較完整的地貌資訊(Mallet and Bretar, 2009)，見圖 2-5。不過相對來說，大足跡的光達在 光跡內涵蓋的地表資訊較為複雜，比較難以精確的分析與辨識地表覆蓋物，且其波形複 雜不易分析，所以目前商用的陸域光達多屬於小足跡的系統。

7

2.1.2 擬合函數模型選擇與求解

全波形光達研究中，波形分析是一個重要的議題，波形分析之目的為了解與獲得全 波形光達數據中所隱含的資訊。因為接收到的訊號是不同地物回訊的疊置，因此可使用 數學函數來擬合及分解波形，數學函數中的參數可視為回波波形的代表因子，分析這些 不同的因子所代表的不同幾何意義，研究其對應到的地物類別與地形表現是一個重要的 議題 圖 2-6 為當光達掃描到不同地物之回波形狀分布情形示意圖，當光達掃描到不同地 物、傾斜、高度差異或粗糙的表面都會影響回波的形狀。而c、d圖則是代表高程差異愈 小，回波就會愈靠近，不易分辨(Jutzi and Stilla, 2006)。

圖 2-6 不同地形地物的回波形狀(Jutzi and Stilla, 2006)

由於大部分的光達發射的脈衝訊號接近於高斯分布，因此假設接收之回波波形也近 似高斯分佈，其函數主要包含振幅(amplitude)與波寬(width)兩個參數，可以用來量化某 一地物波形。藉由使用多個高斯函數疊加擬合波形訊號，可以快速地將資訊用兩個參數 來表示，所以高斯分解法廣泛應用於全波形之波形分析中(Chauve et al., 2007b; Hofton et al., 2000; Jutzi and Stilla, 2006; Wagner et al., 2006)。

另外回波波形由於時間延遲的現象，會有些微不對稱的情形發生，所以波形分解亦 可使用不對稱的函式進行擬合，並探討使用不對稱的函式之成果。Chauve et al.(2007)比 較三種對稱與不對稱函數模型：(1)高斯分佈(Gaussian)、(2)對數常態分佈(Log normal)、 (3)廣義高斯分佈(Generalize Gaussian)，如圖 2-7，三種函數分別代表的是對稱、不對稱

8 與分佈形狀。對數常態分佈與廣義高斯分佈可提升波形擬合成果，但是某些比較不對稱 的回波波形其實是因為兩個波距離太接近疊加而成，所以不對稱的對數常態分佈會使得 萃取出來的點減少。廣義高斯分佈則是對稱的函數，該函數提供形狀參數α代表波形， 在圖中可以看到α=√ 就是高斯分布(紅線)，α愈小形狀愈尖(藍線)，α愈大形狀愈圓(綠線)。 此特性對於分類來說非常的有用，不同的形狀參數代表不同的地物特性，見圖 2-9，從 統計結果可發現建物屋頂(黑色)、柏油路(深灰)與茂密植物(灰色)的分布十分的不同。 圖 2-7 左：高斯分佈與對數常態分佈、右：廣義高斯分佈形狀參數α的影響 (Chauve et al., 2007b) 圖 2-8 α值統計灰度直方圖：建物屋頂(黑色)，柏油路(深灰)，茂密植物(灰色) (Chauve et al., 2007b) Mallet et al. (2009)則是比較廣義高斯分佈(Generalize Gaussian)與其他三種不對稱的 函數模型Weibull、Nakagami、Burr共四種函數的成果，其參數與函式如下表 2-1。此研 究分別使用了模擬波形、測深光達資料與大足跡的衛載實驗光達(SLICER、LVIS)資料、 小足跡的商用陸域光達資料進行實驗。實驗結果顯示，不同的光達資料所適合的擬合函

9

數皆不相同，因此必須視本身資料的特性來選擇適合的函數。

表 2-1Generalize Gaussian、Weibull、Nakagami、Burr之差異 (Mallet et al., 2009)

另外擬合波形亦可使用小波(wavelet)的方式(Laky et al., 2010)，從圖 2-9 可以看到小 波的Level愈高，擬合成果(紅線)愈好，相對的受到雜訊的干擾也會較嚴重。B-spline也 被提出可以有效的對波形進行擬合(Roncat et al., 2011)，其成果與高斯函數相較，可以找 到更多疊合的回波(圖 2-10)，但可能發生過度分割的問題，容易出現不合理的回波參數， 跟小波一樣都要濾除不合理的擬合成果。 圖 2-9 使用不同Level的小波擬合全波形光達--藍：波形，紅：擬合成果 (Laky et al., 2010)

10

圖 2-10 左：B-spline成果；右：高斯成果(Roncat et al., 2011)

雖然擬合回波的函式有很多，也各有優缺點，但過去文獻中的文章中認為高斯分布 是一個普遍且簡單的對稱函式(Wagner et al., 2006)，也可求得較佳的成果(Jutzi and Stilla, 2005)，因此多數波形分析法以高斯模型為基礎。在不對稱函式函數中，Weibull是一個 普遍且相對簡易並具有意義的不對稱函式(Coops et al., 2007)。

由於擬合函數多是非線性函數，需要進行迭代求解。非線性函數的迭代求解方法包 含 最 大 相 似 的 Expectation Maximation 演 算 法 (Dempster et al., 1977) 、 非 線 性 的 Levenberg-Marquardt (Levenberg, 1944; Marquardt, 1963)、Trust-region演算法(Lin et al., 2008) 或Reversible Jump Markov Chain Monte Carlo演算法 (Green, 1995)等，進行擬合函 數參數求解。在非線性函數求解中，初始值的給定是一項重要的因素。

11

2.1.3 全波形光達應用於森林與都市地物分類

許多研究指出全波形光達可提高地物分類的準確度，獲得更細微的地表物參數，最 常被使用在森林區與都市區植被偵測。除此之外，也有研究對於全波形光達資料於地物 分類的幫助進行量化性的評估。

(一)森林區

在森林區域全波形光達所提供的好處是能夠得到較為精準的樹冠層模型(Canopy Height Model, CHM)，圖 2-11 為全波形光達成果與多重回波光達成果的差異，上面為多 重回波光達，下面為全波形光達，左邊是立體側視圖，右邊是俯視圖，可以看到樹木高 度差異愈大愈複雜的區域，其樹冠高度會差愈多(Chauve et al., 2007a)。

圖 2-11 全波形與多重回波光達CHM的差異(Chauve et al., 2007a)

全波形光達對森林地區的樹木參數萃取也有很大的幫助，萃取的參數有波寬、振幅 與背向散射截面參數，這些參數可以用來分出不同樹種，例如：Larch、Oak、Beech (Höfle et al., 2008)，對於針葉林、茂密樹冠與落葉林也有顯著差別(Reitberger et al., 2008)，針 葉林與闊葉林的分類總體精度更可以高達 91%(Heinzel and Koch, 2011)。

12

表示地表粗糙度的參數通常為點雲高程的差異值(dZ)，但研究顯示全波形光達的波 寬參數對粗糙度分析很有幫助，雜訊較dZ少很多。從圖 2-12 之地表粗糙度暈渲圖可以 看到細緻的資訊，例如小徑、傾倒的樹幹、灌木叢等等(Hollaus and Höfle, 2010)。

圖 2-12 地表粗糙度暈渲圖與實際地物(Hollaus and Höfle, 2010)

另外還能對樹木區塊與單顆樹樹幹偵測(Reitberger et al., 2007)，這項研究提出了一 種新的單一樹種檢測，結合表面重建和樹幹檢測。在複雜森林類型的結果表明，可以在 上層和中間層改善單顆樹樹幹偵測，圖 2-13。

13

(二)都市區

都市區的植物亦可使用波形參數萃取樹木參數(Höfle and Hollaus, 2010)，樹木參數 最高點(樹高)、樹冠直徑、樹幹。主要研究工作分別為植物區域偵測、獨立樹木分割及 參數萃取。研究成果顯示樹木與建物的分類精度很高，最大的錯誤率是 14%，且多位於 建物的邊緣。 波形參數對分類成果有顯著提升(Mallet et al., 2008)，圖 2-14 為分類成果的總體精 度，後面三個參數（振幅A、波寬σ及背向反射α）是從全波形光達資料所萃取，增加波 形參數之後總體精度有上升的趨勢。 圖 2-14 分類參數與分類成果的關係(Mallet et al., 2008) Guo et al. (2011)藉由使用可以指出分類特徵參數的分類器(隨機森林)，對都市區進 行分類，波形參數對於不同地物類別的重要性是不同的，見圖 2-15(A, w, σ, α皆為波形 參數)。對於都市區域來說，波形參數的影響量是顯著的，僅次於高程差(Δz)；但對人 造地面來說，波形參數重要性較低。因此，圖上所表達的意義是全波形光達參數對不同 種類的地物其影響大小也不同，有些種類表現突出但是有些則是不顯著。

14

圖 2-15 波形參數對於不同地物類別的重要性(Guo et al., 2011)

Höfle et al (2012)運用全波形光達的參數與不同分類方式(類神經網路artificial neural network、決策樹decision tree)進行都市區域植被偵測，其分別統計比較與探討多重回波、 網格化分割、輻射校正與分類法對於都市區植被偵測結果的影響。其研究成果顯示，全 波形光達資料增加描述物體特徵的條件，可以有效幫助地物分類，但是仍然需要配合物 件導向的網格資料(object-based raster)與三維點雲分析(3D point cloud analysis)才能使其 發揮其最大的優勢(Höfle et al., 2012)。 Mallet et al. (2011)量化評估全波形光達資料在都市區分類當中的影響，並將之與傳 統脈衝式多重回波光達做比較，見圖 2-16，比較了相當多由波形萃取而得的參數。從 圖上可以看到，相較於脈衝式光達，使用全波形光達資料有助於提高地物分類的正確性， 但是也可以觀察到不同種類的擬合函數對於正確率反而降低。因此，此研究對於全波形 光達參數於地物分類的幫助是肯定的，但參數並不一定是愈多成果愈好；不全然有關係 的參數若是一起進到分類器中反而容易造成反效果，不同測區不同地物其最適用的參數 也不盡相同，使用者需要慎選。

15 圖 2-16 全波形與脈衝式光達對於地物分類之精度評估(Mallet et al., 2011)

2.1.4 支持向量機(SVM)、隨機森林(Random Forest)分類器應用

在過去的研究中使用了許多分類器進行地物分類，如類神經網路、人工智慧、支持 向量機(SVM)、隨機森林(Random Forest, RF)等等，分類器的發展是為了提升分類的成 果。而SVM因具有良好的分類類度，廣泛使用於遙測資料的地物分類(Huang et al., 2002; Pal and Mather, 2005), 及全波形光達於都市區域分類(Mallet et al., 2008)。

Random Forest (隨機森林)為Breiman在 2001 年所提出的演算法，是一種以決策樹為 基礎的整合式演算法，此方法擁有良好的資料處理能力。隨機森林的其中一個重要的優 點就是可以分析每一個分類參數所佔整個分類成果的重要性資訊 (Breiman, 2001)。而此 方法可以同時處理大量的變數，分類速度與效能都很好，適合用於多光譜的資料或是多 來源的資料(Gislason et al., 2006; Pal, 2005)。

有研究指出相較於SVM，隨機森林對於都市區域的分類成果精度來的更好，其可使 用光譜資料、粗糙度、全波形光達的回波波寬振幅、背向散射參數、形狀參數等眾多的 分類參數，顯示用隨機森林比SVM分類總體精度來的好，且其精度可達 95.75%。(Chehata et al., 2009)。 振幅 波寬、背向散射 形狀 擬合函數種類、不對稱性

16

2.2 全波形光達儀器介紹

本節對於目前商用空載陸域小足跡全波形光達儀器進行簡單的介紹，包含其功能特 色與基本參數，依照廠商分別有Leica ALS系列、Riegl LMS系列、Optech與Trimble。表 2-2 整理各全波形光達儀器之參數，包含其掃描方式、飛行高度、發射脈衝的波長、寬 度、能量、頻率、雷射光束之發散角、掃描角、足跡大小與其波形紀錄的取樣間隔。

表 2-2 全波形光達儀器參數

ALS60 LMSQ560 LMSQ680i ALTM3100 Harrier56

掃描方式 旋轉鏡 多邊形 多邊形 旋轉鏡 多邊形 飛行高度(m) 200~5000 <1500 <1600 <2500 <1000 發射波長(nm) 1064 1550 1550 1064 1550 脈衝寬度(ns) 5 4 3 8 < 4 脈衝能量(μJ) <0.2(mJ) 8 <200 脈衝頻率(kHz) <50 <100 <400 <200 <240 發散角(mrad) 0.22 0.3 0.5 0.3or0.8 0.3or1.0

掃描角 75 22.5 30 25 45or60

足跡(m@1km) 0.22 0.5 0.5 0.3or0.8 0.6

取樣間隔(ns) 1 1 1 1 1

17

2.2.1 Leica ALS系列

ALS60 為Leica公司所推出的全波形光達儀器，由掃描發射器、紀錄器與操作控制 器三個機器所組成，圖 2-17。其有高速的掃描旋轉鏡，因此在FOV為 75 度下，航高可 達 5000m，紀錄四個回波，可搭載中相幅相機(Leica RCD105)，同步拍攝影像，其影像 大小為 1280x1024 pixels (Leica Geosystems, 2012a)。

Leica公司目前最新推出的全波形光達儀器型號為ALS70，其有三種子機型：CM、 HP、HA(差別在飛行航高)，與ALS60 的差別是紀錄之回波數無限制且其掛載之相機 (Leica RCD30)可拍攝R G B NIR四個波段 (Leica Geosystems, 2012b)。

圖 2-17 Leica ALS60

2.2.2 Riegl LMS系列

LMS-Q560、Q680i為Riegl所推出的眾多全波形光達儀器中的兩種(圖 2-18)，許多全 波形光達研究文獻使用LMS-Q560 及Q680i，兩者皆使用旋轉多邊形掃描(圖 2-19)，Q680i 較不同的地方是脈衝頻率較Q560 高，掃描角度也較大(RIEGL, 2012a; RIEGL, 2012b)。

18 圖 2-18 旋轉多邊形掃描原理 圖 2-19 LMS-Q560(左)、Q680i(右)

2.2.3 Optech ALTM3100

ALTM3100 與前述儀器不同的地方在於其飛行高度較高，而發射之脈衝寬度也較其 它來的長，這會影響到其對地的解析度大小(Optech, 2012)，其儀器如下圖 2-20。 圖 2-20 Optech ALTM3100

19

2.2.4 Trimble Harrier56

Trimble(TopoSys）的Harrier56 也是一可以記錄全波形的光達儀器，其掃描原理與Riegl 系列相同，都是旋轉多邊形，其儀器如下圖 2-21 (Trimble, 2012)。

20

2.3 光達資料格式介紹

以下就本研究所使用的兩種全波形光達資料格式進行簡單介紹，分別為美國攝影測 量與遙感學會(American Society of Photogrammetry and Remote Sensing, ASPRS)所發布 的LAS1.3 格式(包含原始波形資訊)，和Riegl公司所使用的資料記錄格式SDW(未含原始 波形資訊)。

2.3.1 LAS1.3 格式

LAS格式是空載光達系統的標準資料格式，由 2003 年美國攝影測量與遙感學會所 發布的LAS 1.0 為標準格式，其後還有相繼發布 1.1、1.2、1.3、1.4 與 2.0 版。LAS文件 以二進位方式儲存，1.0 版的文件是由(1)公開檔頭區（public header block）、(2)變異長 度紀錄區（variable length records）、(3)點資料記錄區（point data records）三個部分所 組成，如圖 2-22。 公開檔頭區包含了光達資料的ID、文件產生的方式、飛行紀錄日期、坐標系統、坐 標尺度偏移因子、點雲數量、最大最小值等等基本的資訊，見圖 2-23。變異長度紀錄 區則是用來記錄使用者自行定義的訊息，長度不限。而點資料記錄區就是最重要的資訊， 內含光達每個點雲的三維坐標資料、回波反射強度值、回波電壓值與點的類別。 2009 年 7 月 14 日，最新的LAS標準格式 1.3 版正式發布，這也是本研究所使用的 全波形光達資料格式。其與LAS 1.0 最大的不同點在於，從原本資料格式的三個部份多 加了波形數據紀錄的區域(waveform package)，如圖 2-24。此外還有一些細節的增加與 缺點的修改，例如將點類別的部分進行擴充、點雲顏色紀錄、波形數據的參數資訊等等。

21

public header block public header block variable length records variable length records point data records point data records waveform package

圖 2-22 LAS 1.0 與LAS 1.3 圖 2-23 LAS公開檔頭資訊範例 圖 2-24 波形數據紀錄waveform package範例 LAS 1.3 的波形資料紀錄中，除了記錄三維坐標值與波形的電壓值之外，並記錄坐 標參數 (X(t), Y(t), Z(t))，使用者可自行計算出記錄訊號所對應的三維坐標，公式如下 (ASPRS, 2009)。 為初始坐標值，再加上LAS 1.3 的波形資料紀錄的坐標參數乘 以時間間隔t(單位為pico sec)，就可以得到不同時間點的三維坐標(X, Y, Z)，公式(1)如 下。

22 (1)

2.3.2 SDW格式

除了LAS1.3 格式記錄原始波形之外，Riegl公司使用SDF格式儲存波形，但是該格 式不是公開格式。另外Riegl自行制訂一個紀錄處理後波形資訊與參數的檔案格式，其為 一個二進位的檔案，檔案當中紀錄了Range, X, Y, Z, amplitude, echo width, return number, ClassID等資訊(見圖 2-25)，其檔案可以使用Riegl公司的SDCView軟體開啟，其檔案中 所提供的波寬振幅是使用高斯分解所求得(Riegl, 2010)。 由於資料取得的因素，本研究的實驗測區當中所使用的Riegl 公司的全波形光達資 料當中，測區二僅有SDW檔，所以採取直接應用其檔案中所提供之波寬振幅資訊，並未 自行對於波形進行分析處理。 圖 2-25 SDW檔案內容範例

23

第3章 研究方法

3.1 波形分析

本項研究工作目的為比較不同波形分析方法之成果。主要研究工作包含以下幾個部 分：資料前處理、計算初始值、擬合函數計算、檢驗分析。各項工作內容詳述如下，其 研究流程如圖 3-1 所示。 圖 3-1 波形分析研究流程

24

3.1.1 資料前處理

資料前處理目的是雜訊濾除，分為高斯平滑及濾除背景雜訊兩部份。高斯平滑使用 的是高斯函數濾波，目的是降低局部雜訊的干擾，高斯平滑化必須給定罩窗大小及平滑 因子σ兩個參數，當σ值越小的情況下，權值分布會向中央集中，濾波器形狀高且瘦其 平滑效果不彰，因受到平滑因子σ的影響，罩窗範圍開很大，罩窗邊緣的權值會非常接 近零，所以需要設定門檻。圖 3-2 紅色線為全波形光達原始波形，藍色為經過高斯平滑 後的成果，可以看到經過平滑之後的波形在峰值的部分會較原始資料來的低。在此資料 中，平滑因子設定為 5，罩窗大小為 9。 濾除背景雜訊是針對背景雜訊(Background Noise)進行處理，背景雜訊是指反射值 小於光達系統訊號雜訊比(Signal-to-Noise ratio)的雜訊值，因此以訊雜訊比設定門檻，若 反射值小於門檻則去除。圖 3-3 綠色為濾除背景雜訊前後的成果，正常而言在沒有接收 到回波訊號時值應為 0。在此資料中，背景雜訊的大小為 13。 圖 3-2 高斯平滑前後成果

25 圖 3-3 濾除背景雜訊前後成果

3.1.2 計算初始值

由於擬合波形的函數為非線性函數，所以需使用迭代方式求解，因此求解時必需 給 參數 初始值。由經 過 前處理 的 原 始資料 ， 利用 斜率變化 蒐 尋局部最大 值 (Local Maximum)，加上門檻設定以取得初始峰值位置，再依初始峰值的位置和其所對應之波 寬當作是迭代求解之初始值。由這種方法所找到的初始值運算效率佳，可做為高斯函數 求解的良好初始值。而韋伯函數需要比較精確的初始值，所以使用高斯函數所擬合出來 的波形參數作為韋伯函數求解初始值。 圖 3-4 是直接使用斜率變化從平滑後的波形當中所找出來初步初始值，可以看到紅 色箭頭的部分是一些錯誤的點，有可能是因為平滑的成果不好所致。所以加上門檻值設 定(例如水平與垂直距離太接近予以刪除)，刪除這些錯誤的點，在進到波形擬合解算， 就可以求得擬合成果，如圖 3-5 所示，可以萃取到三個回波。

26 圖 3-4 有錯誤的初始峰值 圖 3-5 使用門檻濾除錯誤的初始峰值與擬合成果 Echo 1 Echo 2 Echo 3

27

3.1.3 擬合函數計算

光達的回波波形與地物形狀相關，平坦地的回波寬度小，地物粗糙或是傾斜時的回 波寬度大且較不對稱，因此在波形擬合時，分別使用對稱及非對稱函數進行比較分析。 擬合函數計算分成以下三個部分，分別是(1) 波形分解公式、光達方程式與背向散 射參數推導、(2)對稱與非對稱函數、(3)Trust Region演算法。

3.1.3.1 波形分解、光達方程式與背向散射參數

本研究所使用的高斯分解法(Gaussian decomposition)其精神是將回波的波形，拆解 成數個不同的高斯函數，如圖 3-6 就有 4 個回波；同理非對稱的函數也是將波形分解成 數個不同的韋伯函數。其數學式如下， 是擬合之波形， 為數學函式， 為雜訊 值。 ∑ (2) 圖 3-6 波形分解圖例 時間(ns) 回波訊號 Echo 1 Echo 2 Echo 3 Echo 4

28 首先由雷達方程式進行光達方程式的公式推導，用於說明波形分解及背向散射參數 之求解方法。光達方程式是由雷達方程式改寫而得(Wagner et al., 2006)，其原理如圖 3-7， 由載具上的儀器發射出一個能量 ，經過雷射經距離 之後到地面，經過反射之後的能 量 再經由儀器接收。 能量遇到阻礙時，由於表面很粗糙，所以會發生反射及散射的現象，圖 3-8。其散 射的能量取決於地表的反射係數 與接觸的面積 。 圖 3-7 雷達原理示意圖(Wagner et al., 2006) 圖 3-8 背向散射示意圖(Wagner et al., 2006)

29

由以上的敘述可以推得下面的公式：式(3)為雷射足跡於距離R的面積大小；式(4) 為雷射到地面後的散射能量密度(雷射訊號發射功率除以雷射足跡的面積)；式(5)為雷射 訊號散射功率，與地表接觸面積和反射率有關；式(6)為雷射訊號接收功率，與接收能量 密度和接收孔徑大小有關；式(7)為最後整理完成的雷達方程式(radar equation)；式(8)為 移項後的背向散射參數(Backscatter cross-section coefficient)。

₍₃₎ (4) (5) ₍₆₎ (7) (8) 其中， ：雷射足跡面積 ：雷射發射器到物表的距離 ：光速發散寬度 ：散射能量密度 ：雷射訊號發射功率 ：雷射訊號散射功率 ：反射率 ：接收能量密度 ：反射圓錐角 ：雷射訊號接收功率 ：接收器孔徑

30 而光達方程式是由雷達方程式加上時間的延遲與衰減改寫而成式(9)如下，時間衰減 。當在range一定的範圍內[ ]，可以改寫成積分式(10)如下。當 時，方程式中的 可以提出積分之外，變成 與 的摺積(convolution)式(11)如下。 最後，接收到的能量就會等於N個回波能量疊加，如下，此方程式亦說明把疊置的回波 訊號分解成個別的回波相加的合理性。 ( ) (9) ∫ ( ) (10) ∫ ( ) (11) ∑ (12) ：時間 ：雷射脈衝在大氣中的速率 ：摺積運算子 而式(14)為背向散射參數 (Backscatter cross-section)與回波振幅、回波波寬的關係式， 由式(13)所推導而得(Wagner et al., 2006)。式(13)與式(12)不同的地方在於考慮系統與大 氣對能量衰減的引響。而式(14)當中的 為一校正常數，移項之後整理得式(15)。因此 配合距離、回波振幅、回波波寬、反射率、光速發散寬度就可以將每個回波的背向散射 參數求出來。 在計算 時，嚴密的作法是在地面設標量測實際不同地物的反射率，而在這邊我 們引用文獻，在特定條件與波長下(雷射波長 1550nm，光束直徑 0.15m，距離 1-1.5m)， 柏油路的反射率經由量測而得等於 0.25 (Alexander et al., 2010; Briese et al., 2008) 。因此

31 本研究簡化程序，假設同一航帶大氣條件一樣，使用人工選取柏油路面點求出 ，再 將其應用在整個航帶上。 ̂ (13) ̂ (14) , 其中 ̂ (15) 整合式(8)及式(15) (16) 柏油路 (17) ：系統轉換因子 ：大氣轉換因子 ：振幅 ：標準差 ：校正常數 ：回波振幅 ：回波波寬

3.1.3.2 對稱與非對稱函數

式(18)及(19)分別為高斯分佈數學式與韋伯分佈數學式，藉由此函數可以求得所需 的波形參數，該函數示意圖如圖 3-9。圖 3-10 是使用兩種不同函數進行擬合的示意圖， 此原始波形光達資料有些微的不對稱，使用不對稱的函數對擬合的成果將會有所助益。

32 ( ) ₍₁₈₎ ：時間(nanosec) ：接收訊號(volt) ：x軸位移 ：振幅(amplitude) σ：波寬(echo width) ( ) ( ) ₍₁₉₎ ：時間(nanosec) ：接收訊號(volt) ：振幅(amplitude) κ：尺度參數(scale parameter) λ：形狀參數(shape parameter) 圖 3-9 高斯函數 (左) 與韋伯函數 (右)

33

圖 3-10 對稱與不對稱函數之差異

3.1.3.3 Trust Region algorithm

因為擬合函數為非線性，求解需要進行迭代演算。本研究使用Trust-Region法來求 解非線性函數的未知參數，其方法的主要流程是在區域內找尋極值當作初始值，再移動 到另一區域反覆疊代計算求出最佳解，其主要的數學型式如下，每一次的疊代都要使 愈小(Lin et al., 2008) 。 ₍₂₀₎

f

:目標函數 k

w

:疊代次數 k  :相信區間 k q :二次方程

s

:迭代次數 R eturn S ignal Time [ns]

34

3.1.4 精度分析

波形擬合分析的檢驗分析方法有以下三個指標，分別為SSE、RMSE與。R-squared (R2)。SSE的涵義為殘差的平方合，RMSE為SSE除以樣本數再開根號，這兩個指標都是 愈小愈接近 0 成果愈好。R-squared (R2 )則是統計中預測觀測值好壞的指標，其範圍介在 0 到 1 之間，其值愈接近 1 代表成果愈好。本研究使用這三個指標來決定對稱與不對稱 函數的擬合的成果，各項參數如式(21)至(23)。 殘差的平方和 ∑ ̂ (21) 均方根誤差 √ √∑ ̂ (22) ∑ ̅ (23) : 觀測資料 ̂: 參考資料(函數擬合而得的估計值) : 觀測資料數目 ̅: 觀測資料平均值

35

3.2 地物分類

地物分類使用的全波形光達參數包含三維坐標幾何資訊、高程差、波寬、振幅、回 波數目，這些分類特徵內插成為網格式的資料，進入到分類器中進行分類。選用的分類 法分別有隨機森林與SVM兩種，最後進行精度評估與全波形資料探討，流程圖如圖 3-11。 圖 3-11 地物分類流程圖

3.2.1 分類特徵參數

本研究所使用的分類特徵參數皆來自於光達資料，分別有波寬(echo width)、振幅 (amplitude)、背向散射參數(backscatter cross-section coefficient)、高程值(Z)、高程差(dZ)、 回波數(echo number)、多重回波百分比(echo ratio) 。每個特徵參數萃取出來之後，會經 由內插產生網格式資料圖層，再進行分類。

36

據式(14)，使用人工選取道路點求出 ，再將其應用在整個航帶上，計算出每個回波 的背向散射參數(backscatter cross-section coefficient)。

dZ為一個代表地形粗糙度的參數，其計算方式為在一個固定範圍的面積當中，找尋 光達點雲高程最大與最小之差值(圖 3-12)。由於本研究所測試的區域地形高程起伏不大， 所以不考慮地形的影響，若是在地形起伏大的地方，必須使用的是對地形面鉛垂的高度 差值，如圖 3-13。 echo number則是點雲的回波編號，可以讓我們了解回波的先後順序與地物的關係。 多重回波百分比(echo ratio) 是一個計算多重回波比率的參數，此數值愈高就代表 此區域有較多回波地貌較為複雜，公式(24)如下， 為第一回波的數量， 為 多重回波中非第一也最後回波的數量， 為最後回波， 為單一回波(Höfle et al., 2008)。 圖 3-12 地形粗糙度的參數dZ ` 圖 3-13 考慮地形影響的dZ ( ⁄ ) (24) dZ dZ

37

3.2.2 分類方法

本研究使用兩種分類方法，一是隨機森林分類，二是支持式向量機(SVM)。分類的 一般作法是將訓練資料(由已知類別標記的資料所組成)建立分類模式，分類模式對未知 類別的測試資料進行預測。 隨機森林法是以決策樹分類為基礎，將許多不同決策樹的成果進行投票，選擇出最 佳的分類決策樹。SVM則廣泛使用於影像分類領域，是以統計學習理論提出的一種機器 學習方法，主要特性是在特徵空間中尋求具最大邊界的區分Hyper-plane以區分不同的兩 個類別。

3.2.2.1 隨機森林分類(Random Forests)

隨機森林分類法是一個基於決策樹分類的整合式分類法，而決策樹分類問題及答案 事實上可以用一個包含節點及方向箭頭，且具層級式結構的決策樹來完成，包含根節點、 內部節點、葉節點(圖 3-14)。 圖 3-14 決策樹分類(Tan et al., 2006) 決策樹之特性：決策樹分類的類別及其他屬性並不需要滿足任何機率分布；通常是 使用經驗法則找最佳的決策樹，所以很多決策樹都在大量的假設空間中進行搜尋。在資

38 料量大的情形下建立決策樹不難，其執行的速度也迅速，一旦決策樹建立後，針對測試 資料的分類也會變得非常快。尤其是對較小的樹而言，解釋上相對容易，其正確性也可 以與其他分類法進行比較。決策樹演算法可以處理雜訊問題，尤其是可以避免「過度學 習」的情形。重複的屬性不會影響決策樹的正確性，但是太多重複或是無關的屬性會造 成樹太大的問題，需要刪除修剪。因為大部分的決策樹都是採取由上而下及遞迴分割的 方法來處理，所以在葉節點中，也許會因為資料量太少而無法達到統計的顯著性，所以 可以設定當資料量少於門檻就不能再分割。子樹可以在決策樹中重複多次，這會使得決 策樹變得更為複雜，也許會更不易解釋。 隨機森林是一個特別設計給決策樹分類的整合分類方法，其結合多個決策樹的預測 結果，其由Breiman (2001)所提出，此分類方法的優點為可提供分類特徵資料的重要性， 並且其使用分類成果最佳的方針建立模型，在這樣的前提之下，其分類的精度會再一定 的標準以上。 圖 3-15 隨機森林分類步驟示意圖(Guo et al., 2011) 隨機森林法詳細演算步驟如下，隨機森林中的每棵樹都是根據隨機向量值所建立的， 而隨機向量是依據固定機率分配所產生。流程圖為圖 3-15： (1)決定訓練區資料n組，分類特徵參數m種。

39 (2)於m個輸入分類特徵參數中，選出T個子集合，T<m。 (3)運用bootstrap抽樣，從訓練資料重複抽樣形成InBag。 (4)對每一個InBag，隨機選擇m個分類特徵參數，根據這些分類特徵參數，計算決策樹， 最後進行投票(vote)，選擇出最佳的成果。 (5)注意的是每個決策樹都不需要進行修剪、完全生長。 隨機向量的決定方法：每個節點隨機選取F輸入特徵來進行分割，所要分割的節點 是由選取的F特徵中決定出來，數目持續生長不需修改。如果原始特徵個數太少，很難 隨機選出互相獨立的輸入特徵，所以就要增加特徵空間，從這些隨機結合的新特徵F中 選出最適合分割的節點進行之。在每個節點上隨機從F個最好的分割點中選取一個來產 生隨機樹，演算法必須檢查所有節點上的分割特徵。 隨機森林法需要控制的參數只有兩種：(1)決策樹的數量、(2)進入每個決策樹中的 特徵參數數量。在bootstrap抽樣的過程中，沒有被抽中的樣本為out-of-bag(OOB)，可以 用來估計模型的錯誤率(OOB-error)。 理論上，當樹夠大時，隨機森林可以保證所產生的推論錯誤率之上限公式(25)如下 : 推論錯誤率 ̅ (25) ̅ ：表示樹間的平均相關程度 ：表示樹分類法的分類能力

3.2.2.1 Support Vector Machine (SVM)分類

Support Vector Machine (SVM) (Vapmik (1995),Burges (1998))在航遙測的領域中，被 廣泛利用在影像分類的部分。起源於統計學習理論，已成功用在一些應用問題中，像是 手寫辨識或是文件分類，也可以處理高維度的資料。

其分類法主要的精神是對於一群在特徵空間中的資料，希望能夠在該空間之中找出 一Hyper-plane，並且希望此 Hyper-plane 可以將這群資料切成兩群(例如:群組A、群組

40 B) 。 而 屬 於 群 組 A 的 資 料 均 位 於 Hyper-plane 的 同 側 ， 而 群 組 B 的 資 料 均 位 於 Hyper-plane的另一側，見圖 3-16。而為了可以群組AB明確地分辨出來，所以兩個群組 間的邊界(Margin)愈大愈好分離。 一般可以線性分割的函數，其公式如式(26)及(27)，而處理非線性的資料，線性函 數無法分割，必須進行屬性轉換，把資料轉至更高維度的空間或是特徵空間，才能加以 分辨出來(圖 3-17)，其公式如式(28)及(29)。 而通常轉換函數 是一個複雜且不容易求 得的函數，所以我們會將轉換函數 做內積，得到要一個較為簡單的函數，稱作核函數 (kernel function)。 圖 3-16 SVM原理示意圖(線性分割) 圖 3-17 非線性屬性轉換示意圖

41 線性分割： ‖ ‖ (26) (27) 非線性要進行屬性轉換 ： ‖ ‖ (28) (29) 本研究使用三個不同的kernel 進行分類成果的討論：(1) Polynomial Kernel 、(2) a universal kernel function based on the Pearson VII function (PUK Kernel) 、(3) Radial Basis Function (RBF) Kernel。第一個Polynomial Kernel為一個基礎的多項式核函數，而第二 個PUK Kernel是一個通用的核函數，其具有大的強鈍性(Robustness)，具有很好的屬性 轉換能力，且其成果跟其他核函數相較並不差或是更好，因此常用來替代一般常見的 多項式核函數或是RBF函數(Ü stün et al., 2006)。這三個核函數皆較基礎且較多人使用， 也被常用在SVM分類中做為討論比較的項目(Pal, 2009)。

42

3.2.3 精度評估

此部分的精度評估提供兩項成果，一為使用參考地真資料所產生的正確率、錯誤率、 Kappa值、及RMSE指標。二為隨機樹分類中各類別的重要性評估統計圖。RMSE均方根 誤差同波形分析之精度評估中的式(22)。正確率與錯誤率，公式如下： 正確率 正確預測的個數 總預測個數 (30) 錯誤率 錯誤預測的個數 總預測個數 (31) Kappa值之公式如下： 總體準確度 期望準確度 期望準確度 ∑ ∑ ∑ (32) n:類別 N:檢核點數 :誤差矩陣對角線元素 :列總合 :行總合

43

第4章 研究資料

4.1 測區一

本資料的飛行時間為 2010 年 8 月 2 日，由Leica ALS60 系統所蒐集得，記錄成LAS1.3 格式，地面位置在台灣中部山區(中央山脈)，地物多為樹林。總共有 4085611 筆資料： 第一回波 3568012 筆、第二回波 488833 筆、第三回波 28265 筆、第四回波 501 筆。下 圖 4-1 左邊是測區的全圖，右上則是本研究所使用的測試區域(70m*50m)，右下為其光 達點雲的地表高程分布圖，其高程在兩千一百公尺左右，總點數(peak)有 1981 個點，其 中單一回訊佔了很大的部分，可能的原因是樹林太茂密，導致光達的穿透率不佳，難以 獲取地面點。 本測區資料提供接收器所接收的原始波形電壓值(DN值)之LAS1.3 檔，因此可以對 此進行波形的分析，觀察其對稱與不對稱函數之擬合成果。但因其地表物較難分辨，地 真資料蒐集困難，不易判識，因此並未對此測區進行地物分類。

4.2 測區二

本資料位於成功大學附近之台南都市區，掃描時間是 2010 年 8 月 4 日，由RIEGL LMS-Q680i系統所蒐集，總共有 21553696 筆資料：第一回波 15830714 筆、第二回波 2134079 筆、第三回波 399615 筆、第四回波 56772 筆、第五回波 5268 筆、第六回波 379 筆、第七回波 27 筆、第八回波 11 筆。其記錄的資料格式為SDW與SDC，紀錄內容有 Range, X, Y, Z, amplitude, echo width, return number, Class ID等，但是缺少波形的原始資 料，因此本測區在研究中主要是用來進行地物分類，用資料中所提供的波寬與振幅加上 自行計算的背向散射參數，加以比較其加入波形參數的重要性。

台南成大校區，地物大部分為建築物與人造路面。下圖 4-2 左為全區的光達點雲強 度圖，右上為測試區(160m*100m)影像圖，右下為測試區之光達點雲高程分布圖。圖 4-3

44

則為測區二之參考地物類別資料，總共分為七個類別：樹木、柏油路、高起牆、草地、 水體、建物、人工地。

45

46 樹木 柏油路 高起牆 草地 水體 建物 人工地 圖 4-3 測區二之參考地物類別資料

4.3 測區三

本測區位於桃園的龍潭，其地表覆蓋物型態為半都市半林區，掃描時間是 2010 年 8 月 4 日，由RIEGL LMS-Q680i系統所蒐集，總共有 71667219 筆資料：第一回波 49948499 筆、第二回波 13835419 筆、第三回波 5362194 筆、第四回波 1853375 筆、第五回波 522313 筆、第六回波 145419 筆。其記錄檔有LAS1.3 格式的波形數據原始檔與提供Range, X, Y, Z, amplitude, echo width, return number, Class ID等的SDW檔。

下圖 4-4 左為本測區的全圖，可以看到航帶成細長型，右上為測試區(60m*40m)影 像圖，右下為測試區之光達點雲高程分布圖。圖 4-5 則為測區三之參考地物類別資料， 總共分為五個類別：樹木 1、樹木 2、裸露地、草地、建物(樹木 1 與樹木 2 為不同的樹 種)。

47 圖 4-4 測區三點雲影像及四角坐標 圖 4-5 測區三之參考地物類別資料 樹木 1 樹木 2 建物 裸露地 草地

48 表 4-1 研究資料與實驗內容 為本研究的資料與實驗內容整理表，因為不同測區的 地物型態不同、其原始的資料類型不同，所以各個測區之實驗內容探討也根據資料的差 異而變化。測區一因其地物類型單一，參考地物分類不易取得，所以並未做地物分類； 而測區二的資料並沒有波形資料的原始檔，只有經過處理後的波寬振幅，所以只有使用 來做地物分類的探討；而測區三的資料類型較為完善，所以兩個實驗都有使用此測區的 資料來進行。 表 4-1 研究資料與實驗內容 測區一 測區二 測區三 地點 中部南投山區 台南成大校區 桃園龍潭 地物型態 林區 都市區 半林區半都市

掃描儀器 Leica ALS60 RIEGL LMS-Q680i

實驗內容 波形分析 〇 X 〇

地物分類 X 〇 〇

資料類型 LAS1.3 SDW LAS1.3&SDW

49

第5章 成果分析

本研究的成果分析分為三個部分，第一個是波形模擬，其主要是在觀察與探討不同 狀況下兩個回波的分辨力；第二部分為波形分析，主要為觀察波形與地表物之關係、波 形擬合成果與展示波形參數圖；最後部分則是地物分類的成果，探討不同分類方法以及 分類法中訓練區與參數設定的成果與比較。

5.1 波形模擬

因為在做波形擬合的時候，必須提供初始值來進行，所以波形模擬的目的是在不同 狀況下，量化兩個回波的分辨力，藉此來觀察理想狀況下波形的特性，提供一個合理的 參考標準。此模擬設計不同的振幅、雜訊比下，觀察兩個波形的分辨力與波寬之間的關 係。 進行模擬波形步驟：(1)定義參數，使用高斯分佈模擬兩個波、(2)將兩個波疊加， 並加入不同大小的雜訊、(3)將前面步驟所模擬出來的波形，使用本研究波形分析的前處 理流程，觀察統計初始否可以分辨出兩個相靠近的波。

5.1.1 模擬變數設定

以下三項為波形模擬的變數設定，分別為振幅比例、波寬大小與雜訊大小。舉例如 圖 5-1，振幅比例 0.8(80/100)，波寬 8ns，兩個波峰距離 2.5 倍波寬，雜訊大小 1。圖 5-2 為模擬波形分辨力示意圖，分辨力為兩個波峰可被找出的最接近值。 (1)振幅比例：分為 1.0、0.8、0.6、0.4、0.2 倍五項。 (2)波寬大小：分為 4、8、16、32 (ns)，4ns為目前多數空載全波形光達的系統設定。 (3)雜訊大小：noise deviation大小為 1、2、3，平均值根據Signal-to-Noise ratio給定為 5，

noise deviation的給定標準為統計本研究測區三的波形資料之統計標準差 而得。

50 圖 5-1 模擬波形示意圖 圖 5-2 模擬波形分辨力示意圖(左：可以分辨、右：無法分辨) 0 10 20 30 40 50 60 70 80 90 100 0 20 40 60 80 100 120 Time [ns] R e tu rn S ig n a l Waveform Original Signal Reconstructed Signal Initial Value Peak Extracted 0 10 20 30 40 50 60 70 80 90 100 0 20 40 60 80 100 120 140 Time [ns] R e tu rn S ig n a l Waveform Original Signal Reconstructed Signal Initial Value Peak Extracted w Re tur n Signa l Time [ns]

51

5.1.2 模擬成果分析與討論

表 5 1 為本實驗的成果，表格中的數值為分辨力，其單位為波寬的倍數；表 5 2 是 將表 5 1 的比值換算成地面的距離，這樣較為直觀。當noise deviation為 1、波寬 4ns、 振幅比例 1 時，分辨力為 1.101 倍波寬。代表的就是當兩個一樣大的波形，在彼此距離 靠近到小於 4.404ns時，是無法被分解出來，而換算成地面距離為 0.66 公尺，所以表示 在地面高度差異 0.66 公尺以下的物體沒有辦法被分辨出來。 由實驗的成果顯示：兩個波的振幅相差愈大(一大一小)，其愈難分辨。波寬愈大的時 候，其分辨力的值會愈小，但是乘上波寬大小換算成地距之後，其解析度是較差的。當 雜訊愈小的時候，分辨力的成果會愈一致，雜訊大的成果趨勢線則會較分散，可從圖 5-3、 圖 5-4、圖 5-5 觀察出來。 本實驗的成果顯示一般接近目前光達系統設定的情況下：noise deviation =1、波寬 =4ns時，不同比例的振幅大小下其分辨力約為 1~1.6 倍波寬之間，換算成地面距離則是 60~100 公分。圖 5-6 為不同波寬與noise deviation互相影響之關係圖，可以觀察到當波 寬小的時候(4ns)，noise deviation對分辨力的影響量很大，所以對於波形資料來說，前處 理的雜訊濾除工作必須要做得好，才不會影響接下來的函數擬合。

52 表 5-1 模擬成果之分辨力大小(單位：與波寬之比值) 1 0.8 0.6 0.4 noise deviation=1 4ns 1.100 1.247 1.359 1.506 8ns 1.018 1.160 1.270 1.387 16ns 0.932 1.109 1.216 1.317 32ns 0.716 0.828 1.087 1.171 noise deviation=2 4ns 1.305 1.442 1.578 1.694 8ns 1.060 1.202 1.316 1.437 16ns 0.946 1.114 1.226 1.329 32ns 0.721 0.800 1.063 1.141 noise deviation=3 4ns 1.561 1.704 1.817 1.951 8ns 1.119 1.264 1.380 1.505 16ns 0.979 1.136 1.246 1.350 32ns 0.775 0.884 1.101 1.172 表 5-2 模擬成果之分辨力大小(單位：距離(m)) 1 0.8 0.6 0.4 noise deviation=1 4ns 0.66 0.75 0.82 0.90 8ns 1.22 1.39 1.52 1.66 16ns 2.24 2.66 2.92 3.16 32ns 3.44 3.97 5.22 5.62 noise deviation=2 4ns 0.78 0.87 0.95 1.02 8ns 1.27 1.44 1.58 1.72 16ns 2.27 2.67 2.94 3.19 32ns 3.46 3.84 5.10 5.48 noise deviation=3 4ns 0.94 1.02 1.09 1.17 8ns 1.34 1.52 1.66 1.81 16ns 2.35 2.73 2.99 3.24 32ns 3.72 4.24 5.28 5.63 波寬 振幅比例 波寬 振幅比例

53 圖 5-3 noise deviation=1 圖 5-4 noise deviation=2 圖 5-5 noise deviation=3 (ratio) (ratio) (ratio)

54 圖 5-6 不同波寬與noise deviation相互影響關係圖 (ratio) (ratio) (ratio) (ratio)

55

5.2 波形分析

波形分析分為以下四個部分：第一部分為波形與地物的關係，藉以展示不同地物類 別的波形；第二部分為對稱與不對稱函數擬合波形的成果比較；第三部分為成果圖層展 示；第四部份則是背向散射參數的改正。

5.2.1 波形與地物

當雷射光碰到物體的時候，會散射產生一個回波，所以當其接觸到一光滑平面的時 候只會有一個波形，如圖 5-7 可以看到左邊十字的中心為無植被覆蓋的裸露地，因此就 會產生右邊單一回波的波形。 但是當雷射光掃到不只一個物體的時候，就會產生多個回波，如圖 5-8 可以看到雷 射光穿透植被到達地面時會產生好幾個回波，就可以幫助我們了解地物的垂直分布。但 回波數量除了顯示地物的類型之外，跟光達能量的穿透性也有很大的關係，如果穿透性 不佳的話，還是沒有辦法獲得覆蓋在植物冠層以下的地面回波。 圖 5-7 裸露地與其回波波形 Ret u rn Sig n al Time [ns]

56 圖 5-8 波形與地物的關係(測區一) R etu rn Sig n al Time [ns] R etu rn Sig n al Time [ns] R etu rn Sig n al Time [ns]