應用空載光達資料估計森林樹冠高度模型及 葉面積指數

Volume19, No2, December 2014, pp. 107-123

1國立成功大學測量及空間資訊學系 碩士 收到日期:民國 103 年 03 月 21 日

2國立成功大學測量及空間資訊學系 博士 修改日期:民國 103 年 10 月 13 日

3國立成功大學測量及空間資訊學系 教授 接受日期:民國 103 年 10 月 30 日

4 國立成功大學測量及空間資訊學系 助理教授

＊通訊作者, 電話: 06-2383876 ext. 852, E-mail:s3350839@gmail.com

應用空載光達資料估計森林樹冠高度模型及 葉面積指數

林莉萍

王正楷

曾義星

朱宏杰

摘 要

空載光達系統可直接獲取高精度三維坐標的點雲資料，且其雷射光具有穿透樹葉縫隙的特性，能快 速偵測得森林林分結構的三維空間資訊。透過其蘊涵的森林厚度及林木密度資訊，可用來估計樹冠高度 模型(Canopy Height Model, CHM)和葉面積指數(Leaf Area Index, LAI)。本研究目的是探討如何應用空載光 達資料估計森林地區之樹冠高度模型和葉面積指數。在推估 CHM 方面，本文是以點雲資料產製數值表面 模型(Digital Surface Model, DSM)及數值高程模型(Digital Elevation Model, DEM)，將森林區域的 DSM 減 DEM 即得 CHM。並以三種不同的空載光達系統產生之原始點雲及全波型資料進行實驗，比較六個推估 得的 CHM，實驗結果發現，整體差異均小，少數差異較明顯的地方主要發生在森林表面高度落差較大處。

在 LAI 方面，則是以點雲資料計算五種雷射穿透率指標(Laser Penetration Index, LPI)來推估，LPI₁是地面 點占全部點的比例，LPI₂是地面點強度值和全部點強度值的比例，LPI₃是地面點與所有雷射光束的比值，

LPI₄是改良自 LPI₁，但增加單一回波地面點的權，LPI₅則是以全波形資料計算的雷射光照射地物間面積 與非地面的面積比。五種 LPI 與實際地面量測資料迴歸後的成果顯示，在本研究測區的低航高點雲資料 中，LPI₄都能得到高於 0.5 的 R²值，說明 LPI₄具有較穩定且準確的 LAI 估計能力。而相較於即時解算的 多重回波點雲，使用全波形的光達點雲穿透率指標，可提升對 LAI 的估計，R²可達到 0.8 以上，增加利 用小範圍地面實測資料來推估大範圍森林區資料的可行性。

關鍵詞：樹冠高度模型、葉面積指數、空載光達、雷射穿透率、全波型

1. 前言

森林中的資訊包含了林分高度、林木材積、葉 面積指數、樹種等等，不同的指標都代表著森林生 長的情形。因此利用觀測儀器探測森林的相關參數 指標，有利於我們了解在整個森林環境的變化。為 了能快速且能經常性監測森林生長，除了應用現地 量測的方式外，利用航遙測的技術進行大規模的森 林量測亦趨於普遍 (吳守從及陳永寬，2004) 。常 見的遙測有多光譜影像或可見光之航空影像，然其 缺點為其影像所記錄之觀測物體大多為樹冠層，樹 冠下的資訊則難以得知，而另一遙測技術─空載光 達系統，因其具有穿透樹冠層之能力以及可直接提

供三維資訊的優勢(Danson et al., 2009)，相比於傳 統遙測之多光譜或可見光之航空影像，空載光達之 使用已逐漸普遍應用在森林的測繪應用上，許多文 獻也根據空載光達資料發展出多種點雲處理方法 以推估森林的相關參數指標，例如偵測立木位置、

林木高度、分辨樹種、生物量及葉面積指數等 (Koetz et al., 2006)。

樹冠高度模型包含了森林的高度，以及整體的 森林體積、範圍等，可透過量測樹冠表層及地面的 三維資訊來得到。空載光達系統可提供高密度的掃 描點三維坐標，這些點位多數分佈於樹冠層，但約 有 30%穿透樹冠表層的點，分佈於森林內部或地表

面，經過點雲過濾篩選處理，可產製數值表面模型 (DSM)及數值高程模型(DEM)，將森林區域的 DSM 減 DEM 即得 CHM。本研究應用同地區多種空載 光達資料分別製作 CHM 並比較其差異，以評估應 用空載光達點雲製作 CHM 的穩定性及精度。

本研究另一主要目的是利用空載光達資料估 計森林葉面積指數(LAI)，森林 LAI 的分佈可展示 森林區域內的林木樹冠生長情形，而由於植冠密度 的分佈直接影響森林整體環境，如光合作用、蒸散 作用、呼吸作用等，都與整個地球的水資源循環、

碳含量的變化息息相關。藉由地面直接量測、地面 光學儀器、或光譜影像等方式都能得到森林中某地 點的葉面積指數(Jonckheere et al., 2004；Masona et al., 2012 )，但若需要全面性的觀測，則有賴航遙 測的作法。而由於空載光達具備直接獲取三維的點 資訊及穿透之特性，適合用來獲取森林樹冠及內部 之幾何結構。依據 LAI 的定義，LAI 與樹林葉片的 總面積成正比，而雷射光可穿透樹冠層的比例與葉 片總面積成反比，因此可透過此反比關係將空載光 達的樹林穿透率指標(Laser Penetration Index, LPI) 來估計 LAI。Morsdorf et al. (2008)、彭炳勳(2007)、

及 Solberg(2010)皆曾利用光達點雲的 LPI 來推估 針葉林的 LAI，其成果顯示 LPI 與 LAI 的相關性判 定係數 R²都能達到 0.7 以上。

台灣低海拔森林屬闊葉林，森林型態複雜且明 顯異於針葉林，值得進一步探討空載光達是否同樣 可以用來估計闊葉林的 LAI，本研究嘗試以多種 LPI 推估闊葉林的 LAI，並探討全波型光達對 LAI 的估計是否有幫助。

2. 研究資料及方法流程

2.1 實驗區域

實驗區域位於台灣東南部的南仁山生態保護 區，此區為熱帶季風型雨林，是台灣僅存的低海拔 原始林，其森林型態複雜，孕育非常豐富的生態景 觀。由於此地區地形複雜、樹冠遮蔽嚴重，造成定 位、取樣不易，因此我們實驗區的的選擇乃尋找其

鄰近位置具備可接收 GPS 衛星資料之區域，再利 用導線測量的方式將 GPS 獲取的控制點資訊導入 選定樣區。本研所選取了兩個小範圍區域如圖 1 所示。其中 A 區附近具備良好的空闊地點，可接 收 GPS 衛星訊號(以 e-GPS 方式進行)，因此此區 的實驗資料適合用來進行 LAI 之估計，而 B 區由 於 GPS 訊號仍不穩定，難以解算，因此缺少控制 資訊，無法計算對應的 LAI 資料，而對於 CHM 的 估計，兩區則都有進行。

對於兩區的森林型態，A 區為迎風面的位置，

整體樹高不超過 15 公尺，呈現樹幹較細且立木間 距較近的森林型態，整區其大小約為 300×250m²。 而 B 區位於背風面，其範圍大小約 250×250m²，該 區的樹冠生長茂密，整體的林木高度較 A 區來的 高，其中的優勢樹種的高度約在 15~20 公尺左右，

地形包含了中間低兩側高的溪谷型態。

2.2 空載光達資料

研究中共使用了三種空載光達系統，分別是 Leica ALS60、Riegl LMS-Q680i 和 Optech Pegasus HD400，為了能獲取高密度的地面點點雲來製作數 值高程模型，皆使用低航高(離地面高度 1000 公尺) 且光跡直徑小於 1 公尺的掃描方式進行森林區的 量測(表 1)。另外，在 Leica 系統的部分，則是增 加了高航高(離地面高度 2500 公尺)的資料來進行 估計葉面積指數的比較。其中 Leica ALS60 和 Optech Pegasus HD400 光達系統是利用加裝全波 形資料記錄器的方式，因此可同時提供即時解算的 多 重 回 波 光 達 點 雲 和 全 波 形 資 料 ； 而 Riegel LMS-Q680i 儀器則是提供經全波形資料後處理過 的三維點雲以及原始全波形資料。而本研究對於三 家儀器所提供的全波形資料，利用小波偵測器搭配 高斯擬合的方式，能夠再重新萃取出新的點雲，為 了有所區別，系統所給予的即時解算點雲我們稱為 系統點雲(System LiDAR Points)，而利用小波偵測 器 再 重 新 由 波 型 解 算 之 點 雲 稱 為 全 波 形 點 雲 (Waveform LiDAR Points)。

對於各種低航高光達系統的資料而言，本研究

使用了三個高重疊航帶的 Leica 光達資料，來提升 在單位面積上的點雲數，期望能得到更完整的森林 結構提高計算穿透率指標的準確性，進而提升葉面 積指數的估計。而在相同航帶的資料下，透過全波 形萃取出的點雲總數量可以增加約百分之三十左 右(與系統點雲比較)，地面點的數量也有明顯的增 加(表 2、表 3)； Riegl 的部分相較於多重回波，利 用小波方式萃取出的點雲數量增加的比例較 Leica 的資料少。其原因為 Riegl 原有的資料即是利用全

波形資料萃取得來，在單一條光束的回波訊號沒有 受到點數的限制，因此與本研究中使用的全波形點 雲數量較為一致；雖然 Optech 在單位雷射數量上 的地面點比例較高。但在本樣區只用了單一條航帶 的 Optech 點雲進行樹冠模型的製作，因此地面點 相較於其他兩種資料來的要少。而經全波形資料萃 取出的總點雲雖有些許的增加，但主要增加的部分 都是屬於地表上樹冠層的區域，地面的數量反而有 減少的現象。

(a)

(b) (c) (d)

圖 1 樣區範圍圖(a) 南仁山保護區，框選處為本研究的兩個實驗區；(b) A 區的現地實景圖；(c) B 區現地 實景圖；(d)A 區中取樣 LAI 的樣本分佈

表 1 光達系統掃描參數

Leica ALS60 Riegl LMS-Q680i Optech Pegasus HD400 AGL 1000m 2500 m 1000m 1000m

雷射脈衝頻率 99 kHz 99 kHz 210 kHz 150 kHz

FOV 28° 20° 60° 40°

波長 1064nm 1064nm 1550 nm 1064nm 發散角 0.22 mrad 0.22 mrad ≤ 0.5 mrad 0.25 mrad 光跡大小 0.22 m 0.55 m 0.5 m 0.25 m 掃描日期 2011/10/21

2011/10/24

2011/10/24

2012/01/08 2011/11/3 飛行速度 185 KM/H 185 KM/H 180 KM/H 203.5 KM/H

A B

表 2 A 區光達點雲資料 總 雷 射

光束

雷射密度 (pt/m2)

總點數 點雲密度 (pt/m2)

地面點點 數

地面點雲密 度(pt/m2) Leica 多重回波 547758 7.3 600802 8.0 106975 1.4 Leica 全波形光達 541059 7.2 697642 9.3 115498 1.5 Riegl 多重回波 491364 6.5 657740 8.7 128060 1.7 Riegl 全波形光達 491296 6.5 723564 9.6 121432 1.6 Optech 多重回波 105884 1.4 161215 2.1 38885 0.5 Optech 全波形光達 106040 1.4 163837 2.2 36511 0.5

表 3 B 區光達點雲資料 總雷射

光束

雷射密度 (pt/m²)

總點數 點雲密度 (pt/m²)

地面點點 數

地面點雲密 度(pt/m²) Leica 多重回波 875503 12.3 1027138 14.4 22949 0.32 Leica 全波形光達 875518 12.3 1337504 18.8 38209 0.54 Riegl 多重回波 638711 9.0 930788 13.1 30406 0.43 Riegl 全波形光達 638663 9.0 1026081 14.4 33904 0.48 Optech 多重回波 96269 1.4 168830 2.4 12701 0.18 Optech 全波形光達 96291 1.4 179382 2.5 12413 0.17

2.3 現地觀測資料

對於葉面積估計的現地觀測部分，我們一共 在 A 樣區內選擇了 10 個樣本點進行葉面積指數的 量測。為了確認各樣本點中心的位置，我們先在樣 區附近較空曠的地方架設 E-GPS，進行多次量測，

並利用 E-GPS 量測測區周圍一等水準點 Q029、

Q028 的坐標來進行修正，以得到較準確的控制點 坐標。利用測得的控制點坐標，我們以全站儀進行 閉合導線的量測，於測區中釘定臨時木樁，作為測 區內控制資訊之用，藉以獲取各樣本的中心位置。

在定位完成後，我們以 LAI-2000 進行葉面積指數 的量測。

LAI-2000 的儀器是利用樹冠透光的原理，以 魚眼鏡頭的光感應器水平向上接收樹冠底下的光 線，並同時以另一台儀器接收未受到遮蔽的全天光 進行光量的比較，藉由其中的亮度比值求得葉面積 指數。本研究量測各樣本點中心的方式是在距離樣 本中心點 0~1 公尺的位置，分別朝八個不同方向，

高度約 1 公尺處向上接收樹冠層下的光量(圖 2)。

為了降低量測時受到操作人員或其它鄰近物體之 遮蔽、以及陽光過度曝曬的影響，本研究中使用了 90∘的透鏡遮蓋輔助 LAI-2000 的量測。最後再將 接收的八個平均亮度與全天光下的亮度相比，即可 得到該樣本點的葉面積指數。

圖 2 LAI-2000 測量方式之俯視圖

2.4 樹冠高度模型估計

樹冠高度模型(Canopy Height Model, CHM)是 由數值表面模型(DSM)和數值高程模型(DEM)兩

種模型所相減而得，圖 3 展示森林區典型的光達點 雲剖面圖，以及所萃取得的 DEM 及 DSM 剖面。

CHM 可以用於描述森林區域樹木的形狀和範圍，

諸如林木高度、樹體積、林木位置等。由於森林的 地形複雜，人員不易進入測量該地區的地形變化，

且森林區有茂密的樹冠覆蓋，難以用遙測影像找尋 匹配點進行立體製圖，甚至無法觀測到地形起伏的 情況製作 DEM，因此目前多使用空載光達點雲的 資料來進行樹冠高度模型的製作。Næ sset(1997)以 光達點雲的最後回波值高度，加上胸高直徑來推算 樹的平均高度，最後得到的結果普遍有低估的情形。

其主要原因為當光達系統在進行掃描時，可能會錯 過樹木最高點的位置，或是因為最高點的樹梢回訊 較弱而無法取得，因此取到的樹高通常會比實際樹 高要來的低。

數值高程模型的精度也會影響森林體積的估

算，尤其是在地勢起伏較大的山區，若是能有較準 確的地形模型(DEM)，才能夠得到較好的樹冠高程 模型 (蕭淳伊，2008) 。因此我們在製作 CHM 之 前會針對 DEM 的部份，比較利用全波形光達點雲 和一般光達點雲製作成果的差異。而模型的網格大 小也是個考慮的重點，因為若單位網格內的點雲數 量過少會容易造成取樣不足有產生錯誤(Naesset, 1997)。在彭炳勳 (2007)研究中也顯示，樹冠模型 的網格大小和計算面積對評估林分平均高都會有 影響。在其測試下，以空間解析力為 3 公尺之網格，

搭配 15 公尺×15 公尺的計算面積可以得到最符合 實際樹高的成果。因此本研究利用了三種不同的光 達系統點雲資料，並配合點雲密度切割適合的網格 大小製作樹冠高度模型，來比較其中的差異和製作 之精度評估。

(a) (b)

(c)

圖 3 利用空載光達點雲製作樹冠高程模型示意圖(a)數值高程模型(DEM)；(b)數值表面模型(DSM)；(c)樹 冠高度模型(CHM)

2.5 森林葉面積指數估計

利用光達資料推估葉面積的方式，乃藉由點雲 的分佈，計算出能夠表示林分樹冠與地表間關係的 穿透率指標(Laser penetration index, LPI)。將 LPI 和已知的現地量測葉面積指數資料進行迴歸計算，

可得到 LPI 推估葉面積指數的模型參數，利用此參 數就能透過 LPI 來估計森林葉面積指數。其中，由 於無法確認地真資料量測時所涵蓋的區域，因此以 測試的方式，找出相對應穿透率指標推估葉面積指

數的計算範圍。最後再以最佳的迴歸參數，針對大 範圍的森林區進行葉面積推估。以下之各小節分別 針對點雲穿透率指標，敘述估計葉面積指數作業流 程(圖 4)中各個步驟的詳細做法。

本研究按照不同點雲的特性，並參考研究文獻 中推估葉面積或孔隙率的穿透率指標計算方法，統 整出五種葉面積的推估方法(Solberg et al., 2006;

Hopkinson& Laura, 2009；Arnaud et al., 2011)。其 中，為了要和地面葉面積指數觀測資料進行比較，

將 LAI-2000 地面測量的高度(一公尺)設為門檻值，

將此高度之下的點皆視為地面點雲。

圖 4 利用光達點雲推估葉面積指數流程圖

2.5.1 計算穿透率指標

第一種穿透率指標(LPI₁)為較直觀的算法，是 以地面點占全部點的點數來計算(式 1)。

LPI₁=_ΣNtotal^ΣNg (1) 其中的 Ng 為地面點點數；Ntotal 表示所有點點數。

為多數研究文獻採用的算法(Solberg et al., 2006;

Arnaud et al., 2011)，與植被覆蓋率也有較好的相關 性(黃清美，2007)；第二種穿透率指標(LPI2)也是計 算地面點和全部點的比例，但改以點雲的強度值為 比例的計算(式 2) (Solberg, 2008)。

LPI2=_ΣItotal^ΣIg (2) 式中，Ig=地面點強度值 ; Itotal=所有點強度值。

第三種穿透率指標(LPI3)是全部雷射光束的數 量(Nl)和地面點的點數比(式 3) (黃清美，2007；Zhao

& Popescu, 2009)

LPI3=^ΣNg_ΣNl (3) 第四種穿透率指標(LPI4)是改良自第一種指標，

對單一回波的地面點進行加權。此方法的想法來自 在許多文獻中(Solberg, 2010；Koetz B. et al., 2006；

Arnaud et al., 2011)，除了使用第一種穿透率指標外，

也使用了單一回波的點雲做穿透率指標的計算，以 全部單一回波的地面點和雷射數的關係推算葉面

積指數。由於雷射光經過樹冠層得到的多重回波訊 號較為複雜，難以確定得到的點數與茂密程度呈現 絕對的正比關係，因此本研究中將其合併使用，強 調單一回波地面點在森林中顯示的直接穿透性。

LPI₄如式(4)，其中以 Ng1 代表所有單一回波且為 地面點的點數。

LPI4=^Σ(Ng+Ng1)_ΣNtotal (4)

Solberg(2010) 利用點雲強度值和測試的物體 反射率比值，推算雷射光照射地物間面積的比例。

因此本研究依據面積比例的概念，提出第五種穿透 率指標(LPI₅) ， 主 要 透 過 全 波 形 資 料 中 之 cross section (ς) 指標來推導出雷射照射面積。從 cross section 的定義(Wagner et al., 2006) (式 5)，我們得 知雷射光接觸物體的面積可由物體反射率和 cross section 計算而得(式 6)。

ς =^4𝜋_Ω𝜌𝐴 (5) 式中，ς =cross section。

Ω =物體反射角≈ 2π。

𝐴 =雷射光與物體接觸的面積。

𝜌 = 物體反射率。

A =_4𝜋𝜌^Ω ς =_2𝜌^σ ，Ω = 2π (6)

而穿透率的指標 LPI5即定義為(雷射光到地面 的面積總和)/(雷射光到地面的面積總和+雷射光打 到樹冠的面積總和)，如式 7，式中 k 值乃地面反射 率和植被反射率之比值，由於單一回波的雷射光打 到地物的面積可視為相同大小，因此由(式 8)之定 義得知，利用單一回波的反射強度值可估算雷射光 在森林中地面和非地面的反射率比例(亦即求 k 值)，

再帶入(式 7)即可求出 LPI5之數值。

LPIS= Σ𝐴𝑔 Σ(𝐴𝑔+ 𝐴𝐶)=

Σ𝜎𝑔

𝜌𝑔

Σ (𝜎_𝑔 𝜌𝑔+𝜎𝐶

𝜌𝐶)=

_Σ(𝜍^Σ𝜍𝑔

𝑔+𝑘𝜍_𝐶)，𝑘＝^𝜌_𝜌^𝑔

𝑐 (7) ρ = 𝐼/𝐴 ， 𝐼 =反射強度值 (8)

Waveform LiDAR Points

DEM

LPI Calculation

Regression LAI

Calculation Chose area

System LiDAR Points

Data processing

In-situ data

基本上這五種穿透率指標都是利用穿透過樹 冠的地面點雲，來計算樹冠的穿透性。在這些指標 中不同的地方在於有的使用點數上的比例，有的是 計算強度的比值，或是利用隱含全波形資料中的雷 射光照射地物的面積比去計算。當光達點雲資料均 勻且完整的分布在森林中時，LPI₁可以以簡單的點 數比就能表示森林樹冠的疏密程度。但是在茂密的 樹冠區中，點雲於冠層的個數會受到雷射光的波長 影響萃取時的解析能力，造成較茂密區的冠層點未 必比較稀疏區的冠層點多，造成點數計算上的錯誤。

因此出現 LPI₂和 LPI₃，分別以強度值或直接捨棄 樹冠層點數用雷射數進行比例上的計算。而 LPI₄ 也是以類似的概念，增加直接穿透點雲的比例，來 強調在較稀疏區域的穿透性。最後一個穿透率指標 則是利用了全波形資料中的資訊，計算出雷射光照 射到地面的面積比，期望增加隱含在波形中的特徵 提高對推估葉面積指數的準確性。

2.5.2 尋找最適合計算範圍

本研究使用 LAI-2000 的地面觀測資料為地真 資料，與計算得的光達穿透率指標進行迴歸分析。

由於 LAI-2000 的設計為接收鏡頭上方距離天頂角 0∘~74∘內的光線，而其實際接收範圍會受到整體 森林的種類和型態、樹冠的密度等影響，無法確定 儀器接受的光照範圍(彭炳勳，2007)。對應光達點 雲資料，各範圍內計算出的穿透率指標因點雲的分 佈或是掃描時的密度，都會有所不同。因此為了尋 找本區域最適合計算穿透率指標來推估 LAI 的範 圍大小，本研究測試了距離樣本點中心半徑 1 公尺 至 15 公尺圓範圍內的光達點雲分別計算 LPI(圖 5)，

並透過與現地 LAI 資料迴歸找出最佳的迴歸參 數。

2.5.3 葉面積指數的迴歸模式和評 估方式

由於森林範圍廣大，在進行森林研究時難以實 際量測到所有的資訊，因此常利用迴歸的方式，以 測量的部分的樣本資料來進行整體的推估。本研究

利用現地量測的 LAI 為依變數(y)，LPI 為自變數(x)，

進行線性模式的迴歸(式 9)。利用已知的樣本點，

計算方程式中的參數(β0、β1)。

, (9) 利用迴歸的方式可以了解穿透率指標和葉面 積指數間是否相關、相關的方向與強度，並且可以 透過判定係數來估計迴歸方程式之適合度。而相關 性判定係數（Coefficient of determination, R²）即是 迴歸造成的平方和(SSR)占總平方和(SST)的比例，

利用最小二乘法的概念，找出誤差項平方和、迴歸 項平方和計算得來。(式 10~12)

SST = ∑^𝑛_𝑖=1(𝑦_𝑖− 𝑦)² (10) SSR = ∑^𝑛_𝑖=1(𝑦̂ − 𝑦̅)𝑖 ² (11) 𝑅²=^𝑆𝑆𝑅_𝑆𝑆𝑇=^∑_∑^{𝑛𝑖=1(𝑦}_(𝑦^{̂−𝑦̅)𝑖 2}

𝑖−𝑦̅)²

𝑛𝑖=1 (12) R² 的數值範圍介於 0~1。當 R²愈大，表示估 計的誤差量越小，代表此迴歸模式能夠解釋全體 y_i變異量的比例愈大。因此 R²愈接近 1.0，代表此 模式愈有解釋能力。但當樣本數(n)較少時，對於 式 12 而言，容易會有高估的現象，因此本研究中 利用校正後 R²（式 13）進行評估，減少樣本數量 造成的膨脹效果。最後使用迴歸得到之參數建立模 型，來推估大範圍葉面積指數得分佈。

𝑅𝑎𝑑𝑗2 = 1 −^(1−𝑅_{𝑛−𝑘−1}^2)(𝑛−1)，𝑘 = 解釋變異之個數 (13) 在計算出最佳的迴歸方程式後，即可透過已知 的參數建立推估方程式來推算整區森林的葉面積 指數。

3. 研究成果分析

3.1 樹冠高度模型估計成果

以三種光達系統之多重回波及全波形資料，A 及 B 實驗區皆可獲得六種點雲(三種系統點雲、三 種全波形點雲)資料來分別產製 DEM、DSM 及 CHM。而不同點雲資料所產製的 CHM 都有些許差 異，本研究利用六種 CHM 在相同網格位置的高度

i

i x

y ₀₁ ⁱ^1,...,n

值進行平均值及標準偏差的計算，得到的成果如圖 6，標準偏差之直方圖如圖 7。從圖 6 我們可以明 顯的看出，在森林區的邊緣、植生厚度差異較大的 地方，由不同點雲資料得到的 CHM 高度差異也較 大，導致標準偏差容易會有較大的偏差值產生。

CHM 的資訊除了可以用來推估森林中林木之 樹高(圖 8)，也可以用來計算森林的樹體積。Lefsky et al.,(1999)提出樹體積的定義為包含樹本身的體

積以及樹冠底下的空間，相對於樹冠高度模型來說，

樹體積可利用樹冠表面到地表的距離累積量來得 到。應用本研究中所製作之 CHM 來推估兩個實驗 區整體範圍的樹林體積如表 4，整體的標準偏差占 整體體積的 2%~3%。顯示了用不同光達點雲製作 出的 CHM，在森林樹體積得推算上沒有太大的差 異。

(a) (b)

圖 5 光達穿透率指標計算範圍(a)側視示意圖(b)俯視示意圖

(a) (b)

(c) (d)

圖 6 A、B 區分別之 CHM 網格點的平均值(a)(b)及標準偏差(c)(d)分佈圖 Radius of area

(1~15)m

Radius of area

(a) (b)

圖 7 六種點雲資料產製之 CHM 標準偏差統計直方圖(a) A 區；(b) B 區

(a) (b)

圖 8 CHM 高度統計直方圖(a)A 區；(b)B 區(橫軸為高度(m)；縱軸為累計量) 表 4 CHM 用於樹體積之推估

A 區 (306m*246m)

樹體積(m³)

B 區 (256m*269m)

樹體積(m³) Leica 多重回波 370961 630319

Leica 全波形 366338 676975 Reigl 多重回波 381245 677702 Reigl 全波形 379134 674281 Optech 多重回波 391125 685498 Optech 全波形 373516 665603 平均樹體積 34788 19895 標準偏差(std) 8767.653 19722.94 std 占總體積的比例 2.3% 3.0%

平均植生厚度(m) 4.928 9.1531 (平均植生厚度：dsm 的平均高度，可顯示出在該 區域內平均植生的表面高度)

3.2 森林葉面積指數估計成果

本研究中利用 Leica 和 Riegl 兩種光達系統的 資料來進行葉面積指數的推估，除了比較此兩種系 統外，在 Leica 光達系統上，我們亦有獲取不同航 高的點雲資料，以茲進行不同航高是否會對穿透率 指標之計算有所影響。在進行穿透率指標計算之前，

必須將所有點雲進行分類處理，把點雲分成地面和 非地面點，且利用人工的方式將非地面點中如房屋 等高於一公尺的人工建物去除，避免建物的低穿透 率造成森林葉面積指數的高估。接著將所有點雲值 與產製出的 DEM 相減，讓每個點雲的高程值都為 地面上高度。得到含地面高度的點雲後，我們設定 高程門檻值為 1 公尺，將高程為 1 公尺以下的點雲 設成本研究穿透率指標中的地面點(圖 9)。

經過資料的前處理後，即可利用各個點雲分類 後的屬性，包含回波數、強度值資訊進行穿透率指 標的計算。其中，由於 LPI₅必須要有 cross section

的資訊，因此只適用於全波形點雲，而其計算公式 中的 k 值透過選取單一回波中地面點與非地面點 的比值得到，在各光達資料中使用之 k 值如表 5。

由其中的比值可以發現，地面點與非地面點的比值 在兩種光達系統的資料下並沒有相同的趨勢。其原 因可能是因為在森林中的地表有許多植生覆蓋，加 上本研究中將 1 公尺以下的點都設成非地面點，在 這 1 公尺高度範圍內可能包含低矮植物、草叢等類 似於樹冠層反射的地物，因此得到的結果並沒有一 致。

表 5 單一回波 intensity 統計平均值 地面點 非地面點 k Leica 光達系統

(低航高) 90.0 96.2 0.94 Leica 光達系統

(高航高) 48.6 51.1 0.95 Riegl 光達系統

(低航高) 60.4 51.5 1.17 為了要和地面觀測資料進行比較，LPI 計算範 圍的選取方式是根據地面量測樣本點的中心位置，

切割半徑 1~15 公尺圓柱範圍內的點雲進行計算。

接著將計算之成果與現地測得之葉面積指數做線 性迴歸的分析，得到在不同範圍內的相關性判定係 數 (R2)。利用 LPI 與 LAI 之間的 R2 來判定 LPI 在不同大小的點雲範圍內對於 LAI 解釋的能力，

並選定較好的計算範圍得到 LPI 推估 LAI 的估計 方程式。

由圖 10 可以看出，該區穿透率指標對葉面積

指數最大的解釋能力約在 13、14 公尺的計算範圍。

表示出當利用 13、14 公尺範圍內的點雲計算的穿 透率指標，可能和現地量測資料接受的光量範圍有 較一致的現象。因此在本研究中選取最佳推估半徑 來進行計算，得到的各穿透率指標與葉面積指數的 R2 如表 6 所示。整體而言，在 Riegl 的資料下，各 種穿透率指標對於 LAI 都有較好的解釋能力；而 利用 Leica 多重回波光達點雲計算得到的成果則是 最差。其可能的原因為 Leica 多重回波資料的光跡 較小，且在一個反射回訊中只能得到最多四個點雲，

容易造成弱回訊的地物被遺漏，造成樹冠層點雲和 地面點之間的分佈無法完整表達森林林分的結構，

使得推估的成果較差。

對於不同的穿透率指標而言，LPI₄在各種的光 達點雲下都能得到高於 0.5 的 R²值，顯示出該指 標對於估計葉面積指數上有較穩定且準確的表現，

其與現地觀測葉面積指數的線性迴歸模式展示如 圖 11。表現次之的為利用強度值所獲取的 LPI₂， 此指標尤其應用在全波形點雲所計算的 R²值有較 明顯的提升。顯示利用對點雲之波形振幅大小所賦 予的強度值能有效的推估 LAI，對於兩種的全波形 資料而言，利用強度值計算的穿透率指標(LPI2)在 估計葉面積指數上都有良好的表現。而本研究中利 用面積概念推導出 LPI5，雖然與其他的穿透率指標 相比之下並非最具解釋力的指標，但在兩個低航高 的全波形資料下都能達到 70%以上的預測能力。

圖 9 前處理後的點雲 (橘色點為原分類的地面點；白色點為非地面點；綠色點為經 1 公尺門檻後增加的 地面點)

(a) (b)

(c) (d)

(e) (f)

圖 10 葉面積指數和(a) Leica 低航高的多重回波點雲；(b) Leica 低航高的全波形點雲；(c) Leica 高航高的 多重回波點雲；(d) Leica 高航高的全波形點雲；(e) Riegl 低航高的多重回波點雲；(f) Riegl 低航高的 全波形點雲資料中的各種穿透率指標在不同範圍下的相關性判定係數(R²)

表 6 穿透率指標與葉面積指數的相關性判定係數 R²

LPI₁ LPI₂ LPI₃ LPI₄ LPI₅ Leica 多重回波光達點雲(低) 0.30253 0.43644 0.20938 0.55610

Leica 全波形光達點雲(低) 0.48964 0.77171 0.46812 0.68276 0.72653 Leica 多重回波光達點雲(高) 0.63376 0.51218 0.53831 0.61938

Leica 全波形光達點雲(高) 0.66135 0.51752 0.55851 0.64417 0.47414 Riegl 多重回波光達點雲(低) 0.80765 0.67251 0.80174 0.85793

Riegl 全波形光達點雲(低) 0.73431 0.83006 0.83663 0.82152 0.78721

(a) (b)

(c) (d)

(e) (f)

圖 11 (a) Leica 低航高的多重回波點雲；(b) Leica 低航高的全波形點雲；(c) Leica 高航高的多重回波點雲；

(d) Leica 高航高的全波形點雲；(e) Riegl 低航高的多重回波點雲；(f) Riegl 低航高的全波形點雲資料 中的 LPI4 與現地觀測葉面積指數的直線迴歸模式

除了比較不同光達系統的資料外，本研究也加 入了 Leica 光達系統高航高的資料，利用已知的 DEM 進行穿透率指標的計算來推估葉面積指數。

相較於單獨使用低航高的多重回波點雲資料，利用 點數比計算得到的 LPI₁、LPI₃及 LPI₄，都有明顯 的提升。但對於 LPI₂和 LPI₅而言，由於不同航高 的強度值及雷射光光束照射地物的面積的大小不 同，因此若是低航高和高航高的資料一起做計算，

容易會產生估計上的錯誤，如圖 12。由此結果能 看出，雖然高航高的的掃描方式會因為距離較遠造 成雷射光減弱，而較難得到地面點的回訊。但對於

利用點數比例來計算的穿透率指標而言，對於估計 葉面積指數還是能有不錯的成果。

利用 LPI 與地面觀測之葉面積指數樣本資料 進行迴歸，即可得到推估葉面積指數的估計方程式，

如表 7。本文以低航高的資料中，選取表現較佳的 LPI₄來進行測區附近約 350 公尺×350 公尺範圍的 葉面積指數估計。計算的方式為將測區範圍進行網 格的切割，計算每單位網格內的點雲穿透率指標，

接著利用估計方程式推算各網格內的葉面積指數，

即可得到整區的成果(圖 14)。

(a) (b)

圖 12 葉面積指數和各種穿透率指標在不同範圍下的相關性判定係數(R²)(a) Leica 低航高+高航高的多重 回波點雲；(b) Leica 低航高+高航高的全波形點雲

表 7 各 LPI 推估 LAI 之估計方程式

LPI₁ LPI₂ LPI₃ LPI₄ LPI₅

Leica 一般光達光 達點雲(低航高)

1=-20.0789

=6.6153

1=-29.1139

=5.7236

1=-14.1775

=6.1526

1=-19.0732

=6.9327 Leica 全波形光達

點雲(低航高)

1=-14.0471

=6.8429

1=-26.4962

=6.8142

1=-10.8338

=7.1862

1=-13.4004

=6.9551

1=-20.067

=7.1522 Leica 一般光達光

達點雲(高航高)

1=-14.9156

=5.6945

1=-30.1255

=5.5857

1=-12.1656

=5.5517

1=-14.0729

=5.9115 Leica 全波形光達

點雲(高航高)

1=-10.8053

=5.5327

1=-13.3573

=5.3523

1=-7.2155

=5.2504

1=-8.8458

=5.5073

1=-26.2298

=5.2351 Riegl 一般光達光

達點雲(低航高)

1=-15.7229

=7.2919

1=-13.8335

=6.5266

1=-11.5581

=7.5642

1=-14.95

=7.4346 Riegl 全波形光達

點雲(低航高)

1=-14.5435

=6.7992

1=-14.7289

=6.484

1=-11.6239

=7.407

1=-13.2666

=6.7473

1=-15.2182

=6.3788

圖 13 以不同點雲分別計算所得到的推估成果，

設定的網格大小為 1 平方公尺。由於本研究估計的 是森林的葉面積指數，當得到的葉面積為零時即表 示該區在單位面積上的林木沒有葉子，與正射影像 相較下，可以發現穿透率指標利用點雲在垂直方向 的分佈，能大致區分出樹林和非樹林的區域。

就整個區域而言，以不同光達點雲推估的成果 顯示，利用估計方程式得到的整區葉面積指數彼此 之間的 R²可達到 0.8 以上，標準偏差在 0.9 以下(表 8)。從四種點雲計算得到的標準偏差分佈圖如圖 14 中，可以發現在不同光達點雲下使用穿透率指標進 行大範圍面積的葉面積指數估計時，容易在森林區 與非森林區邊界、草地等，與現地取樣區域樣貌差

異較大的地方產生較大的錯誤。其中最大的葉面積 指數差值可達到 3 以上，由不同光達系統得到的剖 面圖如 15 可以看出，應該是由於時間差造成地物 變化所造成。顯示了利用穿透率指標推估葉面積指 數的方法，應用在不同光達系統點雲下的成果，除 了在森林邊界、樹冠高度差異較大的地方外，都能 夠呈現一致的成果。

表 8 利用 LPI₄得到整區葉面積指數估計之比較 標準偏差 Leica 多重回波 0.8414

Leica 全波形 0.6445 Riegl 多重回波 0.7263 Riegl 全波形 0.6194

i

i x

y ₀ ₁

(a) (b) (c)

(d) (e)

圖 13 利用整區葉面積指數的計算成果圖(a)測區正射影像；(b) Leica 多重回波;(c) Leica 全波形;(d) Riegl 多重回波;(e) Riegl 全波形

圖 14 四種不同點雲資料得到的葉面積指數估計之標準偏差分佈圖，圈選處為差異較大的地方

(a) (b)

圖 15 (a)Leica 和(b)Riegl 光達系統在圖 14 圈選處的點雲的剖面圖

4. 結論

空載光達技術能直接獲取森林結構的三維空 間分布，有助於樹冠高度模型的製作。但其成果的 精度受限於點雲的密度，若得到的地面點過於稀疏，

則難以表示實際的地形起伏情況。因此實際應用上

在做人工編修時，應對於地面點點數較少的區域進 行加強檢視，以防止有較大的誤差產生。本研究中，

三種光達系統的點雲密度在每平方公尺內都有 2 個點以上，在森林範圍內所製作出來的數值高程模 型，其相對精度皆在 50 公分以內。而在 DSM 製 作上，容易在樹冠孔隙、及森林邊界的地方有較明

顯的差異，尤其是在林木較高的森林區中，最大的 高度差異值超過 10 公尺。利用 1 公尺網格的 DSM 及 DEM 產製出的樹冠高度模型的標準偏差約在 1.5 公尺內，顯示利用不同光達系統得到的數值模 型對於森林的應用有足夠的一致性，能提供有效且 快速的量測方法。

對於葉面積指數估計方面，實驗結果發現當使 用高點雲密度、地面點點數比例較高的光達資料，

在推估森林葉面積指數能有較好表現。其中，以單 一回波地面點加權計算的 LPI₄在不同光達點雲資 料下的表現最為穩定，各種資料計算出的成果對於 估計葉面積指數都能達到 50%以上的預估能力。就 單一航高的資料中，利用全波形回訊萃取出的點雲，

在計算以強度值比例的 LPI₂時，相較於多重回波 而言有明顯的提升。顯示經波形萃取出的振幅大小，

在樹冠層及樹冠下地表的比例能有效的反映森林 結構，進而有利於推估葉面積指數。對於不同航高 的 Leica 多重回波光達資料而言，利用高航高的點 雲計算以點數為比例的穿透率指標(LPI₁、LPI₃、 LPI₄)，再進行葉面積指數的估計上有較好的表現。

但若是將不同航高的資料疊加進行穿透率指標的 計算，其中的 LPI2和 LPI5會因為不同航高在掃描 時，光束打到地物的雷射強度及光跡大小有所不同，

造成得到的強度值及面積的比例無法表示實際的 現象，使得與現地資料進行迴歸的相關性有下降的 情形。

由於葉面積為一個動態的森林參數指標，會隨 樹種、樹高的變化而改變，因此應在量測前進行森 林樹種的調查，並按樹林樣貌進行劃分來分別進行 葉面積指數的估計。而本研究中所使用的現地觀測 資料與光達掃描時間相差超過 1 年，且量測季節略 有不同，森林生態環境的改變對於葉面積指數的估 計可能會造成影響，若是能夠使用時間較相近的資 料進行推估，較能獲得準確的估計方程式，來達到 以小範圍的現地實測資料即可推估大範圍及無法 到達之森林區域的資料。目前空載光達能進行飛航 掃描的時間受限許多環境及政策等因素影響，或許 能發展使用衛載光達進行森林葉面積指數的探測，

以達到長時間對森林的進行監控的目的。

致謝

感謝國立屏東科技大學森林系借用 LAI-2000

儀 器 ， 以 及 國 科 會 計 畫

(NSC101-2221-E-006-181-MY3)經費的贊助，讓本 研究能順利進行，特此致謝。

參考文獻

吳守從，陳永寬，2004。應用數位航測技術探討南 仁山生態保護區林分高生長量，航測及遙測學 刊 第九卷 第四期 第 21-34 頁

黃清美，2007。空載光達點雲穿透率探討，國立交 通大學土木工程學系碩士學位論文。

彭炳勳，2007。應用空載光達資料推測林木樹高與 葉面積指數，國立屏東科技大學森林系碩士學 位論文。

蕭淳伊，2008。應用空載光達資料與遙測影像推估 樹林分布及體積，國立成功大學測量及空間資 訊學系碩士學位論文。

Arnaud, Jacob, Yi-Hsing Tseng, Chi-Kuei Wang (2011), “Investigating the penetration capability of airborne LiDAR for high vegetation of subtropical evergreen forest.” The 30th conferencn on Surveying and Geomatics, Taiwan.

Danson, FM, Morsdorf, F and Koetz, B (2009),

“Airborne and terrestrial laser scanning for measuring vegetation canopy structure”, in:

Laser Scanning for the Environmental Sciences , Wiley-Blackwell.

Hopkinson Chris, Laura Chasmer (2009), “Testing LiDAR models of fractional cover across multiple forest ecozones” Remote Sensing of Environment 113, Pages 275–288

Jonckheere Inge, Stefan Fleck, Kris Nackaerts, Bart Muys, Pol Coppin, Marie Weiss, Frédéric Baret (2004)“Review of methods for in situ leaf area index determination Part I. Theories, sensors and hemispherical photography” Agricultural and Forest Meteorology 121, 19–35

Koetz B., F. Morsdorf, G. Sun, K. J. Ranson,K. Itten, B. Allgöwer (2006) ” Inversion of a Lidar Waveform Model for Forest Biophysical Parameter Estimation IEEE Geoscience And Remote Sensing Letters, Vol. 3, No. 1 P.49-P53 Lefsky, M. A., Cohen, W. B., Acker, S.A., Parker, G.

G., Spies, T. A. , Harding, D.(1999), “Lidar remote sensing of the canopy structure and biophysical properties of Douglas-fir western hemlock forests” Remote Sensing of Environment, Vol. 70(3), pp. 339-361

Masona Euan G.,, Mathijs Diepstratenb, Guy L.

Pinjuvc, Jean-Pierre Lasserred (2012)

“Comparison of direct and indirect leaf area index measurements of Pinus radiata D. Don”

Agricultural and Forest Meteorology 166– 167 113– 119

Morsdorf F., O. Frey, E. Meier, K.I. Itten , B.

Allgöwer (2008) “Assessment of the influence of flying altitude and scan angle on biophysical vegetation products derived from airborne laser scanning” International Journal of Remote SensingVolume 29, Issue 5, pages 1387-1406 Næsset E.(1997), “Estimating timber volume of

forest stands using airborne laser scanner data“ Remote Sensing of Environment Volume 61(2), Pages 246–253

Solberg, S., Næ sset, E., Hanssen, K.H. and Christiansen, E. (2006), “Mapping defoliation during a severe insect attack on Scots pine using airborne laser scanning.” Remote Sensing of Environment, 102, pp. 364–376.

Solberg, S. (2008),”Comparing discrete echoes counts and intensity sums from ALS for estimatingforest LAI and gap fraction.” In Proceedings of SilviLaser: 8th International Conference on LiDAR Applications in Forest Assessment and Inventory, R.A. Hill,J. Rosette.

and J. Sua´rez (Eds), September 2008, Edinburgh, UK, pp. 446–455.

Solberg, S. (2010),”Mapping gap fraction, LAI and defoliation using various ALS penetration variables.” International Journal of Remote Sensing, 31(5): 1227-1244.

Wagner, W., A. Ullrich, V. Ducic, T. Melzer and N.

Studnicka, (2006)” Gaussian decomposition and calibration ofa novel small-footprint full-waveform digitising airborne laser scanner”

ISPRS Journal of Photogrammetry and Remote Sensing 60(2), pp. 100-112.

Zhao Kaiguang, Sorin Popescu (2009) “Lidar-based mapping of leaf area index and its use for validating GLOBCARBON satellite LAI product in a temperate forest of the southern USA” Remote Sensing of Environment 113:

1628–1645

1 Master, Department of Geomatics, National Cheng Kung University Received Date: Mar. 21, 2014

2 Ph.D, Department of Geomatics, National Cheng Kung University Revised Date: Oct. 13, 2014

3 Professor, Department of Geomatics, National Cheng Kung University Accepted Date: Oct. 30, 2014

4 Assistant Professor, Department of Geomatics, National Cheng Kung University

*.Corresponding Author, Phone: 886-6-2383876 ext. 852, E-mail:s3350839@gmail.com

Estimation of Forest Canopy Height Model and Leaf Area Index Using Airborne LiDAR data

Li-Ping Lin¹ Cheng-Kai Wang^2* Yi-Hsing Tseng³ Hone-Jay Chu⁴

ABSTRACT

Efficiently obtaining the information in forest region such as forest structure, forest ecosystems is important for forest management. The purpose of this study is to estimate the Canopy Height Model (CHM) and the Leave Area Index (LAI) of a dense forest area by using airborne LiDAR data. CHM is estimated by taking the difference of DSM and DEM derived from LiDAR data. Estimation of LAI is achieved based on the calculation of Laser Penetration Index (LPI). Five calculations of LPI were applied in this paper: (1.) The ratio between the number of ground points and that of all the points; (2.) the ratio between the intensities of ground points and that of all the points; (3) the ratio between the number of ground points and the number of laser beams;

(4) a weighting method modified from index (1); and (5) the ratio between the area of ground points and that of all the points.

The study area is in a natural broadleaf forest of south Taiwan. In this study, we use three sets of airborne LiDAR data acquired with different full waveform LiDAR systems including Leica ALS60, Riegl LMS-Q680i and Optech Pegasus HD400. All of these LiDAR systems are capable of recording full waveform data, then we can get the waveform point clouds by the echo detector to do the comparison. Our experiments results show that the accuracy of CHM by different LiDAR data is about 1.5 meter. And the fourth LPI index has the highest coefficient of determination (about 0.8) and the estimation of LAI can be improved by using the waveform points.

Keywords: