應用NDVI植生指標與平均值 調整影像分割法於崩裸地萃取
—以六龜試驗林地區多期福衛二號影像為例
謝漢欽1,3) 鍾智昕1) 黃俊元2)
摘 要
面對高解析度光學衛星影像用於地表目標物的分類，傳統上以像元為基礎的分類方法，往往效果 不令人滿意；近年來採用影像分割技術，以物件為基礎的影像分類程序可以得到更好的結果。本研究 以高雄市六龜試驗林及其周邊29,400公頃土地為試驗區，選用福衛二號2006年、2009年及2011年3期夏 季影像，經影像輻射校正前期處理及不改變光譜特性的全色影像銳化融合處理，建立各期影像NDVI指 標，可減除地形效應。應用具多項優點的平均值調整影像分割法及簡易物件形狀指標及坡度法則篩選程 序，自動萃取3期試驗區的崩裸地。使用樣區系統取樣法，以各期樣區人工數化的崩裸地當作地真值，
改進崩裸地自動萃取的精度評估。結果3期綜合分類準確度皆達95%以上，Kappa係數達85%以上，顯示 本研究的影像分割程序對多期FS2影像的崩裸地的自動萃取具有好的結果，具可應用性。從3期崩裸地 面積變遷分析結果可得知，2009年莫拉克颱風造成非崩裸地轉變成崩裸地的轉移率是後續兩年期間的 3.5倍，而2009年至2011年間由崩裸地轉為非崩裸地的植生復育率是2006年至2009年的18.7倍。此外將 2009年至2011年新增的崩裸地，透過數值地形衍生的向源河系套疊，可證明新增崩塌地多位於原有崩塌 地往上游之向源侵蝕區域。本研究成果有利於六龜試驗林莫拉克颱風後崩塌地的空間治理。
關鍵詞：影像分割、崩塌地、福衛二號、植生指標。
謝漢欽、鍾智昕、黃俊元。2017。應用NDVI植生指標與平均值調整影像分割法於崩裸地萃取―以六龜 試驗林地區多期福衛二號影像為例。台灣林業科學32(3):203-22。
1)林業試驗所森林經營組，10066台北市南海路53號 Division of Forest Management, Taiwan Forestry Research Institute, 53 Nanhai Rd., Taipei 10066, Taiwan.
2 )林 業 試 驗 所 太 麻 里 研 究 中 心 ，9 6 3 4 1 台 東 縣 太 麻 里 鄉 大 王 村 橋 頭 6 號 Ta i m a l i R e s e a r c h C e n t e r, Taiwan Forestry Research Institute, 6 Bridge Header, Taimali Township, Taitung County 96341, Taiwan.
3)通訊作者 Corresponding author, e-mail:[email protected]
2016年11月送審 2017年1月通過 Received November 2016, Accepted January 2017.
【Summary】
When applying high-spatial-resolution satellite imagery to classify of earth surface targets, traditional pixel-based classification methods often produce unsatisfactory results. In contrast, using object-based image classification with image segmentation approaches over the last decade achieved further improvements. In this study, 29,400 ha of forestland involving the Lioukuei Experimental Forest and the surrounding area in Kaohsiung City was used as the experimental area. Three Formosat-II (FS2) images of the experimental area obtained in the summers of 2006, 2009, and 2011 were respectively selected and preprocessed through radiometric calibration and pan-sharpening fusion without distorting their spectral characteristics. The normalized difference vegetation index (NDVI) of each FS2 image was established for follow-up image segmentation with a mean-shift algorithm. After NDVI segmentation with the mean-shift method, landslide and bare (LSAB) areas for each date were automatically extracted using a simple rule including a shape index and average slope of an object to filter out noise objects and mountain village areas. To improve the accuracy assessment for extracting of LSAB areas for each date, an assessment was performed by comparing to manually digitized sub-images of each pan-sharpened image through a system sampling method. Results showed that the overall accuracies of the 3 dates were > 95%, and their Kappa coefficients were > 85%; thus, the mean shift procedure can successfully be applied to extract LSAB areas from multi-date FS2 images. The area transition rate of Un-LSAB to LSAB areas between 2006 and 2009 after typhoon Morakot was 3.5-times that in the following 2 yr, and from 2009 to 2011, the area transition rate of vegetation regeneration was 18.7-times that between 2006 and 2009. Inspecting the spatial distribution of newly increased LSAB areas from 2009 to 2011 with extended tributaries derived from digital elevation model, transition areas from vegetation cover to LSAB were mostly located in headward areas.
Key words: NDVI, mean shift segmentation, FORMOSAT-2, landslide.
Hsieh HC, Chung CH, Huang CY. 2017. Using the NDVI and mean Shift segmentation to extract landslide areas in the Lioukuei Experimental Forest region with multi-temporal FORMOSAT-2 Images. Taiwan J For Sci 32(3):203-22.
緒 言
從1999年IKONOS衛星影像問世以來，
光 學 資 源 衛 星 影 像 邁 入 高 解 析 度 的 新 紀 元 。 台 灣 也 於2 0 0 4 年 5 月 起 有 了 自 主 的 福 衛 二 號 (FORMOSAT-2，以下稱簡稱FS2)資源探測衛 星 ， 使 得 本 土 的 高 解 析 度 衛 星 影 像 在 時 間 及 空間解析度上皆顯著提升(Liu 2006, Lin et al.
2013)。FS2影像空間解析度在多光譜方面有 8 m，全色光有2 m，兩者透過適用的光譜融 合技術(Zhang 2004, Liu et al. 2011)，可使得 多光譜資料具有2 m高空間解析度。然而隨著 空間解析度的提高，應用高解析度光譜影像進 行林地目標物分類時，如果使用以單一像元為 基礎(pixel-based)的分類方法，通常僅能利用 影 像 的 光 譜 和 紋 理 資 訊 進 行 分 類 ， 忽 略 了 影 像上下文(context)與形狀等空間特徵(Zhou et al. 2009, Blaschke 2010)。對於高空間解析度 影像，如果僅依靠像元的光譜資訊對影像進行 分類，往往會造成嚴重的「胡椒鹽」現象，導 致分類圖的區塊破碎度高，難與現有地理資訊 圖層進行整合；一般還需要透過人機交談的方 式，進行影像後處理及分類圖編修等工作，才 能滿足應用的需求(Blaschke et al. 2014, Kun et al. 2014)。
近 年 來 運 用 各 種 影 像 分 割 技 術(Schiewe 2002, Neubert et al. 2008, Dey et al. 2010)，將 影 像 分 割 後 的 連 續 型 區 塊 當 作 物 件(object)，
運 用 於 以 物 件 為 基 礎 件(object-based)的分類 (Kamagata et al. 2008)。分類時不僅依靠地 物 的 光 譜 特 徵 ， 更 可 加 入 物 件 的 形 狀 、 紋 理 (texture)等幾何特徵和空間關係等資訊，把具 有相似特徵的物件進行分類，可以克服以像元 為基礎的影像分類的缺點(Blaschke et al. 2011, Xiao and Feng 2012, Zhao et al. 2014)。國內外 相關研究，將其應用在高解析度影像的土地覆 蓋類型或土地利用型分類上，往往獲得良好的 結果(Chang and Chen 2006, Wulder et al. 2008, Lei et al. 2009, Boukir et al. 2012, Kun et al.
2014)。此外長期以來由衛星影像的近紅外及 紅光光譜波段的衍生的Normalized Difference Vegetation Index (NDVI)指標，已能有效率
地 分 辨 出 植 生 覆 蓋 與 非 植 生 覆 蓋 的 地 物 或 區 域(Lyon et al. 1998, Hayes and Sader 2001, Nathalie 2013)，可用於崩裸地的偵測與萃取 (Huang et al. 2010, Liu 2011, Hussain et al.
2013, Lin et al. 2013)。
本研究針對2006年至2011年間飽受颱風、
暴雨天然災害的國有林地熱區－高雄市林試所 六龜試驗林轄區及周邊土地，以之當作試驗區 域。選用2009年莫拉克颱風後及其前後各1期 共3期，拍攝日期適當(接近夏季，太陽高度角 較高)地形陰影較少且無雲覆蓋的FS2影像。建 立各期影像的NDVI植生指標，採用效果深受好 評的平均值調整影像分割法(mean shift image segmentation) (Chen 1995, Comaniciu and Meer 2002, Lin et al. 2006, Bo et al. 2009, Boukir et al. 2012, Mohan and Leela 2013, Kun et al.
2014)，從各期NDVI影像自動萃取試驗區的崩 裸地，評估各期崩裸地萃取的準確度。此外，
在應用面上，進行3期影像自動萃取崩裸地變遷 分析，期望從分析結果得知六龜試驗林轄地的 在這2個時段的崩裸地消長及其在空間分布的差 異性，以及證明2009年莫拉克颱風後至2011年 崩裸地增加的區位，多位在原有崩塌地往上游 發展的向源侵蝕的區域，據以提供六龜試驗林 經營管理有關崩塌地治理的空間決策資訊。
材 料 與 方 法
一、研究試驗地與材料
本研究以高雄市林業試驗所六龜試驗林轄 區及其周邊土地共29,427 ha為試驗區(研究試 驗區位置圖如Fig. 1)，採用本區2006年9月28 日、2009年8月19日(2009年8月9日莫拉克颱風 過後)及2011年8月17日拍攝，幾乎無雲，太陽 高度角較高(角度依次為54.5°、58.5°及60.1°) 的3期FS2衛星影像。在地理資訊方面，使用林 業試驗所六龜試驗林地理資訊圖層、高雄市電 子 地 圖(提供山村聚落、道路及溪流等地理空 間資訊)，以及試驗區域5 m的數值高程(digital elevation model)資料(用於影像分割後，區塊平 均坡度計算)。
二、研究方法
(一) 影像前處理與植生指標的建立
依據3期影像所附的詮釋資料(meta data) 將灰度轉換為輻射值，並以近紅外光波段為基 準 ， 選 取1％暗值像元，依各波段波長特性去 除路徑輻射(path radiance)。再將各波段輻射 值除以影像感測當時的太陽入射角餘弦，以減 除 不 同 拍 攝 日 期 太 陽 入 射 角 差 異 的 入 射 效 應 (Schroeder et al. 2006)。以2009年為基準影像，
分別將2006年、2011年對2009年，選取平均分 布30個地面控制點，進行影像間空間幾何再對 位，使得不同時期影像間套疊誤差減小。採用 PCI Geomatics影像處理系統的全色銳化(pan sharpening)影像融合方法(Zhang 2004)，本方 法有別於一般傳統的影像融合的方法(如HIS、
PCA、Brovey、SVR、RE及wavelet等)，採用 統計最小平方法計算方法，以求取個別波段間
灰度值最佳適合值，據以調整個別波段融合時 灰度值的貢獻量，除了可防止顏色的失調外，
還可減低使用不同衛星感測器時，因感測器光 譜域不同所造成融合影像的差異(Zhang 2002, Zhang et al. 2004)。本研究使用此一影像融合 法，分別將3期FS2多光譜影像與全色波譜融合 處理，不僅可維持融合後多光譜波段與原先的 光譜特性，並能具有2 m的高空間解析度。最後 再以融合後影像之近紅外光波段及紅光波段，
分別建立3期FS2的NDVI(公式如第1式)影像，
以利於後續以NDVI進行影像分割。
NDVI = (NIR－R)
(NIR+R) (1)
其中：NIR為近紅外光波段、R為紅外光波 段，NDVI值域為[-1.0, 1.0]
(二) 平均值調整影像分割的演算法
影像分割(image segmentation)是依據影 Fig. 1. Location of the experimental area in Taiwan. Right image shows a Formosat-II image acquired on 19 August, 2009 overlaid with the boundary (yellow lines) of the Lioukuei Experimental Forest.
像像元的灰度、光譜、紋理、空間相關等相似 特徵，把影像分割成若干互不交疊，具有相對 均質特徵的物件的過程。影像分割是從影像處 理到影像分析的關鍵步驟，是達成地物表達及 特徵測量的途徑；同時也是將原始影像轉化成 以 物 件 為 導 向 的 基 礎 ， 可 利 於 促 成 更 高 階 影 像分析與理解(Blaschke et al. 2011, Zhao et al.
2014)。若從遙測影像分割的觀點，影像分割 的方法很多，一般可依據分割所使用的影像特 徵 ， 將 影 像 分 割 法 分 成6大類1.直方圖門檻值 法、2.特徵空間聚集法、3.以區域為基礎的分割 法、4.以邊緣為基礎的分割法、5.以數學形態學 為輔助的綜合分割法，以及6.多解析度影像分 割法(Zhou et al. 2009, Zhao et al. 2014)。各種 影像分割法都有其優缺點，很難說存在一個最 優影像分割法，如何應用端視試驗區域、研究 對象及影像特徵而異(Xiao and Feng 2012)。
本研究採用的平均值調整(mean shift)影 像分割法，是屬於上述第2.類特徵空間聚集法 的 一 種 有 效 方 法 。 本 方 法 最 先 由Fukunaga及 Hostetler (1975)提出，之後由Cheng (1995)應用 於影像分析，直到2002年由Comaniciu與Meer (2002)及Comaniciu與Ramesh (2003)將之推廣 應用於電腦視覺領域；可應用項目包括影像分 割、自適應平滑處理(adaptive smoothing)及目 標追綜等低階視覺(low-level vision)問題處理上 (Christoudias et al. 2002)。
基 本 上 ， 平 均 值 調 整 影 像 分 割 法 是 一 種 從 多 維 的 特 徵 空 間(feature space)自動聚 集(clustering)的演算法，也是屬於一種無母 數(non-parametric)經驗式的機率密度函估計 (probability density function estimation)方法。
其相較於全局式(global)的影像聚集分類演算法 (如K-means)，使用密度估計的方法具有對先驗 知識需求少，不需先預定期望的聚集(cluster or mode)數目，以及須假定聚集組在特徵空間上 需要符合高斯常態分布的固定形狀(會產生聚集 組間在特徵空間上相互重疊)的優點。可完全依 靠輸入的資料進行估計，適用於局部(local)區 域的離散型資料分布，能進行任意形狀密度函 數的估計，對於不同局部結構的資料，在方法
上有很好的收斂性與穩健性(Christoudias et al.
2002, Lin et al. 2006, Zhao et al. 2014)。
在演算上能使特徵空間的每一個點，通過 有效的統計迭代(recursive)，依據密度函數梯 度上升程序「調整」至密度函數分布的局部極 大值的點(平均值迭代移動的示意說明如Fig. 2) (Derpanis 2005)，此靜止點代表一個聚集組的分 布，與此點相關的點當作此聚集組成員。對於 遙測影像而言，經由設定平均值調整演算的不 同參數值，很容易透過平均值濾波(mean filter) 的移動視窗(moving window)的計算基礎，達 成不同精細程度的影像分割(Kun et al. 2014)。
其 數 學 演 算 公 式 及 電 腦 演 算 程 序 簡 述 如 下 (Comaniciu and Meer 2002, Zhao et al. 2014)：
給 定d維歐幾理得空間Rd中 的n個任意點 {xii = 1, ....n}，在x點的機率密度函數f (x)的核密 度函數K (x) (核函數一般採用均向對稱的函數如 Epanechnikov及高斯常態函數，以之進行核密度 估計就能達成滿意的結果)，其估計公式為：
(2)
式中，k (ǁxǁ2)：核函數的剖面 Ck,,d：常係數
h：移動視窗的半徑
為 找 到 機 率 密 度 最 大 值 的 點 ， 令g ( x ) = -k'(x)對第(2)式求導數，令導數值為0，得到第 (3)式：
(3) 第(3)式第二個中括號內的部分稱之為平均 值調整向量(mean shift vector)：
(4) 平均值調整向量總是指向具有局部機率密
度最大增量的方向，當經由多次的迭代計算，
當 平 均 質 調 整 向 量 為0時，機率密度達到最大 值。在實際操作時，使用第(5)式的迭代公式：
(5) 式中，yj的初始值為x，yj +1是第j 次迭代求 得的特徵值。可以證明序列yj能夠收斂至y'，此 時y'點為 x對應的聚集組。
一幅多光譜影像常被表示為具有p維向量的 二維網格資料，當p = 1時代表灰度值影像，當 P = 3代表彩色影像。綜合考量影像本身的空間 資訊和色彩(或灰度值)資訊，可組成一個p+2維 的向量x = (xs, xr)，其中xs表示網格點的座標，
xr表示該網格點的p維光譜(或彩色域)的特徵向 量。由於兩者計量方式不同，為利後續演算分 析需作常態化處理，使之具有歐幾里得特徵空 間，因而將多維核函數定義如第(6)式：
(6) 式中，xs為特徵向量的空間部分，xr為特 徵向量顏色的部分；k (x)在空間及顏色部分都 使用相同的核函數。核函數決定了平均值調整 演算法的收斂性質。其中hs、hr分別代表空間
及色彩空間的移動式窗半徑或稱之為搜尋寬帶 (band width)，C為一常態化常數。因此寬帶參 數(hs、hr)就成為平均值調整分割過程中必須由 使用者決定的輸入參數，其對演算的收斂速度 與準確性有很大的影響。
以平均值調整為基礎的影像分割程序，主 要分成3個步驟：1.影像平均值濾波、2.以濾波 結果進行的影像分割(聚集)，以及3.影像分割結 果的最適化。以數學式說明如下：
以xi和zi (i = 1, …, n) 分別表示原始和濾波 後的影像像元的特徵向量，平均值調整濾波算 法可描述成：
針對每一個像元點
(1) 初始化j = 1，並使yi,1 = xi，其中 {y: (y-x)T(y-x)≦h2}
(2) 運用第(5)式迭代計算yi,j+1直到收斂，記 取收斂後的的值為yi,c = (ysi,c，yri,c) (3) 賦值(assign) zi = (xsi，yri,c)
此一迭代最終成果可視為一個動態移動視 窗與影像卷積(convolution)運算的過程，透過此 一迭代過程，每個像元點皆和一個潛在的收斂 聚集點相聯繫。
因 而 基 於 上 述 平 均 值 調 整 的 濾 波 演 算 法 ， 可 容 易 地 拓 展 以 平 均 值 調 整 的 影 像 分 割 法 。 設Li為 第i個像元的標號，影像分割過程 描述如下；
Fig. 2. Mean shift procedure. Starting at data point xi, the mean shift procedure is run to find stationary points of the density function. Superscripts denote the mean shift iteration, the shaded and gray dots denote input data points and successive window centers, respectively, and the dotted circles denote the density estimation windows, (modified from Depanis, 2005).
(1) 對所有像元進行平均值濾波處理，獲得 所有像元點的收斂值zi(i = 1,…, n)。 (2) 對特徵空間中zi = (xsi，yri,c)，將所有
zsi≦hs且zri≦hr歸為一類，獲得M個聚集 組的中心：M{Cp}p = 1, .. M。
(3) 為每個像元 i (i = 1,….n) 依據收斂聚集 點 所 屬 聚 集 組 中 心 賦 予 物 件 標 號 ， 即 Li = {P|Ziϵ Cp}，所有像元標記結束。
最適化：限制個別聚集區域所包含的像元 個數須≧Q，可用來決定遙測影像分割物件的最 小面積單元(minimum map unit; MMU)，可去 除雜訊(Boukir et al. 2012)。因而平均值調整影 像分割法所需輸入的限定參數(本文稱門檻值)有 3個，即hs、hr與Q。
(三) NDVI影像分割與崩裸地萃取
在3期FS2影像崩裸地自動萃取程序中，先 以平均值調整影像分割法，依據演算法必須輸 入的3個參數(如上節所述的hs、hr與Q)，選定適 合的空間及NDVI的灰度尺度門檻值，及分割區 塊最小單元面積(以像元數計算)，進行影像分 割。再依據分割區塊(object)的NDVI平均值，
依據直方圖分割法，選取最佳門檻值，將植生 覆蓋區塊與非植生覆蓋區塊分離出來，最後將 彼此相鄰接的非植生區塊合併成連續的區塊。
之後在區塊屬性選取上，先計算可以包含個別 區塊的橢圓的長軸與短軸的長度，以長軸與短 軸的比值當作區塊的形狀方形比指數，據以濾 除指標值小、區塊分割形狀接近方形及圓形、
面積很小的非植生區塊(將其歸類為分類不確定 之雜訊)。再以5 m空間解析度的數值高程資料 計算試驗區域的坡度主題圖，進而套疊計算每 一非植生覆蓋區塊(物件)的平均坡度，將平均坡 度低於選定的門檻值的區塊分離出來，其屬於 地形平坦地區之非植生覆蓋區塊，依此可將本 試驗區的屬於原住民村落的區塊分離出來，最 後得到3個時期試區的崩裸地分布主題圖。整個 影像分割及崩裸地萃取流程圖如Fig. 3。
(四) 萃取崩裸地的精度評估
3期崩裸地分布主題圖萃取結果準確度評
估，係以系統樣區取樣法，以每一矩形樣區大 小為2.2×2.2 km，使用起始點逢機的系統取樣 共取得8個樣區，總取樣率為12.2%。在每個樣 區內依據同一時期之全色銳化融合影像，在地 理資訊系統上，採用人為數化方式，為盡量達 到數化結果的客觀性，全部數化過程只由一位 有豐富判釋經驗的研究人員擔任數化工作。於 每 一 期 數 化 結 果 由 另 一 位 資 深 研 究 人 員 以3D 展示進行檢核與挑錯工作，以減低人為數化引 起 的 疏 失 ， 並 能 確 保3期影像崩裸地數化結果 的一致性與品質。依此數化所有樣區內崩裸地 塊，以其當作地真(ground truth)值，將各期自 動萃取的崩裸地區塊與各期人工數化崩裸地相 互比對。在準確度評估上，則以建立混淆矩陣 (confusion matrix or error matrix)、計算Kappa 係數(Kappa coefficient)，以及使用者精度(user accuracy)等，依此進行自動萃取崩裸地的準確 度評估。
(五) 崩裸地變遷分析
將上述步驟所得的3期崩裸地分布主題圖，
以影像主題圖變遷套疊分析法，得到兩個時段 (2006年9月28日至2009年8月19日、2009年8月 19日至2011年8月17日)的崩裸地與非崩裸地的 變遷分析圖，以及兩年期崩裸地面積轉移矩陣 (area transition matrix)。從面積轉移矩陣可得 知兩個時段有關崩裸地與非崩裸地面積的消長 變動及彼此的轉換率。由於3期影像時間間隔不 一，在轉換率計算與比較上，則使用兩年期時 間長度為基準，加權計算兩個時段的轉換率，
以利於相互比較。
(六) 崩裸地向源侵蝕區位分析
本研究所稱的向源侵蝕(headward erosion) 是指在河流的上游及山區的河流，由於河床的 縱坡比降大，因此流速也大、搬運力強，對河 底岩層的磨蝕及沖蝕強烈，故以下蝕為主。這 種下蝕作用使得河床不斷加深，時常形成深而 窄的谷溝，因而不斷的向上游延伸的現象；向 源侵蝕一方面使河流延長，另一方面也使河流 上游的分支越來越多(Pan 2007)。以往的國內
Fig. 3. Flow chart of the mean shift image segmentation and object extraction in this study.
研究(Chen and Ferng 2009, Lin et al. 2013)指出 莫拉克颱風發生後，位在暴雨下的山區林地增 加的崩塌地有向源侵蝕的現象。因此，本研究 以林業試驗所六龜試驗林轄區為範圍(包含在本 研究試驗區與範圍內)，使用試區5 m數值高程 資料，以Strahler’s河川級序萃取原理，先求取 河川流向，後以網格點個累積流量作為河川支 流萃取的門檻值；本研究選取兩個合適的門檻 值，將兩個萃取的向源河川支流套疊在2011年 新增崩塌地，以證明2011年新增的崩塌區位是
否位於向源侵蝕區域。
結 果
一、NDVI影像分割與崩裸地萃取
用於本研究的3期NDVI平均值調整影像分 割的3個參數經反覆測試後，3期影像的都用同樣 的參數門檻值，即hs = 15、hr = 15.5與Q = 64。
其中Q = 64個像元，每個像元面積為4 m2，因 此設定本研究NDVI影像分割後最小區塊面積
(MMU)為0.0256 ha。有關hs輸入的值域為[1, 20]
的整數值，代表影像特徵在空間尺度的近鄰程 度，其值越小搜尋半徑越大，輸出結果呈現平滑 化；使用較高的值對於較小的圖徵較具聚集效 果，所得成果較詳細，本研究為得到較詳細的空 間分割結果，選用15當作空間尺度的門檻值。有 關hr輸入的值域為[1.0, 20.0]的實數值，代表光 譜域的區分程度，其值越高代表分割區塊間光譜 差異性要求越高。本研究要求NDVI的灰度值差 異在中上程度，故選擇15.5當作灰度差異的門檻 值。經由本平均值調整影像分割後，可輸出一個 8位元(bit)的灰度值區塊分割影像，相當於最多 儲存256個聚集組。由於輸入的NDVI灰度值域為 [-1.0, 1.0]，因此影像分割後依據NDVI的大小，
延展成[0, 255]的灰度影像。因此可以在地理資 訊系統上，針對灰度直方圖分割法，採用natural break分級方法(本方法以符合分類組之組間變異 最大，組內變異最小的統計原理進行分級)，將 分級影像與全色銳化融合影像套疊，以交談方式 尋得分離崩裸地與非崩裸地的分割門檻值，低於
某灰度值的分割區塊屬於崩裸地。經由分析結果 2006年、2009年及2011年的們門檻值略有差距，
分別為189、184及188，若小於等於此門檻值的 區塊屬於非植生覆蓋區塊。
以Fig. 4說明本研究運用平均值調整影像 分割NDVI影像，萃取非植生覆蓋區塊的過程。
圖中以2009年8月19日拍攝的FS2全色銳化融合 影像的NDVI局部區域(右圖中黑色框選區域)放 大顯示為例加以說明。從Fig. 4B可看見NDVI 經由本影像分割法，採取上述適當的門檻值分 割後的NDVI均質的區塊，屬於非植生覆蓋區塊 (紫色區域)邊界相當吻合。再經由選取適當的區 塊NDVI平均值小於等於184的灰度門檻值，可 萃取出非植生覆蓋區塊。將相鄰區連續的區塊 加以合併後，經濾除區塊的形狀方形指數小於 等於1.1的面積很小的雜訊區塊後。再濾除區塊 平均坡度小於等於16.5度的區塊，可將座落在 本研究試驗區域的茂林、萬山及多納原住民山 村聚落區塊分離，得到各期的崩裸地主題圖。
Fig. 5顯示以本分析程序完成的3期崩裸地(以紫
Fig. 4 (A). FS2 image acquired on 19 August 2009 as an example and sub-figures to explain the procedure of mean shift image segmentation and object extraction in this study. Using pan-sharpened image (left image) to derive the normalize difference vegetation index (NDVI), the area selected by the black box is magnified in the next sub-figure (B).
Value High:1.0 Low:-1.0
Fig. 4 (B). Top image showing part of the image magnified from (A), middle image showing the segmented image with the mean shift algorithm, and bottom image showing non- vegetation objects using the NDVI mean threshold.
Fig. 5. Thematic map of landslide and bare areas (in purple) on each date output from the mean shift segmentation and object-extraction procedure in this study.
色表示)萃取的主題圖。
二、崩裸地萃取精度評估
Fig. 6顯示在試驗區域以系統取樣法設定的 8個2.2×2.2 km樣區分布情形，以3期影像的全 色銳化後高解析度影像為依據，藉由人為目視 判釋方式，在地理資訊系統在螢幕數化的崩裸 地塊，當作地真值(如Fig. 7)。將其與各期自動 萃取的崩裸地套疊比對。以混淆矩陣進行崩裸 地自動萃取的精度評估，所得結果如Table 1。
其中Table 1.(a)以2009年為例顯示其混淆矩陣；
將3期崩裸地萃取之各種精度評估結果整理後如 Table 1. (b)。分析結果得知，3期綜合分類準確 度(overall accuracy)皆達95%以上，Kappa係數 皆達0.85以上，其中使用者精度亦皆達87%以 上。顯見本研究採用平均值調整影像分割法及
簡單設定形狀與坡度門檻值的區塊過濾程序，
具有可應用性及在研究方法具有再現性。
另以Table 2與Fig. 8顯示3期崩裸地面積占 全研究試驗區域的面積比例，可得知以2009年 8月莫拉克颱風過後，造成本區崩裸地面積急 遽增加，2009年比2006年增加了3598.8 ha，增 加率為12.2%；到了2011年崩裸地面積卻有下 降，比2009年減少1039.3 ha，減少率為3.5%。
三、崩裸地面積變遷分析
Fig. 9顯示3期崩裸地分布主題圖，以影像 主題圖變遷套疊分析法，得到兩個時段(從2006 年9月28日至2009年8月19日，以及從2009年8 月19日至2011年8月17日)的崩裸地與非崩裸地 的變遷分析圖。每一時段的變遷有4種狀況，以 圖例表示於圖中；其中空白色代表有植生覆蓋
Fig. 6. FS2 image acquired on 17 August 2011 as an example to display distributions of 8 (2.2 by 2.2 km) plots determined by the system sampling method. Shaded red patches in each plot represent landslide and bare areas, which were output from the mean shift segmentation and object-extraction procedure in this study.
Fig. 7. Distributions of 8 plots determined by the system sampling method. Shaded red patches in each plot represent landslide and bare areas artificially interpreted and digitized through a geographic information system, and them were used for ground truthing as an accuracy assessment.
Table 1. (A) Accuracy assessment: confusion matrix on 19 August 2009
Class Manual digitization (%)
Non-LSAB LSAB Sum (%)
Non-LSAB 89.86 3.23 20.78
LSAB 10.14 96.77 79.22
Sum (%) 100.00 100.00 100.00
LSAB, landslide and bare area.
Table 1. (B) Accuracy assessment: summary table of different accuracies on 3 dates
Date 2006.09.28 2009.08.19 2011.8.17
Accuracy Producer User Producer User Producer User
Class (%) (%) (%) (%) (%) (%)
Non-LSAB 87.17 87.11 89.86 87.62 98.22 97.30
LSAB 99.18 99.18 96.77 97.41 87.69 91.61
Overall (%) 98.45 95.37 96.31
Kappa 0.863 0.858 0.873
LSAB, landslide and bare areas; Kappa: Kappa coefficient; overall: Overall accuracy.
Table 2. Area occupancies to the entire study area of landslide and bare areas on 3 dates Date / Area LSAB (ha) Rate (%) Non-LSAB (ha) Rate (%) 2006.09.28 1404.48 4.77 28,022.50 95.23 2009.08.19 5003.26 17.00 24,422.55 82.99 2011.08.17 3963.94 13.47 25,461.82 86.53
*Total area of the experimental region: 29,426.98 ha; LSAB, landslide and bare area.
和非崩裸區(如山村聚落地區)，紅色表示崩裸地 增加的區域，綠色表示原為崩裸地轉為有植生 覆蓋的植生復育區，藍色表示兩期都是崩裸地 無變動地區。從3期崩裸地變遷面積分析，可得 到兩個時段的崩裸地與非崩裸地的面積轉移矩 陣(如Table 3)。表中以兩時段的衛星影像拍攝 時間間隔的天數為加權基準，以間隔兩年為基 期，經調整後的非崩裸地與崩裸地之間的面積 轉換率(如Table 3最右欄的兩年變化率)。以直 方圖表示兩時段的非崩裸地與崩裸地之間的面 積轉移率(如Fig. 10)。
從 崩 塌 地 變 遷 分 析 結 果(Table 3及 Fig.
10)，得知2009年由於莫拉克颱風之超大型暴 雨，使得試區的非崩裸地轉移成崩裸地的轉移 率是後續兩年期間的3.56倍；而2009年至2011 年期間，由崩裸地轉移為非崩裸地的復育率是
2006年至2009年的18.69倍，且遠大於崩裸地的 轉移率，其植生復育率是崩裸率的2.44倍。顯 見本區有明顯植生(經現地勘查，多為草本植物 及先驅樹種的幼樹所覆蓋)復育發生，因此崩裸 狀況不再繼續增加。
四、向源侵蝕區位分析
本研究以Strahler’s河川級序萃取原理，以 1500個網格點及1000個網格點，兩個累積流量 作為河川支流萃取的門檻值，將兩個河川支流 套疊在2011年新增崩塌地，可明顯發現2011年 大部分的新增的崩塌往一級河流上游發展，位 於向源侵蝕區位(如Fig. 11)。以Fig. 11A說明 之，其底圖為2011年FS2全色銳化融合影像，紅 色區塊為2009年至2011年增加的崩裸地，深藍 色粗線為常流河，淺藍色支流為以1500個網格 Fig. 8. Bar chart showing the ratio of the total area of landslide and bare areas (red-colored portions) on 3 dates.
Fig. 9. Thematic change maps of landslide and bare areas between 2 periods, from 28 September 2006 to 19 August 2009, and from 19 August 2009 to 17 August 2011. There are 4 categories in the change maps: the white color-shaded areas represent vegetation-covered areas, the red ones represent increased landslide and bare areas, the green ones represent vegetation-recovered areas, and blue ones represent common existing landslide and bare areas between the 2 dates. The left map shows the change result of the first time period, and right one shows the change result of the second period.
點為門檻決定的向源河段、黃色支流為以1000 個網格點為門檻決定且增長的向源河段(可代表 河川向源延伸的方向)。Fig. 11B圖中灰色區區塊 代表2009年莫拉克颱風後的崩塌地區塊，紅色 區塊為2011年新增的崩塌地區塊；藍色粗線為 常流河，藍色細線代表以1500個網格點為門檻 決定的向源河段、紫色細線代表為以1000個網 格點為門檻決定且增長的向源河段(可代表河川 向源延伸的方向)。此外，國內在阿里山集水區 的相關研究(Lin et al. 2013)，也證實了向源侵蝕 的論點。
討 論
一、地形效應的影響
本研究在研究初期從中央大學遙測中心蒐
集了本試驗區從2005年至2016年共7期幾乎無雲 的影像，從中選取莫拉克颱風前後各隔兩年接近 夏季拍攝的FS2影像。其中2007年無夏季可用的 影像，只有當年1月15日拍攝影像品質較佳，由 於拍攝時間處於冬季太陽高度角偏低(36.7°)，在 試驗研究區域造成過多的山區地形陰影區。經本 期影像的NDVI影像分析後，發現許多陰影區仍 然無法消除，為了更客觀進行3個時期崩裸地的 比較，如使用本期影像恐會造過多的誤差。最後 往後再推一年，選擇與2009年8月拍攝月日相近 的影像，找到2006年9月28日拍攝的FS2影像(太 陽高度角54.5°)。發現該影像的NDVI影像，可減 除地形陰影的影響，同時與其他兩期太陽高度角 差異亦不大(2009年8月19日太陽高度角為55.8°、
2011年8月17日太陽高度角為60.1°)。因此本研 究3期影像的NDVI經由平均值調整影像分割法，
Table 3. Area transition matrixes of landslide and bare areas between 2 periods
2006.09.28 2009.08.19 Area (ha) Rate (%) Rate in 2 yr. (%)*
Non-LSAB Non-LSAB 24,514.46 82.64 57.13
Non-LSAB LSAB 3736.98 12.60 8.71
LSAB Non-LSAB 138.76 0.47 0.32
LSAB LSAB 1265.26 4.27 2.95
2009.08.19 2011.08.17 Area (ha) Rate (%) Rate in 2 yr. (%)*
Non-LSAB Non-LSAB 23,934.63 80.69 80.91
Non-LSAB LSAB 724.53 2.44 2.45
LSAB Non-LSAB 1763.85 5.95 5.98
LSAB LSAB 3239.49 5.95 5.98
* time weight: 2006.09.28 → 2009.08.19 1056 d or 2.89 yr;
2009.08.19 → 2011.08.17 728 d or 1.99 yr.
LSAB, landslide and bare area.
Fig. 10. Bar chart indicating transition rates (red-colored portions) of Landslide and bare areas between the 2 periods, from 28 September 2006 to 19 August 2009 and from 19 August 2009 to 17 August 2011.
萃取出來的崩裸地區塊可避免地形陰影的影響。
過去在山區進行的物件導向的崩裸地萃取，使用 不同的影像分割法，難免都會受到地形陰影的影
響，必需將受地形陰影影響的區塊以遮罩方式隔 離(mask out)，不進一步分析(Liu et al. 2011, Lin et al. 2013)。其在進行多期影像比較時，會濾除
Fig. 11 (A). Background picture which is a part of the FS2 pan-sharpened image in 2011, with red parches indicating increased landslide and bare areas from 2009 to 2011. Dark- blue lines are perennial flow streams, light-blue lines stand for headward tributaries extracted by a threshold with 1500 pixels cumulative flow, and short yellow lines show the direction of extended reaches which stand for further headward tributaries extracted by a threshold with 1000 pixels cumulative flow.
更多的陰影區，需使用共同非陰影區進行比較，
在實際應用上有其缺點。
二、同質性影像分割方法的探討
本研究採用的平均值調整影像分割法，在 影像分類法上可歸類為非監督式聚集分類法。
有 關 聚 集 分 類 方 法 上 可 再 分 為3大類型1.分區 式方法(partitioning methods)、2.階層式方法 (hierarchical methods)及3.以密度為基礎的方法 (density-based methods) (Boukir et al. 2012)。在 遙測影像分類領域，分區式方法較為人熟知的如 K-means和對其加以改良的Isodata演算法。階層 式方法可分為逐步分離和逐步聚合兩類，其中逐 步分離法以從上到下(top-down approach)的Ward 統計分析法為代表(Chang and Chen 2006, Kun et al. 2014)。本研究的平均值調整影像分割法則屬 於第3.類方法，與其他兩類方法比較具有較多優 勢。其中分區式聚集方法，必須預先決定期望的 聚集組數目及聚集組平均值收斂的預估迭代次
數，同時還要設定聚集演算的適當的起始點，不 同起始點所得聚集結果會有差異。此外依據多維 影像特性空間上由成員聚集的多個聚集組呈球狀 分布，聚集組間在特徵空間上會有重疊的疑慮。
另一方面，階層式方法一旦成員被歸類成同一類 別後就固定，不可能再歸到另一類別。依據多階 層分類樹結構，判定可分離的聚集組數目，過於 主觀，也由於演算方法的限制對影像的雜訊較為 敏感，多維演算時較分區式方法難行，一般分類 結果不如分區式方法。
以密度為基礎的方法的概念，主要是透過 多維特徵空間以核密度估計方法計算在一定搜 尋半徑(門檻值)的成員；以移動式窗方式依據 機率密度函數的梯度方向快速移動，連續加入 聚集成組直到收斂。方法上先透過平滑化過濾 雜訊、有不需預先設定期望的聚集組數目、以 無母數方法對針對個案的輸入資料進行估計，
以及聚集組的形狀不受限等優點。本研究採用 的平均值調整聚集演算法，方法上只需決定移
Fig. 11 (B). Gray patches indicating landslide and bare areas after Typhoon Morakot. Red patches indicate increased landslide and bare areas from 2009 to 2011. Blue thick lines are perennial flow streams, fine blue lines stand for headward tributaries extracted by a threshold with 1500 pixels cumulative flow, and short purple fine lines show the direction of extended reaches which stand for further headward tributaries extracted by a threshold with 1000 pixels cumulative flow.
動視窗的半徑(空間域及光譜域的門檻值)，適用 於高解析度衛星影像的分割，已成功的應用在 彩色影像及灰度影像的影像分割與分類上。缺 點是計算繁複，面對大量的遙測多譜影像電腦 計算時間長，但近年來投入改良本演算法的相 關研究頗多，大多用來改善本方法電腦計算的 速度(Christoudias et al. 2002, Zhang et al. 2005, Huang and Zhang 2008, Tang et al. 2009, Zhou et al. 2012)。
三、改善崩裸地萃取精度評估
本研究所定義的崩裸地有別於一般所稱的 崩塌地(landslide)，是指從面積大於等於0.0256 ha非植生覆蓋區塊、以區塊方形指標濾除雜訊 後、區塊平均坡度大於16.5度，以及不屬於山 村聚落的區塊(參考Fig. 3流程圖)。因此崩裸地 事實上包含因各次風災豪雨、暴雨後，河道變 寬、崩塌土石滑動後的堆積區域，或某些非自 然災害引起的少數裸露地。因考量不同時期受 多次自然干擾後的崩裸狀況不同，以及在各期
FS2相同感測時間的試驗研究區的航空照片缺 乏下。本研究為了兼顧時間及空間尺度的一致 性，所以使用同一時間的FS2全色銳化影像，
在系統取樣區內以人工數化的崩裸地多邊形當 作地真，再與各期在8個系統樣區內的自動萃取 崩裸地比對，進行分類準確度評估。如此可避 免過去的相關研究都以時間相近、空間解析度 較高的航照影像當作地真(Liu et al. 2011, Lin et al. 2013)，有時間與空間尺度不一致，而進行相 互比較的缺點。
結 論
本 研 究 試 驗 區 域 涵 蓋 將 近30,000 ha的林 地，選用福衛二號2006年、2009年及2011年3期 影像，經影像輻射校正前期處理、及全色影像融 合處理增進空間解析度，不改變多譜波段的特 性，以提高後續分析的精度。以各期影像NDVI 影像，應用平均值調整影像分割法，及簡易物件 形狀指標及坡度法則篩選程序，自動萃取3期試
位置，許多位在原有崩塌地往上游之向源侵蝕部 位。由於本試驗區林相、土地利用型及地物空間 異質性相對單純，因此能獲良好的分析結果，期 望後續研究可以將此分析方法擴展於更崎嶇的地 形與更複雜林地類型的分類研究上。
謝 誌
本研究報告得以順利發表，要感謝2位匿名 專業審查委員悉心閱讀及審查內容，並提供寶 貴修正意見。本研究經費由農委會林試所森林 經營組科技計畫及林業發展計畫支援，在此一 併致謝。
引 用 文 獻
Blaschke T. 2010. Object based image analysis for remote sensing. ISPRS J Photo Remote Sens 65(1):2-16.
Blaschke T, Hay GJ, Kelly M, Lang S, Hofmann P, Addink E, Tiede D. 2014. Geo- graphic object-based image analysis-towards a new paradigm. ISPRS J Photo Remote Sens 87:180-91.
Blaschke T, Johansen K, Tiede D. 2011.
Object-based image analysis for vegetation mapping and monitoring. In: Weng Q, editor.
Advances in environmental remote sensing- sensors, algorithms, and applications. New York: CRC Press. p 241-71.
Bo S, Ding L, Li H, Di F, Zhiu C. 2009.
Mean shift based clustering analysis of
Chang LY, Chen CF. 2006. Road extraction based on watershed segmentation for high resolution QuickBird satellite images. J Photo Remote Sens 11(3):261-67. [in Chinese with English summary].
Chen SC, Ferng JW. 2009. Topographic properties of channel head location based on a digital elevation model: a case study in an up- stream basin of the Takejin River. J Chin Soil Water Conserv 40(2):123-34. [in Chinese with English summary].
Chen Y. 1995. Mean shift, mode seeking, and clustering. IEEE Trans Pattern anal and Machine Intel 17(8):790-9.
Christoudias CM, Georgescu B, Meer P.
2002. Synergism in low level vision. In: Pattern Recognition, 2002. Proceedings. 16th Intl Conf IEEE 4:150-5.
Comaniciu D, Meer P. 2002. Mean shift: a robust approach toward feature space analysis.
IEEE Trans Pattern anal and Machine Intel 24(5):603-19.
Comaniciu D, Ramesh V. 2003. Kernel-based object tracking. IEEE Trans Pattern Anal and Machine Intel 25(5):564-77.
Derpanis KG. 2005. Mean shift clustering.
Lecture Notes. Available at http://www.cse.
yorku.ca/~kosta/CompVis_Notes/mean_shift.
pdf. 15 August 2005.
Dey V, Zhang Y, Zhong M. 2010. A review on image segmentation techniques with remote sensing perspective. In: Wanger W, Sezkely B,
editors. ISPRIS, Vienna, Austria, 5 July, 2010.
IASPRS XXXVIII (Part 7A):31-42.
Fukunaga K, Hostetler L. 1975. The estima- tion of gradient of a density function, with applications. in pattern recognition. IEEE Trans Inform Theory 21(1):32-40.
Hayes DJ, Sader SA. 2001. Comparison of change-detection techniques for monitoring tropical forest clearing and vegetation regrowth in a time series. Photo Engin Remote Sens 67(9):1067-75.
Huang X, Zhang L. 2008. An adaptive mean- shift analysis approach for object extraction and classification from urban hyperspectral imagery. IEEE Trans Geosci Remote Sens 46(12):4173-85.
Hussain M, Chen D, Cheng A, Wei H, Stanley D. 2013. Change detection from remotely sensed images: from pixel-based to object-based approaches. ISPRS J Photo Remote Sens 80:91-106.
Kamagata N, Hara K, Mori M, Akamatsu Y, Li Y, Hoshino Y. 2008. Object-based classification of IKONOS data for vegetation mapping in central Japan. In: Blaschke T, Lang S, Hay GL, editors. Object-based image analysis-spatial concepts for knowledge- driven remote sensing applications. Berlin, Heidelberg: Springer-Verlag Press. p 459-75.
Kun FE, Jan JF, Shao YC, Li MY, Yeh KS, Chen LH. 2014. Object-based classification for detecting landslides using high resolution aerial images. J Photo Remote Sens 18(4):267-81.
[in Chinese with English summary].
Lin CY, Jhang SW, Fu GL. 2013. Mechanism of headwater landslides in Alisan Creek watershed. J Chinese Soil and Water Conserv 45(2):653-62. [in Chinese with English summary].
Lin EJ, Liu CC, Chang CH, Cheng IF, Ko MH. 2013. Using the Formosat-2 high spatial and temporal resolution multispectral image for
analysis and interpretation landslide disaster in Taiwan. J Photo Remote Sens 17(1):31-51.
[in Chinese with English summary].
Lin H, Mo DK, Xiong YJ, Sun H. 2006.
The method of segmentation high-resolution remote sensing images based on mean shift. J Center South For Univ. 26(4):85-8. [in simple Chinese with English summary].
Liu CC. 2006. Processing of FORMOSAT-2 daily revisit imagery for site surveillance. IEEE Trans Geosci Remote Sens 44(11):3206-14.
Liu CC, Shieh CL, Lin J C, Wu AM. 2011.
Classification of non-vegetated areas using Formosat-2 high spatiotemporal imagery: the case of Tseng-Wen Reservoir catchment area (Taiwan). Int J Remote Sens 32(23):8519-40.
Lyon JG, Yuan D, Lunetta RS, Elvidge CD.
1998. A change detection experiment using vegetation indices. Photo Engin Remote Sens 64(2):143-50.
Mohan SS, Leela S. 2013. Importance of mean shift in remote sensing segmentation.
IOSR J Comp Engin 14(is6):80-3.
Nathalie P. 2013. The normalized difference vegetation index. New York: Oxford University Press. 194 p.
Neubert M, Herold H, Meinel G. 2008.
Assessing image segmentation quality con- cepts, methods and application. In: Blaschke T, Lang S, Hay GL, editors. Object-based image analysis-spatial concepts for knowledge- driven remote sensing applications. Berlin, Heidelberg: Springer-Verlag Press. p 769-84.
Pan KL, editor. 2007. Slope geology analysis and effective disaster prevention. Taipei, Taiwan: Technology-Book Press. 279 p.
Schiewe J. 2002. Segmentation of high- resolution remotely sensed data-concepts, applications and problems. Int Arch Photo Remote Sens Spat Inform Sci 34(4):380-5.
Schroeder TA, Cohen WB, Song C, Canty MJ, Yang Z. 2006. Radiometric correction
forest inventory polygons on high spatial reso- lution satellite imagery. For Chro 84(2):221-30.
Xiao PF, Feng XZ, editors. 2012. Segmen- tation and information extraction for high resolution remote sensing imagery. Beijing:
Science Press. p 314. [in simple Chinese with English summary].
Zhang K, Tang M, Kwok JT. 2005. Applying neighborhood consistency for fast clustering and kernel density estimation. In: 2005 IEEE Comp Scien Conf on Computer Vision and Pat- tern Recognition (CVPR’05). IEEE 2:1001-7.
Zhang Y. 2002. Problems in the fusion of commercial high-resolution satelitte as well as Landsat 7 images and initial solutions. Int
Beijing: Science Press. p 290. [in simple Chinese with English summary].
Zhou CH, Luo JC. 2009. High resolution satellite imagery in computing for geogra- phy. Beijing: Science Press. p 301. [in simple Chinese with English summary].
Zhou JX, Zhu JJ, Ma HY, Mei XM. 2012.
Improved mean-shift segmentation approach for remote sensing images. J Center South Univ. (Scien Tech) 43(2):620-5. [in simple Chinese with English summary].
Zhou W, Troy A. 2008. An object-oriented approach for analyzing and characterizing urban landscape at the parcel level. Int J Remote Sens 29(11):3119-35.
