中華民國一○六年十二月 DOI: 10.6574/JPRS
Journal of Photogrammetry and Remote Sensing
Volume 22 No.4 December 2017
Published by Chinese Society of Photogrammetry and Remote Sensing
Journal of Photogrammetry and Remote Sensing
發行人：史天元
出版者：中華民國航空測量及遙感探測學會 地址：台北市文山區羅斯福路五段 113 號三樓 電子信件：[email protected]
網址：http://www.csprs.org.tw
PUBLISHER: Peter Tian-Yuan Shih
PUBLISHED BY: Chinese Society of Photogrammetry and Remote Sensing
Address: 3F, No.113, Sec.5, Roosevelt Road, Taipei, Taiwan E-mail: [email protected]
Web Site: http://www.csprs.org.tw 總編輯：楊明德 EDITOR-IN-CHIEF:Ming-Der Yang
國立中興大學土木工程學系 Department of Civil Engineering, National Chung Hsing University 電 話：886-4-2285-2168 Tel: 886-4-2285-2168
電子信件：[email protected] E-Mail: [email protected] 編輯委員(依中文姓氏筆劃排列) EDITORIAL BOARD
王素芬 (國立彰化師範大) S. F. Wang (National Changhua University of Education) 王聖鐸 (國立臺灣師範大學) S. D. Wang (National Taiwan Normal University) 江凱偉 (國立成功大學) K. W. Chiang (National Cheng Kung University) 何宗儒 (國立臺灣海洋大學) C. R. Ho (National Taiwan Ocean University) 吳治達 (國立嘉義大學) C. D. Wu (National Chiayi University) 周天穎 (逢甲大學) T. Y. Chou (Feng Chia University) 林士淵 (國立政治大學) S. Y. Lin (National Chengchi University) 林老生 (國立政治大學) L. S. Lin (National Chengchi University) 林昭宏 (國立成功大學) C. H. Lin (National Cheng Kung University) 林唐煌 (國立中央大學) T. H. Lin (National Central University) 邱式鴻 (國立政治大學) S. H. Chio (National Chengchi University) 洪榮宏 (國立成功大學) J. H. Hong (National Cheng Kung University) 徐百輝 (國立臺灣大學) P. H. Hsu (National Taiwan University) 張中白 (國立中央大學) C. P. Chang (National Central University) 張智安 (國立交通大學) T. A. Teo (National Chiao Tung University) 莊永忠 (逢甲大學) Y. C. Chuang (Feng Chia University) 郭俊麟 (國立東華大學) C. L. Kuo (National Dong Hwa University)
陳朝圳 (國立屏東科技大學) C. T. Chen (National Pingtung University of Science and Technology) 曾義星 (國立成功大學) Y. H. Tseng (National Cheng Kung University)
黃金聰 (國立臺北大學) J. T. Hwang (National Taipei University) 詹進發 (國立政治大學) J. F. Jan (National Chengchi University) 雷祖強 (逢甲大學) T. C. Lei (Feng Chia University) 趙鍵哲 (國立臺灣大學) J. J. Jaw (National Taiwan University) 蔡富安 (國立中央大學) F. H. Tsai (National Central University) 蔡榮得 (國立中興大學) J. D. Tsai (National Chung Hsing University)
封面照片說明 About the Cover
西藏地區衛星影像雲偵測成果比較。左欄為原始影像，接續自左至右依次為李冠毅與林昭宏所提之方 法、Fmask(預設門檻值)、Fmask(最佳門檻值)、參考影像(人眼分辨)等方法辨識之成果。分類目標物包含雲
（紅）、水（藍）、雪（綠）、陰影（深紅）、其他（白）。
DOI：10.6574/JPRS.2017.22(4).1
1國立成功大學測量及空間資訊學系 碩士生 收到日期：民國 106 年 02 月 22 日
2國立成功大學測量及空間資訊學系 教授 修改日期：民國 106 年 06 月 23 日
*通訊作者, 電話: 0910336123 , E-mail: [email protected] 接受日期：民國 106 年 08 月 04 日
Landsat 8 衛星影像支持向量機雲偵測演算法
李冠毅
1*林昭宏
2摘要
光學遙感探測衛星影像中普遍有雲覆蓋地表的問題，同時限制處理影像的方法。多數先前研究中，
使用門檻值是最普遍的方法，然而，因地制宜的門檻值通常只適合該研究地區；地球環境不斷推移變化，
如持續用同一門檻值，勢必有失效情況發生。根據支持向量機雲偵測演算法，可避免上述問題。故本研 究利用統計模式建立分類基準，避免使用相對主觀的門檻值，依照雲和其他目標物對於不同波段之物理 特性相互調整出合適光譜特徵；紋理特徵部分則使用 Hotelling transform 再經過共現矩陣產生之紋理影像。
經實驗後，結果顯示本研究提出之方法其整體準確度介於 93%至 97%。
關鍵字：雲偵測、分類、支持向量機
1. 前言
雲偵測演算法是光學式衛星影像的重要議題，
因為雲可能造成許多相關影像分析和處理計算上 的誤差，故通常被視為需要被去除的目標物。因此，
光學式衛星影像的雲偵測法是許多研究前處理的 必要步驟。
從多數前人研究中發現 (Irish, 2000；Irish et al., 2006 ； Zhu and Woodcock, 2012 ； Zhu and Woodcock, 2014；Zhu et al., 2015) ，門檻值法是最 常被使用的方法，但決定門檻值的標準仍有許多問 題存在。地球環境因地球自轉和公轉而週期性地推 移變化，於其環境下若依照使用經驗訂定固定的門 檻值則較不恰當，例如：水會因為不同的溫度而有 三態變化，尤其高山或高緯度地區特別明顯，若使 用同一門檻值，勢必無法兼顧不同物質的物理變化。
故本研究提出一個不使用任何門檻值的方法進行 雲偵測，運用支持向量機相對客觀的統計模式，避 免主觀門檻值所造成的問題。
相較於僅有可見光和近紅外光的衛星影像，具 有短波近紅外光及熱輻射波段之衛星影像，如 Landsat 7 和 8 衛星影像，其優勢在於較能區分雪、
岩石或沙漠和建物等低溫或高反射率之地物。自動 雲 層 覆 蓋 評 估 法 (Automated Cloud Cover Assessment, ACCA) Irish (2000)及 Irish et al. (2006) 則是針對 Landsat 7 衛星影像設計，結合多項門檻 值搭配光譜特性和分類機制進行雲偵測，同理，
Fmask (Zhu and Woodcock, 2012)結合物件導向式 基礎推算雲陰影位置並可同時將其方法應用於 Landsat 4、7、8 或其他衛星影像上，其流程如圖 1 所示，經許多光譜門檻值篩選出雲、雲陰影和雪之 潛在區域再估計雲底層距離地面高度，疊代計算後 得雲、雲陰影和雪之最後確認像元。但 2015 年之 版本(Zhu et al., 2015)提到：「若水在寒冷環境及高 反射率的條件下，其中某一門檻值可能會失效。」
這也是門檻值法普遍難以克服的瓶頸之一。此外，
將同一方法應用於不同衛星影像上也為一大議題，
雖然 Fmask 可應用於多種衛星感測器，但仍需針 對該衛星影像之使用經驗設定不同門檻值，每一衛 星影像皆有特定之門檻值，故較難同時應用於其他 衛星影像，然而本研究提出之方法可避免類似問題，
雖不同區域仍須訓練不同模型進行分類，但不需依 照經驗設定門檻值，故較易使用同一作業流程進行 雲偵測。
雖本 針對雲，更 類目標，
水及其他 標物，因 如山及建 生之陰影 結果顯示 再由太陽 較可能獲 陰影位置
2.1 支
支持 決定一個 其中邊緣 短距離，
and Lin, 2 持向量(s 為非線性 可利用非 資料轉到 免此問題 加權總和 邊緣，此特 維度分類
遙測 研究使用 (Radial B
本研究著重雲偵 更能同時針對
因此本研究之 他。其中雲陰影
陰影可能由不 建物等，但本研 影，故此議題存 示，偵測雲陰影 陽、雲和雲陰影 獲得正確的雲 置勢必先做好
2. 研
支持向量機
持向量機(supp 個最大邊緣(m 緣代表兩組資 而邊緣的中央 2011)，因此最 support vector 性，最佳區分超 非線性之核函 到特徵域(featu 題。核函式的原 和，故可使用訓 特性讓支持向 類問題(Chen a 測影像運用於
用非線性之核 Basis Function
偵測，但仍期 對易與雲混淆 之目標物有四 影為雲偵測領 不同物體經光 研究欲分類出 存在排他性。
影要素時應先 影之光的直射 陰影位置。換 雲偵測。
研究方法
機
port vector m maximal marg 資料和超平面 央即為最佳區 最接近超平面 r)，如圖 2 所 超平面通常無 函式(kernel fu
ure space)上再 原理為使用支
訓練樣本(trai 向量機適合處 and Lin, 2006
分類上應為非 核函式，稱為 n, RBF)如下式
期望分類器能不 淆的雪和水作為 四大類：雲、雪 領域中一大分類 光線遮蔽所造成 出之目標為雲所 經本研究之實 先偵測出雲之位 射路徑條件下計 換句話說，欲知
法
machine)的原理 in)區分兩組資 (hyperplane)的 區分超平面(Ch 面的資料點稱為 所示。若資料來 無法完整運作
unction)把非線 再做分類，即可 支持向量所求出 ning sample)求 處理遙測領域中
；Chen, 2007 非線性問題，故
為徑向基底函 式：
圖 不只 為分 雪、
類目 成，
所產 實驗 位置，
計算 知雲
理為 資料，
的最 hang 為支 來源
，但 線性 可避 出的 求出 中高 7)。
故本 函式
平方 驗證
2.2
3 所 (20 頂 ( (rad 度部
200 不同 作為 同時 之目
階值 (int
和紅 於雲 別此 之數 水及
1 Fmask 流程
, exp
其中 方距離，而 k 證獲得。
2 光譜特
本研究使用 所示，並依照 15)將原始數位 (Top of Atm diance)再轉為 部分則轉為感 本研究參考 06)和 Fmask ( 同波段和目標 為區分雲、雪、
時採用盒鬚圖 目標物之數值 亮度特徵：雲 值，因此使用
ensity)區分像
其中 、
紅光波段之 T 雲和雪在此特 此二目標物，
數值；橫坐標為 及其他。
程圖
可視為兩 是任意設定之
特徵
用 Landsat 8 OL 照 Landsat 8
位數值(digita mosphere, TO 為 TOA 的反射 感測器(at-senso 考 ACCA (Iris (Zhu and Woo 標物的物理特性
、水和其他類別 圖分辨不同特 值成果，其特徵
雲普遍在衛星 用藍、綠和紅 像元亮度，如下
和 TOA 反射率 特徵皆為白色 如圖 4，縱坐 為欲分出之目
for 0
兩個向量間的 之參數，通常
LI 和 TIRS 之 Data Users H al number)轉
OA) 的 輻 射 射率(reflectanc
or)溫度。
sh, 2000；Iri odcock, 2012) 性，以其為依 別的支持向量 特徵下對於四
徵如下列所示 星影像中呈現 紅波段之平均 下式所示：
3
⁄ ...
分別代表藍
；I 代表輸入
，故於此特徵 坐標為亮度特 目標物，分別為
... (1)
歐基里德 常可由交叉
之資料如圖 Handbook
為大氣層 射 強 度 值 ce)；而溫
ish et al., )門檻值法 依據將公式 量機特徵，
種欲區分 示：
現較亮的灰 色彩強度
... (2)
藍光、綠光 入影像。由 徵中較難判 特徵所產生 為雲、雪、
雪特 (Normali 公式。依 性造成在 的差異，
特徵：欲區分 zed Differenc 依據物理特性 在此波段中反 其公式如下
分雲和雪，常 ce Snow Inde 性，雪會吸收短 反射率低於雲
：
圖 2
圖 3 Landsat
圖 常態化差異雪 ex, NDSI)為常
短波近紅外光
，利於凸顯雲
....
支持向量機示
at 7 和 Landsa
圖 4 亮度特徵 雪指標
常用的 光之特 雲和雪
. (3)
四 徵 故
依 段
示意圖
at 8 波段介紹
徵盒鬚圖 其中 四種目標物產生 徵下可明顯區分 故於此特徵中較
溫度特徵 依據此物理特 段難以區分的雲
代表短波近紅 生之數值如圖 分，但水和其他 較不能明顯區
：雲相對於其 性有助於分類 雲和雪，如圖
紅外光波段。
圖 5 所示，雲 他類有太多例 區分。
其他目標物之 類過程中判別 圖 6 所示，縱
。此特徵對於 雲和雪在此特 例外的情況，
之溫度較低，
別在可見光波 縱坐標為溫度 於 特
波 度
特徵數值 於此特徵 較低，有 輕易區分 寒冷 標物，使 和熱輻射
1
其中 和雪於此 和其他類 植物 量，故綠 吸收近紅 相對地，
率較低，
物有衰敗 率增加，
以上物理
其中 態化差異 Index, ND 如圖 6 所 物上仍有 有不同的 此三項特 屬副熱帶 果。
沙漠 外光和近
⁄ 因本
值，橫坐標為 徵中能夠明顯 有一定程度的 分雲和雪。
冷地表物特徵 使用雪特徵目標 射波段特性，
中 代表
此特徵較不相 類則有太多例外 物特徵：因葉 綠光反射率於 紅外光波段使 葉綠素吸收
，其公式為：
敗之現象發生 其公式為：
理特性，植物特
⁄
中 為近
異植披指數(N DVI)。以上三 所示，大致上 有些許不同。
的植物生長情 特徵。因觀察 帶及溫帶地區
漠地表特徵：
近紅外光反射
⁄ ...
本研究之研究
為欲分類之類別 顯區分，但熱輻 的影響，並不能
：針對凍原和 標物中短波近 其公式如下
...
表熱輻射波段 相同，有利於分
外情況發生不 葉綠素幾乎不吸 於此波段較高 使植物於近紅外 收大量紅光波段
⁄ 生，葉綠素會降
⁄ 特徵還包含第
近紅外光波段 Normalized Di 三種公式所計算 上數值成果類似 同時，因不同 況 發生，故本 察植物、卷雲及
，故圖 8 至
依據沙漠和岩 率較高之特性
...
究區域無沙漠地
別。雖四種目 輻射波段之解 能只依靠此特
和較寒冷之地 近紅外光波段
：
...
。如圖 7 所示 分類器計算，
不利於分類成 吸收綠光波段
；同時，葉綠 外光反射率很 段能量讓紅光
。同時，
降低導致綠光
。另外，
第三個公式：
...
，該公式又稱 fference Vege 算而得之特徵 似，但於不同 同植披地形可 本研究仍同時 及陰影之實驗 圖 11 並無雪
岩石具有短波 性，其公式為
...
地形，但此特 目標物 解析度 特徵就
地表目 段特性
. (4)
示，雲
，但水 成果。
段之能 綠素不 很高；
光反射
，若植 光反射
，總和
. (5)
稱為常 etation 徵數值 同目標 可能會 時使用 驗地區 雪之成
波近紅 為：
. (6) 特徵對
於 其
以 波 分 地 能 進
應 的 因 影 但 有
於區分雪有一定 其他類則易於此
卷雲特徵
以 表
波段，其特徵數 分效果。卷雲 地物則會影響 能有許多例外 進行分類。
陰影特徵 應較強；相對 的輻射擴散現 因此，在近紅 影的陰暗效應會 但陰影區域下 有其必要性。
定的幫助，如 此特徵中混淆
：卷雲波段為 示，針對卷雲 數值如圖 10 所
屬於較薄之雲 此特徵成果 情況會發生
：在較短波長 地，在較長波 象則較弱，造 外光波段和短 會更明顯。雖然
仍包含其他地
圖 5 雲特徵
圖 6 溫度
如圖 9 所示，
淆。
為 Landsat 8 最 雲特性設計出
所示，對於雲 雲層，若薄雲
，且變化萬千
，故並不能只
長的波段中，
波長的波段中 造成陰影較其 短波近紅外光 然本研究並無 地表物，故使
徵盒鬚圖
度特徵盒鬚圖
，但雲、水和
最新的波段，
出此近紅外光 雲有明顯的區 雲底下有其他 千的地球中可 只依靠此特徵
大氣散射效 中，陰影區域 其他地物暗，
光波段裡，陰 無分類陰影，
使用此特徵仍 和
光 區 他 可 徵
效 域
陰
仍
圖 7 寒冷
(a)
(b)
圖 8 植
冷地表物特徵盒
⁄
⁄
(c) NDVI 植物特徵盒鬚
盒鬚圖
鬚圖
2
薄 紋 應 如 (u 擷 上
圖
2.3 紋理特
薄雲之透 薄雲底下之地 紋理特性有必 應用許多共現矩 如：相關性(co
uniformity)、同 擷取不同波段 上相似的波段
圖 9 沙漠地表
圖 10 卷雲特
圖 11 陰影特
特徵
明特性造成其 表物，為考慮 要加入分類器 矩陣(co-occu orrelation)、對 同質性(homo
之紋理。然而 會產生許多高
表特徵盒鬚圖
特徵盒鬚圖
特徵盒鬚圖
其反射率也包 慮此特性所造 器特徵。於先 urrence matrix
對比性(contra ogeneity)和亂
而，若擷取可 高度相關的紋
包含了一部分 造成之影響，
先前研究中，
x)之紋理指標 act)、一致性 度(entropy)，
可見光等紋理 紋理特徵，故 分
標，
性
理 故
本研究使 式為y 形成 n 維 方差矩陣 其相對應 徵值所對 因分類器 紋理特性 性，其中 階值降為
∑ ∑
圖 12 Ho 值影
圖 13 共
使用 Hotelling
x ，
維向量x 陣 Cx計算而得 應影像之資訊 對應的三張影像 器中同性質的 性應較為均勻 中共現矩陣所使 為 32 階，其公
⁄ ...
otelling transfo
影像
共現矩陣成果所
圖 14
transform 解 假設 n 組向
, , … ， 得，因此特徵值 訊量則越多。本
像輸入共現矩 的特徵不宜太多 勻和一致，故取 使用的核罩窗 公式如下所示
...
orm 成果所對
所對應之前三
研究區域之
決這項問題，
量有相同的長
，而 A 矩陣是 值(eigenvalue)
本研究將前三 矩陣，如圖 12
多，且考慮到 取用之指標為 窗為 9 × 9 並 示：
...
對應之前三大特
三大特徵值影像
之 Landsat 8 衛
，其公 長度 K 是由協 )越大，
三大特 所示，
到雲的 為一致 並將灰
. (7)
特徵
像
3
8 20 月 由 中 仍 像 包 別 物 供 約 之
研 研
P
景 影
衛星影像。左
3.
3.1 研究區
本研究以美 OLI 和 TIRS 013 年 8 月 6 月 29 日並分別 由於氣候關係 中與雪相似，
仍視該目標物 像挑選三塊 10 包含各式目標物 別自三個研究 物(雲、雪、水 供給支持向量機 約為五萬像元 之詳細資料如表
表1 研 研究區域 研究區域 位置 Path-row
坐標系 景物特徵 農 影像獲得
日期 2
：美國加州。
實驗結
區域
美國加州、西 影像作為案例 6 日、2015 年 別包含農地、
，西藏地區之 但其物理特性 為水。實驗過 000 × 1000 像 物作為測試資 區域之衛星影 和其他)並盡量 機的訓練樣本
，如圖 15、1 表 1 所示。
研究區域之La 美國加州
38.9N, 121.7W
44/33 農地、植披
2013/08/06
中：西藏。
結果與討
西藏和臺灣地 例，影像取得 年 3 月 18 日和 雪地和海島之 之水已結為冰 性仍與雪不同 過程中分別從 像元的正方形 資料(testing d
影像中隨機挑 量包含各種樣 本(training dat
16、17 所示，
andsat 8影像資 西藏 31.7N,
81.7E
144/38 高山雪地
2015/03/18
右：臺灣
討論
地區之 Landsat 得日期分別為 和 2014 年 11 之地表特性。
冰且在可見光 同，故本研究 從三幅衛星影 形區域並盡量 ata)，同時分 挑選四類目標 樣本，作為提 ta)，其樣本數
，而研究區域
資料 臺灣 25.6N, 120.4E
117/43 海島地形
2014/11/29 t 為
光 究 影 量 分 標 提 數 域
圖15 美國 本研
圖16 西藏 究所
圖17 臺灣 究所
國加州之支持 研究選取之訓
藏之支持向量 所選取之訓練
灣之支持向量 所選取之訓練
持向量機訓練 訓練樣本
量機訓練樣本 練樣本
量機訓練樣本 練樣本
練樣本。黃色區
本。黃色區塊為
本。黃色區塊為 區塊為
為本研
為本研
3
行 用 類 像 稱 之 果 色 有 故
算 果 地 為 Fm 較 在
設 生 檻 和 則 情 結 提 性 料
結 誤 較 夾 究
3.2 實驗結
本研究將 行門檻值調整讓 用預設門檻值 類(最佳門檻值 像元為基準使用 稱為參考影像 之比較基準。表 果，分類類別 色分別為紅、綠 有雲陰影類別 故進行數值比較
表 2 加州 算法皆有考慮 果都有一定水 地區之第三張 為敏感，只要 mask 分類結 較寬，經調整門 在水體表現上
表 3 西藏 設門檻值，分 生，進而驗證 檻值可能會失效 和門檻值分類 則較為正常，
情況發生。此 結冰的情況發生 提出之方法仍 性，若將結成 料的一致性及應
表 4 臺灣 結果中發現 Fm 誤情況，經調 較複雜之地區 夾雜海洋和陸 究有些許分類錯
結果
Fmask 作為比 讓準確度提升 及調整過後取 值)。另外，為求
用 Photoshop
，以作為本研 表 2、3、4 為 為雲、雪、水 綠、藍、深紅
，而本研究並 較時將雲陰影 州地區主要觀察
光譜對於葉綠 準。另特別觀 影像可發現本 地表類似薄霧 結果中，使用預 門檻值後分類成
，兩者成果差 地區之影像中 類成果明顯有 水在寒冷環境 效(Zhu et al., 法的瓶頸；經 但仍有把雲陰 外，因氣候因 生，Fmask 將 將其分類為水 冰之水分類為 應用有所影響 灣地區實驗目標
mask 預設門檻 整後未見改善
，除雲之組成 地，另外，薄 錯誤情況發生
比較對象，因 升，故 Fmask 取最高準確度
求準確度分析 p 進行人眼判釋 研究所提出之方
為不同研究地 水、雲陰影和 紅和白色。由
並無針對雲陰 影併入其他類
察農地和植披 綠素之特性，
觀察薄雲分類 本研究之結果 霧的情況下皆 預設門檻值成 成果較接近於 差異不大。
中可以發現 F 有整體分類錯 境及高反射率 , 2015)之門檻 經門檻值調整 陰影分類為水 因素，水在此 將冰視為雪，
水，因考慮水 為雪，將會對 響。
標物為海島與 檻值在都市部 善。此地區第 成非常離散外 薄雲也較難界 生。
因 Fmask 可進 k 成果分為使 度之門檻值兩 析，本研究以 釋，其成果名 方法和 Fmask 地區之分類結 和其他，而顏 由於 Fmask 含 陰影做分類，
類別。
披地形，兩演 故此地區成 類成果，由此 果對於薄雲較 皆分類成雲；
成果明顯範圍 於參考影像。
Fmask 使用預 錯誤的情況發 率的條件下門 檻值設計缺點 整後分類成果 水的錯誤分類 此測試區域有
，而本研究所 水和雪的差異 對於後處理資
與都市，可於 部分有分類錯 第二張影像為 外，雲之中還 界定，故本研 進 使 兩 以 名 k 結 顏 含
演 成 此 較
圍
預 發 門 點 果 類 有 所 異 資
於 錯 為 還 研
Fmask 在最後的步驟中使用形態學(morphology) 之膨脹(dilation)與侵蝕(erosion)填補其分類成果之 小洞，故可發現其成果與參考影像比較起來較為粗 略。本研究所提出之方法則較不易有小洞需填補之 問題，故成果圖較相似於參考影像。另外，也可以 發現本研究所提出之方法中薄雲較容易分類錯誤。
由實驗結果可知，Fmask 經過調整門檻值後，其成 果皆有較佳之表現；而本研究雖需特別挑選訓練區 計算分類模型，但可同一步驟應用於不同地區。
3.3 成果檢核
本研究之成果與Fmask分別對於參考影像進行 檢 核 ， 將 參 考 影 像 作 為 地 真 資 料 並 使 用 準 確 率 (precision)和召回率(recall)計算F度量值(F-measure)，
其 公 式 分 別 為 ⁄ 、
⁄ 和 2 ⁄ ，其中 、
和 分別代表正確正例(true positive)、錯誤正例 (false positive)和錯誤負例(false negative)。整體準確 度 (overall accuracy) 則 定 義 為
⁄ ，其中 代表正
確 負 例 (true negative) 。 另 外 ， Kappa 定 義 為 / 1 ，其中 為觀測一致性機率，
為期望一致性機率。
對應於表2、3、4，研究地區檢核資料如表5、6、
7所示，本研究所提出之方法整體準確度大約於93%
至97%之間，皆高於Fmask使用最佳門檻值，其中準 確度浮動較大之原因可能在於SVM分類成果對於訓 練樣本挑選區域較為敏感，若訓練樣本無包含測試 地區之情況，分類準確度會直接受影響。而F度量值 方面，則易受薄雲所影響，如表6西藏地區第三張測 試影像中厚雲較薄雲量少，而薄雲因含有地物之反 射率，易干擾支持向量機之訓練模型。另因許多衛 星影像之光譜並未涵蓋短波紅外光和熱輻射波段，
若研究地區為雪地及高反射率之沙漠，可見光及紅 外光之反射率與厚雲接近，將不易分辨，以表3、6 之西藏地區為例，因選取之特徵包含雪特徵及寒冷 地表物特徵，本研究成果顯示雪的分類情況無明顯 分類錯誤，但Fmask預設門檻值成果如第1章所述：
「若水在寒冷環境及高反射率的條件下，其中某一 門檻值可能會失效。」(Zhu et al., 2015)，其成果有 明顯大範圍錯誤。
Fmask在本研究成果中，準確度皆較低於Zhu and Woodcock (2012)所顯示之成果，其可能原因在 於參考影像是以像元為基礎，又在Fmask中，為填補 最後生成雲之小洞，使用形態學(morphology)之膨脹 (dilation)與侵蝕(erosion)做後處理，因此Fmask之分 類成果皆高估參考影像之範圍，造成準確度上有所 影響，且Zhu and Woodcock(2012)之成果為平均整體 準確度，於某些地區準確度的確不高。但本研究之 研究地區數量不如Fmask，代表未能將本研究實驗於 各種可能的情況中。
為幫助判斷特徵之成效，本研究使用特徵選取 法，從原有的特徵中挑選最佳部分特徵，透過特徵 選取法而得之分數(F-score)鑑別能力較好的特徵，以 利瞭解本研究分類問題的因果關係。如圖 18 所示，
藍色、橘色和灰色線分別代表特徵對於不同研究地 區之分數，而黃色線則代表分數之加總，因此可知 亮度特徵為最佳特徵，其次為 NDVI，第三為 NIR，
但紋理特徵分數皆較低，其可能因共現矩陣之指標 和核罩窗所導致。考慮分類整體準確度是否會被某 一特徵所拖累，本研究使用分數前五名和前八名之 特徵進行訓練而得之模型測試同一資料，其整體準 確度皆不如使用所有特徵，故可知在分類器中，雖 特徵之特徵選取法分數較低，其仍有一定之幫助。
表 2 美
表 原始
原始
美國加州地區
表 3 西藏地區 始影像
始影像
區雲偵測成果
區雲偵測成果 本研究所提
本研究所提
果。雲、雪、水
。雲、雪、水 出之方法 Fm
出之方法 Fm
水、陰影和其
水、陰影和其 美國加州 mask (預設門
西藏 mask (預設門
其他之類別顏
其他之類別顏色 州
門檻值) Fmas
門檻值) Fmas
色分別為紅、
色分別為紅、
sk (最佳門檻值
sk (最佳門檻值
、綠、藍、深
綠、藍、深
值) 參
值) 參
深紅和白色
紅和白色 參考影像
參考影像
表
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
F-scor e
原始
表 4 臺灣地區
始影像
區雲偵測成果
Dataset I:
本研究所提
。雲、雪、水
圖1
California, USA
出之方法 Fm
水、陰影和其
18 特徵選取
Datas
臺灣 mask (預設門
其他之類別顏色
法成果
set II: Tibet, Chin
門檻值) Fmas
色分別為紅、
a Data
sk (最佳門檻值
綠、藍、深
aset III: Taiwan
值) 參
紅和白色
All
參考影像
表 5 美國 研究
美國加州
國加州地區檢 究所提出之方
原始影像種類
本研究所提出之方法
P(%
R(%
F(%
A(%
K(%
Fmask
( 預設
門檻值)
P(%
R(%
F(%
A(%
K(%
Fmask
( 最佳
門檻值)
P(%
R(%
F(%
A(%
K(%
檢核資料。以準 法、Fmask 使
其他
%) 98.12
%) 97.58
%) 97.85
%)
%)
%) 95.10
%) 73.72
%) 83.05
%)
%)
%) 85.72
%) 92.22
%) 88.85
%)
%)
準確度(P)、召 使用預設門檻
雲 水
90.66 98 94.70 92 92.63 95 96.75 91.76 58.26 70 89.70 66 70.64 68 75.72 55.73 84.76 66 57.86 78 68.77 72 81.81 57.24
召回率(R)、F 檻值和 Fmask 使
水 其他
8.19 98.08 2.73 98.21 5.38 98.14
0.42 96.09 6.80 60.62 8.56 74.34
6.75 83.76 8.37 95.12 2.09 89.08
F 度量值(F)、
使用最佳門檻
他 雲
8 90.67 1 89.62 4 90.14 96.88 88.30 9 38.89 2 87.46 4 53.84 64.35 35.42 6 72.61 2 32.09 8 44.51 80.48 33.69
整體準確度(
檻值
水 其
80.68 99 99.05 93 88.92 96
55.82 95 65.39 80 60.23 87
49.49 90 65.39 97 56.34 93
(A)和 Kappa
其他 雲
9.08 53.47 3.74 94.26 6.34 68.23 93.48 68.47 5.27 30.23 0.20 71.86 7.09 42.56 77.54 32.96 0.62 71.83 7.03 20.29 3.71 31.65 88.19 29.83
值(K)評估本
水
7 76.12 6 82.95 3 79.39 8
7
3 39.64 6 48.73 6 43.72 4
6
3 34.74 9 52.66 5 43.97 9
3
本
表 6 西藏 所提
西藏
藏地區檢核資 提出之方法、F
原始影像種類
本研究所提出之方法
P(%
R(%
F(%
A(%
K(%
Fmask
( 預設
門檻值)
P(%
R(%
F(%
A(%
K(%
Fmask
( 最佳
門檻值)
P(%
R(%
F(%
A(%
K(%
資料。以準確度 Fmask 使用預
其他 雲
%) 93.98 97
%) 99.71 73
%) 96.76 83
%)
%)
%) 0 12
%) 0 99
%) 0 22
%)
%)
%) 87.36 93
%) 87.87 44
%) 87.61 60
%)
%)
度(P)、召回率 預設門檻值和
雲 雪
7.78 91.51 9 .61 97.45 8 .99 94.39 9
94.56 87.35 2.50 0 9.59 0 2.21 0
13.61 3.66 .61 14.43 4.17 62.02 0.02 23.42
73.88 38.76
率(R)、F 度量 和 Fmask 使用最
水 其他
98.10 93.25 83.26 99.65 90.07 96.34
0 0 0 0 0 0
1.77 84.99 0.62 85.1 0.92 85.05
量值(F)、整體準 最佳門檻值
雲 雪
98.11 90.2 65.47 97.7 78.54 93.8 93.65 78.85 10.88 0 98.89 0 19.6 0 13.45
5.58 83.5 16.94 23.27 67.2 36.39 27.0 72.35 26.64
準確度(A) 和
水 其
7 94.83 93.
6 41.91 99.
6 58.13 96.
0 0 0 0 0 0
4 0.10 92.
3 0.68 94.
7 0.18 93.
和 Kappa 值(K
其他 雲
.95 99.79 9 .93 46.86 9 .85 63.77 9 94.2 62.88 0 5.47 0 99.60 0 10.37 9.20 4.39 .66 89.10 3 .73 39.47 2 .68 54.7 6
86.56 35.29
K)評估本研究
雪 水
92.19 89.62 97.41 51.21 94.73 65.18
1 8
0 0 0 0 0 0 0
9
3.50 0.21 23.61 2.59 6.10 0.38 6
9
究
表 7 臺灣 所提
臺灣
本研 使用支持 判斷標準
灣地區檢核資 提出之方法、F
原始影像種類
本研究所提出之方法
P(%) R(%) F(%) A(%) K(%
Fmask
( 預設
門檻值)
P(%) R(%) F(%) A(%) K(%
Fmask
( 最佳
門檻值)
P(%) R(%) F(%) A(%) K(%
4. 結
研究主要貢獻 持向量機代替 準，並應用多光
資料。以準確度 Fmask 使用預
其他 ) 96.27 ) 98.71 ) 97.48 )
)
) 94.75 ) 75.81 ) 84.23 )
)
) 97.19 ) 86.99 ) 91.81 )
)
論與建
在於提出免門 替經驗門檻值
光譜資訊於支
度(P)、召回率 預設門檻值和
雲 水
97.65 98 90.80 99 94.10 98 97.33 95.83 70.97 92 86.44 96 77.95 94 85.48 77.68 97.03 93 90.38 99 93.59 96 92.36 88.35
議
門檻值雲偵測 避免過於主觀 支持向量機，但
率(R)、F 度量 和 Fmask 使用最
水 其他
8.30 93.08 9.45 97.47 8.87 95.23
2.21 92.28 6.54 94.88 4.33 88.42
3.92 97.93 9.11 81.45 6.45 88.93
測法，
觀之 但針
對每 持向 影響 物的 考慮
量值(F)、整體準 最佳門檻值
雲 8 95.48 7 91.04 3 93.21 93.94 88.31 8 88.18 8 89.05 2 88.61 84.62 71.98 3 86.84 5 97.07 3 91.67 86.45 75.67
每張影像仍須 向量機上主要 響，大多僅使 的物理特性本 慮環境因素加
準確度(A) 和
水 其 90.10 96 73.46 91 80.93 93
92.62 96 31.23 70 46.71 81
88.11 95 42.29 77 57.15 58
須訓練分類模型 要探討不同核
用可見光及紅 本身無法突破雲 加入短波紅外光
和 Kappa 值(K
其他 雲 6.93 89.25 1.16 93.39 3.96 91.27 94.94 92.29 6.83 83.86 0.63 92.63 1.68 88.03 86.72 80.49 5.67 90.73 7.89 88.18 8.87 89.44 88.65 83.16
型。多數先前 核函式對於分 紅外光波段，
雲偵測問題，
光和熱輻射波
K)評估本研究
水 5 97.13 9 99.09 7 98.10 4
9
6 96.81 3 96.41 3 96.61 2
9
3 95.2 8 97.94 4 96.55 5
6
前研究於支 類結果的 導致目標 然而，若 波段，則可
究
突破障礙並增加分類準確度，例如：使用熱輻射波 段讓雲本身溫度較低之物理特性發揮作用，讓分類 器較能分辨其他目標物。
本研究所提出之方法為結合多光譜和紋理資 訊成為支持向量機之特徵，透過訓練樣本計算出分 類模型代替經驗門檻值，使分類成果更客觀與穩定，
其中光譜特徵包含 Landsat 8 OLI 和 TIRS 之紅、綠、
藍可見光波段、近紅外光波段、短波長近紅外光波 段和熱輻射波段，而紋理特徵則是使用 Hotelling transform 計算所有波段的大部分資訊後再使用灰 階共現矩陣獲得紋理資訊，雖紋理在特徵選取法結 果中分數較低，但經實驗證明仍需此特徵。研究區 域包含美國加州、西藏和臺灣之 Landsat 8 衛星影 像，挑選因素主要在於該地區分別涵蓋植披、高山 雪地和海島地形，而測試資料每張為 1000 平方像 元；訓練樣本約為五萬像元。實驗結果整體準確度 大約為 93%至 97%、Kappa 值大約落於 70%至 90%，
比 Fmask 更穩定，另 Fmask 於西藏地區中明顯出 現門檻值法失效的情況，經調整門檻值後準確度明 顯改善，故證實固定的經驗門檻值不一定適合不同 實驗資料，而本研究方法可用同一步驟應用於不同 地區，避免訂定過於主觀之門檻值。偵測薄雲為本 研究主要挑戰之一，因反射率涵蓋一部分地表資訊，
造成人眼判斷困難，且易被人為判斷因素所影響，
實驗成果顯示含有薄雲之影像易因訓練樣本是否 完整而影響分類結果，若訓練樣本完善則分類精度 可達 90%以上。
支持向量機能夠以多類別為基礎進行分類，故 未來應可加入雲陰影之類別，並應考慮太陽、衛星 及 雲 之 幾 何 關 係 並 以 數 值 高 程 模 型 (Digital Elevation Model, DEM)為地面計算基礎，使支持向 量機雲偵測演算法更為完整。
參考文獻
Chen, Y.W., and Lin, C.J., 2006. Combining SVMs with various feature selection strategies, In Feature Extraction, pp. 315-324, Springer Berlin Heidelberg.
Chen, G., 2007. Support vector machines for cloud detection over ice-snow areas, Geo-spatial Information Science, 10(2): 117-120.
Chang, C.C., and Lin, C.J., 2011. LIBSVM: A library for support vector machines, ACM Transactions on Intelligent Systems and Technology (TIST), 2(3): 27:1–27:27.
Irish, R., 2000. Landsat-7 automatic cloud cover assessment algorithms for multispectral, hyperspectral, and ultraspectral imagery, The International Society for Optical Engineering, 4049:348–355.
Irish, R.R., Barker, J.L., Goward, S.N., and Arvidson, T., 2006. Characterization of the landsat-7 ETM+ automated cloud-cover assessment (ACCA) algorithm, Photogrammetric Engineering and Remote Sensing, 72(10):
1179-1188.
USGS, 2015. LANDSAT 8 (L8) Data Users Handbook, Department of the Interior US Geological Survey, LSDS-1574 Version 1.0, 105.
Zhu, Z., and Woodcock, C.E., 2012. Object-based cloud and cloud shadow detection in landsat imagery, Remote Sensing of Environment, 118(2012):83–94.
Zhu, Z., and Woodcock, C.E., 2014. Automated cloud, cloud shadow, and snow detection in multitemporal Landsat data: An algorithm designed specifically for monitoring land cover change, Remote Sensing of Environment, 152(2014): 217-234.
Zhu, Z., Wang, S., and Woodcock, C.E., 2015.
Improvement and expansion of the Fmask algorithm: Cloud, cloud shadow, and snow detection for Landsats 4–7, 8, and Sentinel 2 images, Remote Sensing of Environment, 159(2015): 269-277.
1Master, Department of Geomatics, National Cheng-Kung University Received Date: Feb. 22, 2017
2Professor, Department of Geomatics, National Cheng-Kung University Revised Date: Jun. 23, 2017
* Corresponding Author,Tel: 886-910336123 , E-mail: [email protected] Accepted Date: Aug. 04, 2017
Cloud Detection Based on Support Vector Machine for Landsat 8 Imagery
Kuan-Yi Lee 1*
Chao-Hung Lin2
Abstract
Cloud covers are generally present in optical remote-sensing images, which limit the usage of acquired images. In previous studies, thresholding is a common and rapid method in cloud detection. However, a selected threshold is usually suitable for local study areas, and it may be failed in other cases. Besides, there are many exceptions to control, and the environment is changed dynamically. In this study, a threshold-free method based on Support Vector Machine (SVM) is proposed, which can avoid the abovementioned problems. A statistical model is adopted to detect clouds instead of a subjective thresholding-based method. According to the physical characteristics of clouds and other objects, the spectral features are appropriately designed for classification.
Spatial and temporal information are also important for satellite images. Consequently, co-occurrence matrix of the Hotelling transform is used in proposed method. Experiment results demonstrate the detection accuracy of the proposed method is about 93% to 97%.
Keywords：
Cloud Detection, Classification, Support Vector Machine
DOI 10.6574/JPRS.2017.22(4).2
1國立交通大學土木工程學系 教授 收到日期：民國 106 年 04 月 12 日
2國立交通大學土木工程學系 碩士 修改日期：民國 106 年 06 月 30 日
3國立交通大學土木工程學系 博士後研究員 接受日期：民國 106 年 08 月 15 日
*通訊作者, 電話: 03-5131487, E-mail: [email protected]
應用衛星影像的水線辨識於外傘頂洲的灘線變遷
張憲國
1*賴羿齊
2陳蔚瑋
3摘要
本文提出擷取衛星影像的水線，平移至灘線的方法，並探討外傘頂洲的灘線變遷。研究方法的影像 處理技巧包括，IHS 轉換法、NDWI、影像強化、形態學及 Canny 法。另外，修正 NAO99b 潮位數值模 式的潮位資料融合法，來解決外傘頂洲無觀測潮汐資料的問題。
自 1994 至 2015 年的灘線及陸上面積的長期變化分析後發現，外傘頂洲的海側灘線在北端以-30.5 m/yr，
而南端以-119 m/yr 的速度往內陸內縮；陸上面積被侵蝕速率近年有增加的趨勢，且長期的變化有夏季增 加，但冬季減少的季節特性。若以 2009~2015 年的陸地面積被侵蝕速率，可預測外傘頂洲的陸地將於 2060 年完全低於平均潮位而淪為潛沒沙洲。
關鍵字：外傘頂洲、灘線變遷、影像處理、水線辨識
1. 前言
外傘頂洲位於雲林嘉義海岸之外側，為臺灣目 前最大的離岸沙洲，此沙洲除可保護雲嘉地區海岸，
免於被外海波浪直接衝擊外，且與臺灣間形成廣大 的溼地，可庇護多種生物棲息於此。因此，外傘頂 洲具有保護雲嘉海岸，及提供優良的生態棲地的功 能。但是，外傘頂洲鄰近的濁水溪及北港溪近年來 的輸砂量減少，鄰近海埔地及離島工業區的開發改 變海岸地形，與超抽地下水導致地層下陷等複雜因 素，導致雲嘉南地區海岸一直存在有嚴重的侵蝕問 題。黃建維等(2015)指出外傘頂洲的陸上高度下降 及面積的減少，造成更嚴重的海岸侵蝕及衝擊當地 牡蠣養殖。根據蔡雨農等(2015)的研究，預估外傘 頂洲將於 2028 年完全沒入平均海平面以下。探討 外傘頂洲地形變遷特性及提出合適的保護工法已 是刻不容緩。
往昔研究海岸灘線的變遷，大多以現場調查海 岸地形的方式來研究，此方法不但耗時且探討範圍 有限。若欲分析灘線長期及季節的變遷，則需要長
期且每年多次的測量資料，才能獲得可靠的結果，
但是，在臺灣海岸工程規劃的經費限制及操作時程，
很難有充足的地測量資料來分析地形及灘線的變 化。相較於光學或聲學的現地測量，衛星影像具有 大範圍、便宜及省時之優點，且近年來，因衛星影 像的解析度一直提升，在國內外的學者也提出以衛 星影像來監測灘線的技術。因此，在外傘頂洲無充 足的地測量資料下，本文藉由 SPOT5、SPOT6 及 SPOT7 衛星影像處理技巧及海岸工程的知識基礎 下，探討外傘頂洲的長期及季節的灘線變化，即零 米線的變化。
以衛星影像的水線探討灘線變遷的往昔研究，
如 陳 良 健 及 徐 啟 崇 (1999) 、 Kurosawa and Tanaka(2001)、Gilvear and Gratiot (2004)、Ryu et al.
(2002)、王秀雯等(2007)、蕭國鑫等(2007)、楊勤儀 (2007)及張憲國等(2012)。但是，灘線的位置會受 到不同時間之潮汐或波浪而有變化，因此，Boak and Turner (2005)指出，若直接採用瞬時擷取出之 水線位置，來分析灘線的變遷是不可靠的。Gardel and Gratiot (2005)提出，若比較不同時間衛星影像
的水線，需考慮影像拍攝當時潮位及底床的前灘坡 度(beachface slope)。吳哲榮及吳啟南(2003)以潮位 接近之影像來分析臺灣西海岸之變遷。然而，若不 同影像拍攝時間之潮位差距過大或前灘坡度較緩，
評估的海岸變遷就會產生很大的誤差。張憲國及陳 蔚瑋(2005)使用灘線平移修正法(One-line Shifting Method, OSM)計算前灘坡度，並平移不同潮位的 水線至平均灘線，再以不同時間的平移後灘線探討 灘線變遷。吳政杰(2010)以實測坡度平移擷取水線 至平均灘線，探討臺南黃金海岸的灘線變遷。
平移衛星影像的水線至實際灘線的方法，需 要有潮位資料才能進行，但外傘頂洲無實測潮位站，
所以無法進行水線至灘線的修正。面對無潮位站的 問題，雖然可用衛星測高資料或海洋模式來推算潮 汐，如 Cartwright and Ray(1990)的正交潮法，
Matsumoto et al. (1995)的 NAO99b 模式。但是，在 臺灣西岸雲林嘉義一帶的淺水區潮位，陳蔚瑋及張 憲國(2008)發現 NAO99b 模式的推算值就與實測值 就有頗大的誤差，而提出潮位資料融合法 (data blending)來改善。張憲國及黃金維(2001)及顏沛華 等(2010)以臺灣海域十個潮位站的實測潮汐，評估 NAO99b 模式推估值的精度，在工程上是尚可接受 的。在外傘頂洲無潮位站的條件下，本文採用 NAO99b 模式初步推算外傘頂洲的潮位，再以陳蔚 瑋及張憲國(2008)的資料融合法修正 NAO99b 的初 估潮位。
往昔對外傘頂洲海岸變遷之研究，吳啟南等 (1992)利用外傘頂洲 1973 年至 1990 年間之衛星多 光譜影像，分析海灘變遷後發現，於高潮位時外傘 頂洲的沙洲面積急遽減少。蔡金吉及黃清和(1992) 由 1991 及 1992 年的實測數據分析外傘頂洲的地形 變遷後發現，在東北季風期間，北端海灘屬於堆積 斷面型態(normal beach)，海岸線往海側前進，但在 外海遠灘並無堆積成沙洲(sand bar)，而南端海岸在 碎波點附近產生沿岸沙洲，屬於漸變型(transient type)海岸。分析 1940~1990 年間地地形圖資，海岸 線平均往東後退 85 m。但分析 1992 年資料的一年 間，外傘頂洲會向東南方向內移，最大的後退約
110 m，平 均約為 70 m。林 宏仁等 (2003) 蒐集 1900~2002 年雲嘉海岸地形資料，分析外傘頂洲侵 淤量及海岸線變化後發現，2000 年以後沙洲往南 延伸，平均變遷率為每年 60 m。林宏仁等(2003) 提出造成沙洲地形變遷的主要兩個因素分別為，自 然力影響及人為開發。前者包括波潮流作用、颱風 及河川輸沙等自然力，後者包括超抽地下水、海埔 地開發、港灣及水庫興建等。張憲國及陳蔚瑋(2005) 分析 1993 到 2004 年間的外傘頂洲衛星影像發現，
近 11 年的外傘頂洲的陸地面積，以每年 0.35 km2 增加，且沙洲整體逐漸向臺灣本島靠近，有陸化的 趨勢。蔡雨農等(2015)以外傘頂洲往昔圖資、衛星 影像及 1993~2014 年間地形測量資料，分析外傘頂 洲灘線、陸上面積及體積變化，由結果顯示，沙洲 有東移現象，北部的移動速度小於南部的移動速度，
此沙洲變動特性說明，沙洲是以逆時針旋轉的，且 推估 2028 年以後沙洲將沒入平均海水面。本文在 往昔研究外傘頂洲變遷的基礎上，以衛星影像的圖 資來探討外傘頂洲灘線及陸上面積的長期及季節 性變化。
2. 影像處理及水線擷取
本章節主要介紹研究基地、圖資料來源、影像 前置處理及海陸交接水線之擷取方法。
2.1 研究基地
本研究選擇外傘頂洲為研究基地，其地理位置 於臺灣西海岸雲林縣口湖鄉海岸之外側，鄰近北港 溪出海口，示於 Google Earth 的圖 1。外傘頂洲為 臺灣最大的離岸沙洲，範圍從北港溪出海口往南延 伸至東石漁港一帶，呈東北至西南走向，海岸線長 度約為 13 km、寬度介於 1~3 km。
圖 1 研究基地之地理位置圖(引自 Google Earth)
2.2 影像來源
考慮衛星資料取得之方便性及空間解析度，使 用中央大學太空及遙測研究中心能提供之 SPOT-5、
SPOT-6 及 SPOT-7 的影像。上述三種衛星資料的 基本特性示如表 1：
表 1 衛星資料之基本特性
衛星代號 光譜模式 光譜範圍 (μm) 空間解析度
SPOT-5
多頻譜態
0.50~0.59(綠)
10 m 0.61~0.68(紅)
0.79~0.89(近紅外) 1.58~1.75 (短波紅外)
全色態 0.48~0.71 5 m
超解像模式 2.5 m
SPOT-6 多頻譜態
0.455~0.525(藍) 0.53~0.59(綠) 6 m 0.625~0.695(紅) 0.76~0.89(近紅外)
全色態 0.455~0.745 1.5 m
SPOT-7 多頻譜態
0.455~0.525(藍) 0.53~0.59(綠) 6 m 0.625~0.695(紅) 0.76~0.89(近紅外)
全色態 0.455~0.745 1.5 m
SPOT-5 號衛星於 2002 年 5 月發射升空，其多 光譜及全色態空間解析度分別為 10m 及 5m。因其 搭載多組感測器，可同時拍攝影像再經影像融合處 理後，空間解析度可提升至 2.5m，稱超解像模式 (Supermode)之全色態影像。SPOT-6 及 SPOT-7 分 別為 2012 年 9 月、 2014 年 6 月發射升空，兩者 的 多 光 譜 及 全 色 態 空 間 解 析 度 分 別 為 6m 及 1.5m。
選取的衛星影像的資訊，衛星影像地拍射時間 及衛星種類，如表 2 所示。因 SPOT-6 及 SPOT-7 為近年才發射升空的衛星，可選取的影像並不多，
故大多數的衛星影像為 SPOT-5 系列。三種衛星影 像之大地基準座標系統皆為 TWD97 座標系統。另 外，本研究亦收集以光學雷達(Light Detection and
Ranging, LiDAR) 航拍外傘頂洲影像的平均灘線，
當為檢驗衛星影像獲得灘線的基準，來評估本方法 的可行性。此航拍時間為 2013 年 7 月 29 日 11 點 08 分，航拍圖之空間解析度為 1 m，大地基準座標 系統亦為 TWD97。
表 2 外傘頂洲衛星影像圖資訊 日期(年/月/日) 時間(時:分:秒) 衛星種類
2007/11/04 02:22:25 SPOT-5 2008/06/30 02:23:19 SPOT-5 2009/03/17 02:23:23 SPOT-5 2009/10/15 02:51:21 SPOT-5 2010/02/23 02:33:43 SPOT-5 2010/08/14 02:27:16 SPOT-5 2011/07/27 02:52:45 SPOT-5 2011/09/18 02:32:41 SPOT-5 2012/08/21 02:24:21 SPOT-5 2013/02/18 02:37:03 SPOT-5 2014/02/18 01:58:42 SPOT-5 2014/12/30 02:19:41 SPOT-7 2015/08/03 02:05:53 SPOT-6 2015/12/28 02:26:33 SPOT-6
2.3 影像前置處理
由於每張衛星影像的拍攝角度、反射率等特 性不同，須先進行影像前置處理，方能獲得可用的 水線。影像前置處理之流程示如圖 2：
圖 2 影像前置處理的流程圖
因多頻譜衛星影像空間解析度較全色態影像 低，練智恒等(2010)、陳繼藩及陳敏新(2010)、鄭 義霖(2014)建議可用全色態影像的高解析度之空 間分布特性提升多頻譜態特性的空間解析度，如以 IHS 轉換法(Intensity-Hue-Saturation fusion)置換多 頻譜態特性的空間分布，再反轉換至原多頻譜態特 性。本研究使用往昔的影像處理技巧，計算多頻譜 的 反 射 率資 訊成 NDWI (Normalized Difference Water Index) 的圖像、提高海域及陸域兩者間灰階 值差異的影像強化，及以形態學來消除 NDWI 圖 像中的雜訊及填補圖像中的空洞等技巧。以下簡述 個技巧的原理及處理後的效果：
Step 1. IHS 轉換法
多頻譜影像比全色態影像有較多的空間資訊，
但空間解析度卻較低。因此，採用一般能提高空間 解析度的 IHS 影像融合法(Carper et al., 1990)，來 改善多頻譜影像的低解析度問題。IHS 融合的原理 為高解析度全色態影像的單層空間資訊融合於多 頻譜中一個影像空間資訊，再透過 IHS 反轉換後 的全色態影像，全色態影像就保有原彩色多層空間 資訊且具有高的空間解析度。
若紅光、綠光及藍光所組成的 RGB 影像轉換 至色相(Hue)、飽合度(Saturation)及亮度(Intensity) 的 IHS 色彩空間，色相及飽和度為色彩資訊，亮 度為空間資訊。以多頻譜影像的紅光、綠光及近紅 外光波段亮度值置換原本 RGB 影像中三個分量，
進行 IHS 融合。IHS 融合原理的詳細介紹可見 Gonzalez and Woods (2008)。
圖 3 為 IHS 融合前後外傘頂洲衛星影像。比 較圖 3 融合前後的左圖及右圖影像可知，右圖比左 圖能較清晰地顯示出地貌之輪廓。若再取 IHS 融 合前後，各波段亮度值的斷面分布來說明，IHS 融 合前後的亮度值差異，示如圖 4。由圖 4 中可知，
IHS 融合前後之各波段亮度值的分布特性相似。圖 3 及圖 4 的結果證實，融合後的多頻譜影像除保有 各波段之分布特性外，還可提高空間解析度。
(a) IHS 融合前 (b) IHS 融合後 圖 3 外傘頂洲衛星影像圖 IHS 融合前後之比較
(a) IHS 融合前
(b) IHS 融合後
圖 4 IHS 融合前後在 AB 斷面上各波段亮度值分 布的比較
Step 2.常態化差異水體指數
Lillesand and Kiefer(1994)提出在乾燥土地的 反射率最大為近紅外光，而最小為綠光，但是，在 清澈水體的反射率則相反。因此，Mcfeeters (1996) 定義常態化差異水體指數為近紅外光及綠光之差 和比，並用此指標來判斷植物表葉面之含水量多寡。
NDWI 的計算公式為：
Green Nir
Green NDWI Nir
... (1)
其中 Nir 為近紅外光波段反射率、Green 為綠 光波段反射率。近紅外光及綠光波段的反射率可由 影像亮度值轉換之(Gao, 1996)。計算公式為：
0 0
2
cos ) ) (
(
E d
L
...
(2a)其中
Bias ND Gain
L
... (2b)其中，ρ(λ)為無單位的行星反射率(Unitless planetary reflectance) ， Lλ 為 光 譜 輻 射 (Spectral radiance)，Gain 為增益值，ND 為影像亮度值，Bias 為 偏 差 ， d 為 日 地 距 離 (Earth-Sun distance in astronomical unit)，E0 (λ)為各種波段的太陽光譜輻 射量(solar exoatmospheric irradiance)，θ0為太陽天 頂角(solar zenith angle)。Gain、Bias 及 θ0列在衛星 影像資訊中。綠光、紅光及近紅外光的 E0分別為 1826、1554 及 1036(見韋玉春及黄家柱，2006)。
劉宇光等(2005)指出 d 值介於 0.983~1.017 之間，
若計算精度要求不高，可直接取 1。因此，採用劉 宇光等(2005)的建議，d=1。
圖 5 為轉換外傘頂洲的多頻譜影像成 NDWI 的灰階圖。由圖 5 可看出陸域的 NDWI 值約介於 0 到 1 間，而海域約介於 0 到-1 間。因此，NDWI 接近於 0 可當為界線來區別海域及陸域，但是，辨 識水線的位置是利用影像邊緣偵測的 Canny 法來 決定的。
圖 5 多頻譜影像融合後再轉換的 NDWI 圖
Step 3. 影像強化
當影像的灰階值接近時，影像明亮度的差異就 不大。倘若能修改明亮對比，可忽略非重點部分，
且加強重點部分，就可凸顯主題。影像強化原理為，
原影像強度值經一個映射函數轉換到新的輸出值。
影像強化效果與映射函數的設定有關，即輸入範圍、
輸出範圍及轉換函數曲率。影像強化乃影像處理的 基本技巧，可見陳繼藩及陳敏新(2010)。一般可採 用 gamma 轉換函數，且根據原影像明亮度的差異 性，適當調整 gamma 值約在 0.1 至 1.0 之間。本文 強化前後的外傘頂洲衛星影像比較示如圖 6。由圖 6 顯示，右圖影像強化後的海陸交界的明亮度差異 較強化前的左圖大，所以更可清楚顯示出外傘頂洲 的陸域範圍。
(a)影像強化前 (b)影像強化後 圖 6 衛星影像強化前後之比較
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
Step 4. 形態學
上述圖 6(b)的外傘頂洲南端外含有可能碎波 產生的白點雜訊，或者陸域內有些黑點雜訊，此雜 訊會造成在後續的水線辨識，及陸域面積計算上的 困擾。採用形態學(morphology)的技巧可消除不必 要的雜訊問題。形態學透過不同結構元素，進行膨 脹(dilation)或侵蝕(erosion)的組合運算，來分割影 像、濾除雜訊及填補孔洞。其中，斷開(opening) 為目標影像經侵蝕後再膨脹，具有平滑物體輪廓，
且消除較小之影像雜訊等效果。圖 7 為外傘頂洲的 斷開前後的比較。由圖 7 中可看到紅框內與外傘頂 洲無關之水花雜訊已被消除。
區 域 填 補 為 藉 由 膨 脹 及 適 當 的 補 集 (complement)或交集(intersection)等演算，填補影像 中的小洞。圖 8 為外傘頂洲區域填補前後的比較。
圖 8 中紅色框內顯示，原始影像中白色陸地所包圍 的黑色小洞，經區域填補後就完全消失，轉為白色，
而可被辨識為陸地。以處理後影像的白色像素面積 總合當為陸地面積時，區域填補可避免原本為陸地 而誤判成水域的問題，而能獲得正確的陸地面積。
(a) 斷開前 (b) 斷開後 圖 7 外傘頂洲斷開前後的比較
(a) 區域填補前 (b) 區域填補後 圖 8 區域填補前後實例圖
Step 5. Canny 的邊緣偵測法
海水的高或低潮位會在沙灘上造成不同位置 的乾濕交界的水線。因各波段在陸地及海水的反射 率不同，在水線處剛好是各波段的影像灰階差異較 大之處，所以，可利用此特性來擷取水線。本文選 擇邊緣偵測(edge detector)演算法中常用的 Canny 法(1986)，來擷取上述影像前置處理後的水線。
Canny 演 算法 包括梯 度值計 算，非 最大值 抑制 (non-maximum suppression)的細化，及適當決定門 檻(threshold)的二值化。本文採用 Fang et al. (2009) 的建議，以 Otsu 演算法來決定門檻。Canny 演算 法的偵測邊緣的能力已被證實較優於近似導數的 方法(Gonzalez and Woods, 2008)。
圖 9 為套疊前置處理後影像的擷取水線與原 衛星影像。由圖 9 可知，水線與衛星影像中的海陸 交界相當吻合。此結果說明使用的影像處理方法是 適用的擷取外傘頂洲的水線。
圖 9 擷取前置處理後影像的水線與原衛星影像的 套疊
3. 修正水線至灘線的平移法
在相同的高程基準下，比較不同時間的灘線位 置方能正確地分析出灘線的變化 (Chang et al., 2015)。本節說明修正所擷取影像的水線至灘線的 方法，包括方法所須潮位的資料融合、實測潮位、
平移方法及模式精度評估。
3.1 潮位推估與資料統計
外傘頂洲的平均坡度大約 1/100，而當地的潮 差約 3m，若影像選取時的潮位高度不同，辨識出 水線在高潮位或低潮位，造成水平差異可達 300m。
若要平移影像的水線至平均水位上，就必須先知道 影像拍攝當時的潮位。然而，外傘頂洲現地因無設 置潮位站，所以沒有實測潮位資料可以引用。若使 用如 Matsumoto et al. (2000)的 NAO99b 數值模式 來 推 算 此 處 的 潮 位 ， 外 傘 頂 洲 卻 在 淺 灘 區 ， NAO99b 模式所推算的潮位與實際值會有大的偏 差(見陳蔚瑋及張憲國，2008)。因此，再利用陳蔚 瑋及張憲國(2008)的潮位資料融合法，來精準推估 潮位。
潮位資料融合法乃以 NAO99b 模式初估鄰近 潮位站的潮位與實測值的差值，再空間內插至目標 點後，並反饋修正原推算目標點的潮位。NAO99b 全球潮汐模式為 Matsumoto et al. (2000)所發展，可 推算臺灣近岸海域任一位置的逐時潮位。為證實潮 位資料融合法在外傘頂洲的適用性，選擇外傘頂洲 週邊的箔仔寮(Bozihliao)、東石(Dongshih)、將軍 (Jiangjyng)、澎湖(Penghu)四個潮汐測站，來說明 塭港(Wengang)測站的推算精度。外傘頂洲及週邊 測站之地理位置，如圖 10 所示。
圖 10 外傘頂洲與週邊測站之地理位置圖
塭港測站之資料融合的計算公式為
ℎ = ℎ + ∆ℎ ℎ − ℎ , ℎ − ℎ , ℎ − ℎ , ℎ − ℎ
...
(3)其中，h 為潮位高，上標 n 代表 NAO99b 模式 所推算潮位，上標 o 代表潮位站的實測數據。下標 則為地點，塭港縮寫為 W、箔仔寮縮寫為 B、東石 縮寫為 D、將軍縮寫為 J、澎湖縮寫為 P。∆h 為二 維空間內插運算子。
NAO99b 模式推算原歷時潮位、實測值及資料 融合法的修正潮位，分別示如圖 11 的藍色虛線、
黑色點及紅色實線。由圖 11 可看出，模式推側值 與實測值有明顯的差異，但是，資料融合法後的修 正值幾乎吻合於實測值。
圖 11 模式推算塭港的潮位、資料融合修正值及實 測值的時序圖
一般以均方根誤差(Root Mean Squared Error, RMSE)評估模式的精度。均方根誤差的定義為：
RMSE = ∑ ( )
...
(4) 其中 hpi 及 hoi分別為時間 ti 的估計值及實測 值。RMSE 越大代表兩者偏離越多，估計值與實測 值越不吻合；反之，RMSE 越小代表兩者偏離越少，估計值與實測值越吻合。
為量化圖 11 中推側值與實測值之差異大小，
計算塭港 2007 至 2015 年的推算潮位與實測值之每 年的均方根誤差(表為 RMSEn)，及資料融合後的修 正值與實測值之 RMSEc，分別示如表 3 之第二及 第三行。由表 3 的第二行可知，在 2007 至 2015 年間， NAO99b 模式原始推算每年的潮位與實測 值的 RMSE 在 38.62 cm 至 41.34 cm 間，其平均值
119oE 20' 40' 120oE 20' 40' 121oE
23oN 12' 24' 36' 48' 24oN
箔子寮 澎湖馬公
東石 塭港
將軍 外傘頂洲
及標準差分別為 40.24 cm 及 1.04 cm。第三行的 RMSEc在 3.01 cm 至 8.52 cm 間，其平均值及標準 差分別為 5.44 cm 及 2.00 cm。比較各年的 RMSEn
與 RMSEc的結果，明顯降低 30 cm 以上。
為 瞭 解 降 低 百 分 比 ， 再 計 算 (RMSEn-RMSEc)/RMSEn的相對 RMSE，並列於表 3 的第四行。當此值接近於 100%，表示原始推算 潮位與修正潮位的差異很大，且修正潮位接近於實 測潮位。反之，若相對均方根誤差接近於 0%，表 示修正潮位接近於原始推算潮位，修正量不明顯。
由表 3 的第四行可知相對 RMSE 的可達 78.24%到 92.70%，此高相對 RMSE 表示資料融合法的修正 值可大幅改善原始模式推估值，且更加吻合於實測 值。
表 3 推算塭港測站自 2007 至 2015 年的每年潮位 與實測值之 RMSE
年份 RMSEn(cm) RMSEc(cm) 修正比(%) 2007 38.62 7.81 79.78 2008 39.20 6.87 82.47 2009 39.85 4.92 87.65 2010 40.85 3.76 90.80 2011 41.23 3.01 92.70 2012 39.16 8.52 78.24 2013 41.14 4.11 90.01 2014 40.81 3.66 91.03 2015 41.34 6.26 84.86
因塭港測站與外傘頂洲距離僅約 5 km，故資 料融合法可應用於外傘頂洲的潮位修正。外傘頂洲 的資料融合修正公式，如下式(5)所示
h = ℎ + ∆ℎ ℎ − ℎ , ℎ − ℎ , ℎ − ℎ , ℎ − ℎ
...
(5) 其中外傘頂洲潮位的地點下標縮寫為 Wa。3.2 灘線的平移法
為證實擷取影像的水線，再平移至平均灘線的
修正方法是可行的，以 2013 年現地的實測灘線為 基準，來評估所決定灘線位置的誤差。所收集地形 資料的測量日期為 2013 年 7 月 28 日至 31 日，為 避免比較不同時間的灘線已有位置上的差異，所以，
選擇與實測地形的時間須相近且的拍攝時間，即為 7 月 29 日 11 點 08 分的影像。
探討海岸的近岸水流(wave-induced current)及 輸砂等特性，多分解成沿岸及離岸方向兩個分量來 分析。因此，為方便後續探討灘線的變化，將原本 TWD97 的座標旋轉至以沿岸為橫軸且其法線為縱 軸的座標系統。新座標系統透過順時針旋轉 41.017 度及再平移(-1812804 , -1861112)來轉換原 TWD97 的座標。旋轉後外傘頂洲的水線如圖 12 所示。圖 12 上方曲線為外傘頂洲面向臺灣海峽的海側水線，
下方曲線為面向臺灣本島的陸側水線。海側及陸側 的水線的縱軸值差異較原 TWD97 座標小，而易於 後續分析及理解海側水線及陸側水線的變化特 性。
為了方便探討在某剖面的底床坡度或灘線在 不同時間的變化特性，每張影像所擷取的水線像素 點須要在相同橫軸值。因此，把所有類似圖 12 的 水線像素點內插至固定的橫軸上。橫軸自 3900 m 到 16700 m 間，每 100 m 的間隔進行線性內插，因 此，海側及陸側共有 129 個內插點。
圖 12 沿岸及離岸座標下外傘頂洲的水線
使用三角幾何的正切來平移水線至灘線，示如 圖 13 的示意圖。此方法須計算海岸法線的前灘坡 度，tanβ，再配合推算當時的潮位，△h=WL-MSL，
其中 WL 為影像拍攝時間的當時潮位，MSL (Mean Sea Level)為平均潮位。當 tanβ 及△h 已知，計算
△y=△h /tanβ 後，圖 13 中水線的離岸座標 y1為已 知，而未知的平均水位的灘線座標，可由 y0 =y1-△y 計算而獲得。
