光學式衛星影像雲層處理之研究

Volume17, No2, August 2013, pp. 115-134

1

國立台灣大學地理環境資源學系 博士 收到日期:民國 101 年 02 月 21 日

2

國立台灣大學地理環境資源學系 教授 修改日期:民國 101 年 05 月 02 日

＊

通訊作者, 電話:

02-33665830 ext.12, E-mail: [email protected]

接受日期:民國 102 年 01 月 21 日

光學式衛星影像雲層處理之研究

徐逸祥

¹

朱子豪

^2*

摘 要

利用光學式衛星影像進行土地利用判釋或農作物產量估測時，雲層覆蓋是無法避免的干擾之一。就 具有厚雲層的影像而言，本研究以單時期影像及區域增長 (region growing) 之方式偵測並切除無法還原地 物資訊的厚雲層及其雲影。具有薄雲的影像則以傅利葉 (Fourier) 分析建立薄雲的數學模式，再以此模型 薄雲並還原薄雲底下的地物光譜資訊，雖然在模式建立階段需兩時期影像，但建立後模式對其它影像進 行去雲處理時則僅需單時期資訊。研究成果顯示，厚雲及雲影偵測之整體精度皆可達到 90%以上。薄雲 去除方面，薄雲過濾器提升了約4%的分類精度，亦減輕薄雲對正規化差異植被指數 (normalized difference vegetation index, NDVI) 的影響，改善程度在統計上皆達到顯著性 (p < 0.01)。本研究成果可應用在土地 利用判釋和農作物產量估測中的影像前處理程序，除減少去除雲層的人力，亦可增加衛星影像的利用度。

關鍵詞：去雲、光學式衛星影像、地物判釋、區域增長、傅利葉分析

1. 前言

臺灣地處亞熱帶，又位於大陸與海洋的交界 處，為名符其實的多雲多雨區。根據林務局農林航 空測量所於2006 年所提供的資料顯示，全台各地 適合航照拍攝之天數，北部一年為50 天，中部一 年為 70 天，南部為一年 50 天，東部一年則僅有 20 天。雖然中央大學太空與遙測中心的衛星接收 站每年接收大量SPOT-2、-4 及-5 等衛星影像，國 家太空中心的接收站亦可每日接收福爾摩沙衛星 二號影像，但實際上可供使用的無雲影像卻很少。

以SPOT 影像的接收狀況來說，三顆衛星同時接收 至少要二個月才有完全乾淨無雲的影像；以福衛二 號影像的接收狀況來看，亦要三個月才有無雲影 像。而以美國太空總署 (National Aeronautics and Space Administration, NASA) 發射的 MODIS Terra 衛星影像來看，雖然每天皆能拍攝臺灣地區，但依 據2010 年的統計顯示平均雲覆率高達約 80%，全 年僅約10 天雲覆率在 10%以下。由此顯示了臺灣 地區雲覆量的嚴重性。

國內學術界及業界在以遙測進行土地利用的

調查時，最常使用且較易分析的影像仍為光學式影 像。以內政部為例，其最新的國土利用調查計畫，

整體計畫期程為2006 至 2015 年，並已於 2006 至 2008 年完成了全國土地利用資料的建置 (內政部 國土測繪中心，2008)。而自 2001 年開始的營建署 國土監測計畫，主要核心為利用衛星影像進行土地 變遷偵測，週期性的掌握國土利用變遷資訊，遏止 非法土地使用 (內政部營建署，2008)。以上所提 及的調查計畫，皆以遙測影像判釋為主體資料，其 餘再以地面調查進行補強，而未來各項國土利用調 查或監測更將朝向高頻率、高效率的機制邁進，對 遙測影像的需求亦將愈來愈高。

除了國土利用方面的應用外，以遙測影像估算 生物量 (biomass) 並應用在全球環境變遷以及農 作物產量估算，也是近年備受重視的課題。目前遙 測估算生物量主要是以衛星影像的可見光及近紅 外光 (near infrared, NIR) 波段計算各種植被指數 (vegetation indices)，再以迴歸等方式校估出植被指 數和淨初級生產量 (net primary productivity, NPP) 之間的關係，進而推估植物的生物量；了解植物的 生物量，即可推算碳儲存量，進而能了解碳循環的

機制以深入探討各種全球環境變遷的議題 (Dong

et al., 2003; Tan et al. 2007)。而在農作物應用上，

則可間接進行農作物產量的預估 (Doraiswamy et

al., 2004; Bsaibes et al., 2009; Wang et al., 2010)，因

此可協助政府掌握糧食的生產狀況及安全存量，避 免因糧食短缺而造成糧價波動或是饑荒發生。然而 在多雲多雨區可見光及近紅外光波段常會受雲層 的干擾，因而影響生物量估算的準確性，也使得相 關研究人員需花費大量成本進行去雲的前處理工 作。因此就上述的相關研究而言，光學式影像中的 高自動化及高效率雲層處理發展亦勢在必行。

過去研究中所發展的雲層偵測或去除方法，大 致 上 可 分 為 空 間 域 (spatial domain) 及 頻 率 域 (frequency domain) 兩大類。空間域的雲層處理是 指直接對影像平面上的像元灰階值 (gray-level) 進行計算，以得到雲層偵測或去除的結果；頻率域 則 是 利 用 如 離 散 傅 利 葉 轉 換 (discrete Fourier transform) 和 離 散 小 波 轉 換 (discrete wavelet transform) 將空間域的影像轉成頻率域 (frequency domain) 後再進行雲層處理。而頻率域處理的相關 演算法主要是用在薄雲的還原上。

目前已有許多研究證實雲層對於特定的紅外 波段有相當特殊的反射及吸收特性，並以此偵測雲 層 (Bantges et al., 1999; Mannozzi et al., 1999;

Ahmad & Hashim, 2002; Han, 2004)。然而許多文獻 是以經驗閾值來偵測和去除厚雲，而其使用的閾值 常常僅適用於特定的衛星影像或是特定的地區。隨 著 資 訊 工 程 發 展 出 自 動 閾 值 選 取 演 算 法 (automated threshold selection algorithm, ATSA)，厚 雲處理的研究便引用來偵測並去除厚雲 (Yang et

al., 2007)，也有研究以多變量統計為基礎的監督式

(supervised) 和非監督式 (unsupervised) 分類法亦 常用來進行厚雲的偵測 (Irish et al., 2006; 張立雨 等，2007；蔡博閎，2009)，或是以多時期及不同 來源的影像來進行偵測 (Hagolle et al., 2010)。而具 有短波近紅外光 (shortwave infrared, SWIR) 以及 熱紅外光 (thermal infrared) 的衛星影像，相較於 僅有可見光及近紅外光的衛星影像而言，其優勢在 於較能區別厚雲和雪、明亮土壤和建物等地物，如

Irish et al. (2006) 即是利用 Landsat-7 ETM+其中 5 個波段，結合非監督式分類法以及閾值設定，並利 用已知的雲層、雪、明亮的土壤、植被以及水的光 譜 特性 ，發展 了自 動雲層 覆蓋 評估 (automated cloud-cover assessment, ACCA) 的方法進行雲層的 偵測。然而僅有可見光及近紅外波段的影像來說，

則常需要多時期且無雲的參考影像，來進行厚雲的 偵測，才能有較佳的效果 (Hagolle et al., 2010)。

在厚雲處理的研究中，以影像鑲嵌或是影像復 原 (restoration) 等方式來填補厚雲挖除後影像空 缺的部分也是重要的議題，對於無法還原地面資訊 的厚雲層遮蔽區，可利用多時期影像以多重解析度 離 散 小 波 轉 換 方 法 (multiscale discrete wavelet transform method) (曾筱婷，2005)、相對性輻射校 正 (徐逸祥，2006；徐逸祥等，2006)、支援向量 機 (support vector machines, SVM) (Melgani, 2006)、直方圖匹配 (histogram matching) (葉精國，

2007)、帕森方程式 (poisson equation) (蔡博閎，

2009) 等方式，拼接成一張完全無雲的影像，以供 製圖相關應用。而許多影像復原的研究雖然不一定 和 雲 層 處 理 有 直 接 相 關 性 ， 但 其 以 影 像 融 合 (fusion) (Roy et al., 2008) 以及區域最小平方套合 (local least squares fitting) (Rakwatin et al., 2009) 等方式對受雜訊干擾或是缺失的影像進行重建的 概念，亦可用於回復因雲層干擾而損失的地物資訊 之議題上。

較薄的雲層由於未完全遮蔽地面資訊，許多研 究認為可將此類雲層去除並還原地物的資訊，因此 發展了ATCOR (Atmospheric Correction) (Richter, 1996; Richter, 1997) 、 優 化 薄 雲 霧 轉 換 (haze optimized transformation, HOT) (Zhang et al., 2002;

Moro and Halounova, 2007)、傅利葉轉換 (賴格英 等，2003; Feng et al., 2004; Chen et al., 2009; Zhong

et al., 2010) 及小波轉換等 (Du et al. ,2002; Wang et al., 2005)。上述傅利葉和小波轉換兩種方法是以

頻率域的概念，假設薄雲在空間域中具有慢變化的 特性，此即頻率域中的低頻部分，而雲下的地物則 是屬於高頻部分 (Du et al. ,2002; 賴格英等，2003;

Wang et al., 2005)，因此可利用傅利葉和小波轉換

將空間域影像轉為頻率域影像後，再將低頻的部分 進行濾除，即可達到薄雲去除和地物資訊還原的目 的。

總結以上文獻，雖然可知已有許多雲層處理的 研究，但仍有一些難以突破之瓶頸。其一，雖有許 多研究成功以中、紅外光段進行雲層之偵測，但如 SPOT、福衛二號及 Quickbird 等僅有可見光及近紅 外光段，因此對多光譜影像使用者而言去雲處理仍 有困難；其二，ATCOR 及 HOT 等方法所需的輔助 資料太多，在真實世界中，這些參考資訊可能難以 取得；其三，研究方法多需設定經驗閾值，如紅外 光段的雲層偵測、傅利葉轉換的濾波器 (filter)，其 閾值可能只適用於特定衛星影像或是特定地區；其 四，對於去雲結果的優劣，通常是以非量化的方式 來進行視覺化評估，欠缺客觀性，且去雲過程通常 也會破壞原本的地物資訊，因此去雲影像對土地利 用分類或是生物量估算等相關研究的使用者來說 用途可能有限。

因此，本研究著重於改進上述之瓶頸，以多數 資源衛星影像所具有的綠、紅及近紅外光段做為研 究對象，並依前人研究歸納厚薄雲的特性，以此知 識為基礎，再以傅利葉分析建立薄雲模型和濾波 器，確立厚薄雲去除流程，最後再發展量化的檢核 機制，如此除可驗證去雲流程的成效性，也可提供 使用者去雲影像的可靠性資訊。

2. 材料與方法

2.1 研究材料

依據前人研究所探討，影響雲層偵測或地物光 譜 值 還 原 處 理 的 因 子 主 要 可 歸 納 為 地 形 效 應 (Richter, 1997)、地物均質程度 (徐逸祥，2006)、

以及兩時期影像的地物變遷，其中地物變遷會影響 薄雲模型的建立，因此為去雲的優劣的決定性因 素。而在雲層處理後，所產生的誤差以及資訊扭 曲，則會傳遞到後端應用，影響後端使用者應用去 雲影像進行生物量估算、自動化分類或製圖展示的 效果，以上的雲層處理因子及因雲層處理而影響後

端應用的關係可以圖1 的核心架構表示。依據此架 構，本研究以地形效應及地物均質程度兩因子做為 討論的變項，決定以下四種不同的研究區進行不同 地表環境條件對去雲成果的影響探討，各區的特徵 描述如下：

(1) 起伏小之均質區：以沙質海岸區為代表；

(2) 起伏小之非均質區：以聚落和農地分佈之平原 地區為代表；

(3) 起伏大之均質區：以較高海拔的林區為代表；

(4) 起伏大之非均質區：以包含聚落和農地分佈之 中低海拔山坡地為代表。

而本研究即依據上述不同的地表特性進行影 像資料收集。目前所收集的影像主要分為兩個群 組，一為用來建立薄雲影像資料庫，以做為薄雲模 型和濾波器建立的訓練資料，另一個群組則用來進 行厚雲及薄雲去除的驗證。薄雲影像資料庫所需的 影像是在以上四種研究區中各挑選10 組代表性影 像，每一組分別包含薄雲覆蓋和無雲覆蓋的影像，

其中的16 張影像如表 1 所示。厚雲及薄雲去除驗 證影像方面，本研究則在以上四種研究區中各挑選 三張代表性影像做為一組，一張為厚雲覆蓋的影 像，一張為薄雲覆蓋的影像，一張為無雲且與薄雲 影像時間相近的影像，如表2 所示。以上每張影像 皆只有綠、紅及近紅外三個波段，而挑選時間相近 的影像的目的即為使地物變遷的因子成為控制的 變項。

2.2 研究方法

依據上述的研究對象，本研究的方法大致可分 為厚雲層之去除、薄雲層模型建立及地物光譜值還 原、量化檢核機制之建立等三大部分。厚雲層的部 分是取單張影像，以自動閾值選取及區域增長的方 式進行偵測並同時切除雲層。薄雲層的部分本研究 以傅利葉分析建立薄雲模型。關於檢核機制的建立 方面，對於厚雲層的偵測結果，本研究將尋求遙測 影像判釋的專家先劃定影像中厚雲層的範圍，以此 範圍做為地真資料，以檢核本研究之偵測結果；對 於薄雲層移除後的影像，則是以影像分類法和正規 化差異植被指數 (normalized difference vegetation

index, NDVI) 檢核還原地物光譜值的效果。所有流 程如圖2 所示。茲詳述如下：

2.2.1 影像前處理

由於影像的灰階值會受到影像取得時的太陽 入射強度影響，而使得同一地物在不同時間拍攝時 會有不同的反射強度，因此在進行雲層的偵測和去 除處理前，本研究先將影像由數位數值 (digital number, DN)轉 為 大 氣 層頂 (top of atmosphere, TOA) 的輻射強度值 (radiance)，再轉為 TOA 的反 射率 (reflectance) (Chander and Markham, 2003)。

公式如下：

  ₍₁₎

其中 L^λ為某光段之 TOA 光譜輻射強度，單位為 Wm^-2sr^-1μm^-1；Qcalmin 為最小量化校正像元值 (即 DN 值為 0)，與 LMIN^λ相對應；Qcalmax為最大量化

校正像元值 (以 Landsat-5 TM 影像而言即 DN 值為 255)，與 LMAX^λ相對應；LMIN^λ為衛星接收的最 小光譜輻射強度，單位為Wm^-2sr^-1μm^-1；LMAX^λ 為衛星接收的最大光譜輻射強度，單位為Wm^-2sr^-1 μm^-1；Qcal為量化校正像元值，其值為DN。而求 出L^λ後，即可利用太陽入射的相關參數，算出TOA 的反射率 (reflectance)，如下式所示：

(2) 其中ρP 為無單位的 TOA 反射率 (reflectance)；Π 為圓周率；d 為太陽和地球的距離，單位為天文單 位；ESUN^λ為TOA 的太陽入射強度 (irradiance)；

θ_s為太陽天頂角，單位為度。而為方便厚薄雲的 處理，本研究將ρP 換算為 0 到 255 的 8bit 影像。

地物均質度 地形效應

兩時期地物變遷 影響雲層處理之因子

厚雲層及其雲影 之偵測去除 薄雲層偵測 及影像還原

製圖展示 雲層處理程序

影像自動化分類 後端應用

影響

生物量估算

誤差或處理後產生之資訊扭曲效應

圖1 研究核心架構圖

雲層覆蓋 反射率影像

影像加強

、ATSA及區域增長

統計歸納 厚雲層

及雲影

薄雲去除 及地物光譜特性還原

薄雲層

薄雲去除影像

薄雲覆蓋影像 資料庫

無雲層影像 資料庫

傅利葉轉換

薄雲過濾器

去雲成效 之量化檢核結果

迴歸、影像分類、

NDVI檢核 厚雲及雲影去除影像

專家法檢核

雲層覆蓋影像

初步大氣及輻射校正

圖2 研究流程圖

表1 薄雲影像資料庫之影像列表 (部分)

影像編號 影像類型 地表類型 雲覆類型 厚/薄雲/無雲百分比 影像中心座標

經度 / 緯度 影像日期 1 Landsat-5 TM 起伏小之均質區 薄雲 0 / 40 / 60 -114.6 / 39.4 2005/04/17 2 Landsat-5 TM 起伏小之均質區 無雲 0 / 0 / 100 -114.6 / 39.4 2005/04/01 3 Landsat-5 TM 起伏小之均質區 薄雲 0 / 50 / 50 -115.3 / 39.7 2005/07/13 4 Landsat-5 TM 起伏小之均質區 無雲 0 / 0 / 100 -115.3 / 39.7 2005/07/29 5 Landsat-5 TM 起伏小之均質區 薄雲 10 / 50 / 40 -116.5 / 38.5 2005/07/13 6 Landsat-5 TM 起伏小之均質區 無雲 0 / 0 / 100 -115.3 / 37.3 2005/07/29 7 Landsat-5 TM 起伏小之均質區 薄雲 5 / 25 / 70 120.8 / 23.4 2005/01/02 8 Landsat-5 TM 起伏小之均質區 無雲 0 / 0 / 100 120.8 / 23.4 2005/01/18 9 Landsat-5 TM 起伏小之均質區 薄雲 10 / 70 / 20 121.0 / 24.5 2005/01/18 10 Landsat-5 TM 起伏小之均質區 無雲 0 / 0 / 100 121.0 / 24.5 2005/01/02 11 Landsat-5 TM 起伏小之非均質區 薄雲 0 / 20 / 80 -115.6 / 39.5 2005/01/18 12 Landsat-5 TM 起伏小之非均質區 無雲 0 / 0 / 100 -115.6 / 39.5 2005/01/02 13 Landsat-5 TM 起伏小之非均質區 薄雲 0 / 10 / 90 120.8 / 23.5 2008/10/25 14 Landsat-5 TM 起伏小之非均質區 無雲 0 / 0 / 100 120.8 / 23.5 2008/10/09 15 Landsat-5 TM 起伏小之非均質區 薄雲 15 / 30 / 55 -116.5 / 38.5 2008/10/09 16 Landsat-5 TM 起伏小之非均質區 無雲 0 / 0 / 100 -115.3 / 37.3 2008/10/25

表2 厚雲、雲影及薄雲去除驗證組之影像列表

2.2.2 厚雲層之去除

由於可見光和近紅外光段的光皆無法穿透厚 雲層，因此除了遮蔽地物的反射光之外，厚雲層 亦會反射大量可見光及近紅外光段，使厚雲層區 域形成強烈的白光。由於厚雲層具有此特殊的光 譜特性，以前後期影像做變遷偵測必可準確描繪 出厚雲層的區界，但基於多時期影像需較高的影 像成本，因此在雲層的偵測方面，本研究仍以單 時段影像，結合三種方式進行處理，一為影像加 強 (image enhancement) ， 二 為 自 動 閾 值 選 取 (ATSA)，三為區域增長 (region growing)。此三種 演算法的結合可改進前人研究在厚雲處理時對較

薄的厚雲易發生漏判 (omission error)，而對都市 和裸露地等高反射地區易發生誤判 (commission error) 的問題，使偵測結果更為精準。

2.2.2.1 影像加強

利用厚雲層的高反射特性，可將影像中可見 光及近紅外光段皆高的區域界定為厚雲層分佈的 區域。在歸納出厚雲層的特性後，本研究可為影 像的各個光段定義出一閾值，將高於此閾值之部 分切除，即可對厚雲層的本體做初步描述。但由 於大氣狀況及太陽入射角度等因素影響，造成即 使為同一種地物，在不同時間感應器所感應之光 譜強度有所不同的結果，因此較不易以影像原始

影像編號 影像類型 地表類型 雲覆類型 厚/薄雲/無雲百分比 影像中心座標

經度 / 緯度 影像日期 1 Landsat-5 TM 起伏小之均質區 厚雲及雲影 60 / 0 / 40 -115.6 / 39.5 2006/05/06 2 Landsat-5 TM 起伏小之非均質區 厚雲及雲影 30 / 0 / 70 -116.5 / 38.5 2005/07/13 3 Landsat-5 TM 起伏小之非均質區 厚雲及雲影 0 / 0 / 100 -115.3 / 37.3 2005/07/29 4 Landsat-5 TM 起伏小之非均質區 薄雲 0 / 40 / 100 -100.3 / 38.1 2005/04/24 5 Landsat-5 TM 起伏小之非均質區 無雲 0 / 0 / 100 -100.3 / 38.1 2005/05/10 6 Landsat-5 TM 起伏大之均質及非均質區 厚雲及雲影 50 / 5 / 45 121.5 / 24.9 2002/03/09

的 DN 值來定義出固定之閾值。有鑑於此，本研 究在訂定閾值前，先將欲處理之影像進行標準差 延 伸 加 強 (standard deviation stretch enhancement)。其概念為先將影像的原始光譜反射 值分佈直方圖繪出，並統計其平均數及標準差，

並依據標準差將極端值消去，通常取平均值正負 二倍標準差之值做為消去及保留之閾值，被保留 之光譜反射值則延伸為0 到 255，成為一個新的分 佈直方圖。

2.2.2.2 自動閾值選取

Otsu (1979) 所發展的 ATSA 為目前最常使用 的閾值自動化選取的演算法，其假定每張影像的 直方圖是由物件和背景的直方圖所組成，以本研 究而言，物件即為厚雲，背景即為無雲區，厚雲 的直方圖會分佈在值較高的部分，而背景則是較 低的部分，兩者的直方圖會形成一個雙峰的直方 圖，因此若能找到一個閾值在誤差最小的情形下 在雙峰之間做切割，則就能將物件和背景做切 割。首先令 hi 為影像中每個反射率值的累計次 數，

i





1,...,

n 

，n 為影像的最大反射率值 (以本 研究而言將最大反射率換算為 255)。因此可計算 每個反射率值發生的機率如下：

  (3)

其中

p

i>=0 且 pi值的加總為1。若在 i 之間找到其

中一個值

k 做為暫定的閾值，由此可計算直方圖

第

m 階的累積動差分別為：

(4)

(5) ATSA 的目標即是要選擇最好的閾值 t，

 n 

t

1,..., ，使厚雲和背景能有最佳的切割。而 Otsu (1979) 的方法是以統計中的判別分析為基 礎，目的是要找到能將群組間的變異數_B²達到最 大化的閾值。_B²的公式如下：

  (6)

而當下式成立時，即可找出最佳的閾值

t：

(7)

2.2.2.3 區域增長

區域增長 (region growing) 先於某一個像元 為種子 (seed)，再選擇增長之方式為四方位或八 方 位 法 ， 再 給 定 一 個 閾 值 ， 以 光 譜 歐 氏 距 離 (spectral Euclidean distance) 來計算鄰近像元之光 譜歐氏距離是否在閾值之內，若是，則此鄰近之 像元就為增長之一部分，若否，則此鄰近之像元 就不為增長之一部分。光譜歐氏距離以下式計算。

(8) 其中

D 為光譜歐氏距離；n 為波段數 (number of

bands)；i 為某一波段；di為

i 波段中 d 像元的值；

e

i為

i 波段中 e 像元的值。透過自動閾值選取所產

生的閾值過濾，可偵測大部分厚雲本體，而厚雲 周圍較薄之厚雲則可以區域增長進行偵測，得到 更準確之結果，最後將偵測結果的範圍皆去除 後，即完成厚雲去除之程序。

2.2.3 厚雲層雲影之偵測

本研究雲影之偵測是修改Le Hégarat-Mascle and André (2009) 的方法，先建立考慮的太陽、衛 星、厚雲層以及雲影之幾何關係，再以加入影像 加強、Otsu (1979) 所發展的 ATSA 以及區域增長 等與厚雲偵測相同的程序，以改善雲影周圍較亮 但卻仍會影響地物光譜特性的區域，雲影偵測之 流程如圖 4 所示。而太陽、衛星、厚雲層以及雲 影之幾何關係公式如下所示。

(9) 其中為厚雲與雲影之關係之方位角 (以北方為 0 度，逆時針方向為正)； 為太陽方位角 (sun _ azimuth)； 為太陽天頂角 (sun zenith angle)；_ _

為衛星軌道與赤道之夾角； 為衛星感應器視角 _ (looking angle)。

資料來源：Le Hégarat-Mascle and André, 2009

圖3 太陽、衛星、厚雲層以及雲影之幾何關係示 意圖

厚雲偵測結果

依雲和陰影之幾何關係平移厚 雲偵測結果

依每一平移後之厚雲區塊範圍 進行ATSA找出陰影主體範圍

依ATSA結果進行區域增長

原始衛星影像

標準差延伸加強

厚雲雲影偵測結果

圖4 厚雲層雲影偵測流程圖

2.2.4 薄雲模型建立及地物光譜值 還原

模型的建構是以時間相近並為同一衛星感測 器以同一拍攝角度所攝得的無雲及薄雲兩時期影 像做為素材，並將兩時期影像進行傅利葉轉換，

再以變遷分析歸納薄雲在傅利葉頻譜中的特性。

由於兩時期影像時間相近，所以地物可視為無變 遷；使用同一衛星感測器攝得的影像，可避免因 感測器特性不同而造成同一地物具有不同的光譜 特性；使用同一拍攝角度的影像，則可避免地物 反射之輻射因經過大氣層路徑長不同而造成同一 地物具有不同的光譜特性。在以上所有影響影像 變遷的變因後皆控制後，唯一使影像變遷之原因 僅剩薄雲覆蓋，即可進行薄雲模型的建立。

在前人研究中，雖已知薄雲在空間域中具有

慢變化的特性，並發展以高通濾波器進行薄雲的 去除 (賴格英等，2003)。因此可假設所有資訊的 反射分量都存在反射率的變化，並且是頻域中的 高頻成分，而遙測中照射分量在整幅影像上除個 別陰影區外一般差異較小，表現出慢變化的特 徵，與低頻相聯繫，如下式所示：

其中，f (x, y) 為感測器接收到的綜合資訊，r

(x, y) 為地物或地面反射率，代表信號，t (x, y) 為

雲層的透射率，代表雜訊，L 為太陽輻射強度，a 為太陽輻射在大氣傳輸過程中的衰減係數，並且

r (x, y)、t (x, y) 和 a 的數值位於 0 與 1 之間。上式

可以通過限定條件，簡化成下式：

(11) 由公式(10)及(11)可知，只要利用高通濾波器 即可保留地物資訊，將薄雲的資訊去除。但其濾 波器並非經由實際的薄雲特性分析推導產生，因 此常造成非薄雲覆蓋區的光譜特性被改變甚至破 壞，如此可能降低去雲影像的實用性。由於傅利 葉分析對處理週期性雜訊的能力較強，因此本研 究在測試過程中，依據經驗暫定將無雲及薄雲兩 時期影像對分別分割為25×25像元大小的影像單 元，使得每個影像分割內的薄雲厚薄程度接近一 致，若將薄雲視為雜訊，則厚薄愈一致，代表此 愈接近週期性雜訊。將每一影像分割轉為傅利葉 頻譜圖後，再以變遷分析之概念將兩時期頻譜圖 相除，即可觀察薄雲在傅利葉頻譜圖中所呈現的 特性。由圖5可知，薄雲影響下兩時期傅利葉頻譜 相除後的影像會呈現一字形或是十字形的特性，

因此可以歸納此特性，得到薄雲之模型，以此模 型可產生濾波器，即可對所有的薄雲覆蓋影像進 行薄雲濾除，得到去雲影像。本研究以上述方式 將表1所有影像對所產生的25×25像元大小之頻譜 圖相除結果平均後，即可得到薄雲在傅利葉頻譜 圖中之模型。

在建立薄雲模型前，先透過傅利葉轉換將薄 雲和無雲影像的每個波段，分別由空間域轉換至

(10)

頻率域。而當薄雲和無雲影像都轉換至頻率域 後，再轉換為傅利葉頻譜 (spectrum) 圖。將薄雲 和無雲各個相同的波段之頻譜圖相除後，即可得 到衛星影像每一波段受薄雲影響後頻譜圖被增抑 的部分。由目前測試結果亦發現，薄雲不一定僅 影響頻譜圖中的低頻部分 (即接近中央部分)，在 高頻部分亦會受到增揚或是降低，由此也證明前 人研究中單純以高通濾波器進行處理會造成地物 資訊被破壞的原因 (賴格英等，2003; Feng et al., 2004; Chen et al., 2009; Zhong et al., 2010)。而將薄 雲影像資料庫裡所有的薄雲和無雲影像對都經過 此程序處理後，即可得到薄雲的模型。而依照此 模型，即可設計一個濾波器降低被增揚部分的頻 率，而提高被壓抑部分的頻率。其濾波器設計概 念如下：

(12)

其中

m

f為薄雲濾波器；Afi為第

i 個 25×25 像

元大小的單元中；無雲樣區的傅利葉頻譜強度；

A

fi‘為第

i 個單元中薄雲樣區的傅利葉頻譜強度；n

即為所有用來設計濾波器的單元數目。再來，即 可用此濾波器，對所有表 2 的薄雲驗證影像進行 處理。而影像經過濾波器處理後，即可再由頻率 域轉回為空間域的影像。由以上傅利葉轉換工具 所歸納分析出來之薄雲濾波器在應用於驗證影像 進行去雲處理時僅需單時期的影像即可進行。本 研究也將對所有影像的去雲成效進行檢核，以評 估薄雲濾波器進行去雲處理後，應用在土地利用 判釋以及生物量推估上的可行性。

2.2.5 量化檢核機制之建立

以本研究目前所參考之前人文獻中，對於雲 霧之偵測結果，皆僅以視覺化之方式對偵測結果 之優劣進行文字敘述性之評估，而未提出一個量 化的指標來檢核偵測之精度。有鑑於此，本研究 提出以專家法及影像分類法進行偵測精度評估，

以使評估具較高之客觀性。

2.2.5.1 專家法

專家法的目的是檢核厚雲層及雲影之去除成 效，在去雲研究中常用來做為主要的檢核工具 (Stowe, 1984; Berendes et al., 2004)，也有研究結合 人 工 判 釋 及 監 督 性 分 類 法 進 行 人 機 互 動 (man-machine interaction) 判釋，以此判釋的結果 做為厚雲的地真資料 (Yang et al., 2007)。在本研 究中之專家法，則是將衛星影像縮放至固定的比 例尺後，先以經驗將肉眼可見之厚雲霧區數化，

之後再尋求三位專家，檢核編修本研究之數化結 果。本研究將請三位專家在二倍標準差延伸加強 及無延伸加強兩種展示模式下協助確認厚雲霧之 判釋及數化，最後並取三者圈選結果之交集做為 本研究之地真資料。

2.2.5.2 影像分類法

在本研究中，影像分類法用於檢核薄雲霧去 除之成果以評估雲下地物光譜資訊的還原程度，

以及無雲區的地物光譜資訊被破壞的程度。首先 將原始薄雲覆蓋影像以及去雲影像分類進行影像 分類，再以土地利用資料進行分類結果之檢核，

檢核時會分別分析雲覆區以及無雲區的分類精 度。此方法之優點除能提供量化評估數據外，也 可探討各種工具所產生的薄雲模型是否會嚴重破 壞原本的地物資訊，以及各種地物資訊被影響的 程度差異。如此即可了解去雲影像應用於自動化 分類之可行性，提供土地利用研究使用者更多資 訊。

本研究所使用的分類演算法為非監督式分類 中 的 反 覆 自 我 組 織 資 料 分 析 法 (iterative self-organizing data analysis techniques algorithm, ISODATA)。首先將影像進行非監督分類為100類 後，再以人工判釋各類應歸屬的類別並同時進行 併類，最後即可得到分類之成果 (Peterson et al., 2009)。此方法之優點在於可先純化具複雜光譜特 性的地物類別，並能保有人工判釋之高準確性。

另一方面，McNemar 檢定則用以檢驗薄雲霧 處理前後的分類精度差異是否具有統計上的顯著 性。McNemar 檢定是以標準化常態檢定統計為基

礎，常用在不同影像分類演算法精度比較的相關 文獻中 (Foody, 2004; Ratle et al., 2010)，若|z| >

1.96，則代表兩種分類演算法精度上的差異具有統 計上的顯著性，而 z 值的正負號則代表兩種演算 法哪一種的表現較佳。z 值可以由公式(13)計算而 得：

(13)

在本研究中，

f

GP 代表薄雲處理後影像分類正確而 處理前影像分類錯誤的像元個數；fPG 則代表薄雲 處理前影像分類正確而處理後影像分類錯誤的像 元個數。

2.2.5.3 NDVI 法

常態化差異植被指數 (normalized difference vegetation index, NDVI) 一直是遙測間接估算生 物量的重要指標之一 (Dong et al., 2003; Tan et al., 2007)，雖然 NDVI 以比率法的方式進行計算，已 可減低大氣的干擾，但干擾還是無法完全消除 (Kaufman and Tanre, 1992)，尤其是當薄雲層愈厚 時，NDVI 愈難反應真實的地面狀況，估算的生物 量愈不準確。因此若能將薄雲覆蓋地區還原成無 雲時的真實地面狀況，則就能解決此一方面研究

的難題。為驗證本研究解決此一難題的效果，因 此將無雲、薄雲和薄雲去除後的影像轉換為NDVI 影像，並分別與無雲NDVI 影像相減，得到 NDVI 差值影像，以比較薄雲去除後是否確實能改善薄 雲干擾的狀況。

3. 結果與討論

非均質區包含許多土地利用類型，大部分土 地利用研究者的研究區皆為非均質，因此在本文 中選擇幾個較重要的案例影像進行討論，並會較 著重於呈現及討論非均質研究區之雲層處理成 果。

3.1 厚雲層及雲影偵測

本研究仍以單時段影像，結合三種方式進行 厚雲層的處理，一為影像加強，二為自動閾值選 取演算法，三為區域增長。此三種演算法的結合 可改進前人研究在厚雲處理時對較薄的厚雲易發 生漏判，而對都市和裸露地等高反射地區易發生 誤判的問題，使偵測結果更為精準。以下探討厚 雲偵測在各種研究區下之成效。

兩時期影像經傅利葉轉換 後之頻譜強度商值影像 原始無雲影像

建立薄雲過濾器 歸納薄雲霧在傅利葉 頻譜圖中影響的區域

原始薄雲霧覆蓋影像

圖5 由傅利葉轉換並分析產生薄雲過濾器之流程示意圖

3.1.1 起伏小之均質區

本研究以美國紐約州平坦的長島其及附近海 域之Landsat-5 TM 影像做為起伏小之均質區的案 例，影像日期為2006 年 5 月 6 日，此為以厚雲為 主體，周圍略帶薄雲分布之案例。以專家數化之 厚雲霧範圍做為地真資料來檢核厚雲霧處理的結 果，整體精度為95.10%，kappa 值為 0.743 (圖 6)。

由以上結果可知，雖然整體精度和 kappa 值所呈 現的結果有差異，主要是因為海面上厚雲漏判比 率高所造成的 (約達 29%)，雖然漏判的像元數不 多，對整體精度影響小，但以 kappa 值來評估正 確性時即可呈現出厚雲漏判的問題。陸地上的厚 雲偵測效果較海面上的為佳，可能原因為海面過 於均質，因此專家在肉眼判釋雲層時較難判定雲 本身是否仍帶有部分海面的資訊，而將色調較白 的區域直接判釋為厚雲覆蓋區，然而本研究之演 算法仍將其判別為非厚雲區，因此造成漏判的發 生。另一方面，部分反射率較高的都市地區以及 裸露地，則因為在綠、紅及近紅外波段皆有較高 的反射率，因此成為誤判發生的區域，但整體來 說誤判之比率較低 (約為 16%)。

圖6 起伏小之均質區厚雲霧處理之成果。(A)厚雲 霧偵測去除之結果 (黑色區域)；(B)專家劃定 之厚雲霧分佈圖 (黃色區域)。

3.1.2 起伏小之非均質區

本研究以美國內布拉斯加州農業及聚落分佈 區之Landsat-5 TM影像做為起伏小之非均質區的 案例，以2005年7月13日之影像做為案例一、以 2005年7月29日之影像做為案例二素材，分別代表 零碎厚雲以及大片厚雲之案例。以專家數化之厚 雲霧範圍做為地真資料來檢核厚雲霧處理的結

果，案例一影像之整體精度可達99.58%，kappa值 為0.945，案例二則為94.75%以及0.883 (圖7)

由以上結果可知，零碎厚雲的成果較佳，而 大片厚雲處理成果較差是因為大片厚雲中間及周 邊區域摻雜部分較薄的雲層，這些較薄的雲層帶 有部分地面資訊，卻已無法經由肉眼判別出地物 的種類，因此專家判別其為厚雲覆蓋區，但本研 究之演算法仍將其判別為非厚雲區，因此造成漏 判的發生。另一方面，部分反射率較高的都市地 區以及裸露地，則因為在綠、紅及近紅外波段皆 有較高的反射率，因此成為誤判發生的區域。

而厚雲雲影的偵測結果則如圖 8 所示。以專 家數化之厚雲雲影範圍做為地真資料來檢核厚雲 雲 影 處 理 的 結 果 ， 案 例 一 影 像 之 整 體 精 度 為 80.75%，kappa 為 0.703，案例二則為 83.73%以及 0.753。厚雲雲影的本體大致上皆能準確偵測，但 與厚雲偵測結果類似的是，在雲影邊緣部分的效 果較差。另外由於雲影的偵測是以厚雲的偵測結 果為基礎，因此其誤判及漏判會繼承下來而影響 雲影的偵測成效。

3.1.3 起伏大之均質及非均質混合區

本研究以臺灣新北市烏來山區之Landsat-5 TM影像做為起伏大之均質及非均質混合區的案 例，影像日期為2002年3月9日，在影像中同時有 零碎以及較大片之厚雲，亦有少部分的薄雲存 在。以專家數化之厚雲霧範圍做為地真資料來檢 核厚雲霧處理的結果，整體精度可達94.23%，

kappa值為0.802 (圖9)。與起伏小之非均質區中案 例二偵測結果類似的是，由於厚雲中間及周邊區 域摻雜部分較薄的雲層，這些較薄的雲層帶有部 分地面資訊，卻已無法經由肉眼判別出地物的種 類，因此專家判別其為厚雲覆蓋區，但本研究之 演算法仍將其判別為非厚雲區，故造成漏判的發 生，也是此研究區主要的判釋誤差來源。在誤判 方面，此山區由於人為開發較少，僅少數零星的 高反射率都市地區以及裸露地，因此此區的誤判 發生率較其它研究區明顯降低。綜合不同特性研 究區的結果來看，影響本研究厚雲處理優劣的最

(B) (A)

大因素是在地物組成特性，而地形起伏的影響則 相對較低。

厚雲雲影的偵測結果則如圖 9(C)所示。以專 家數化之厚雲雲影範圍做為地真資料來檢核厚雲 雲影處理的結果，整體精度為 70.73%，kappa 為 0.503。與 3.1.2 節案例相似的是，厚雲雲影的本體 大致上皆能準確偵測，但與厚雲偵測結果類似的 是，在雲影邊緣部分的效果較差；然而在地形起 伏造成的陰影區較明顯的地區 (如圖 9 各子圖的 左方中央部分)，則會發生較嚴重的誤判情況，也 因此造成結果不如起伏小的研究區理想。雖然本 研究的陰影偵測演算法在區別雲影及地形造成的 陰影成效有待加強，但若使用者的需求是要排除 所有陰影影響的地區，則本研究發展之技術仍有 助於減少人工篩選資料的人力成本。

3.2 薄雲偵測及去除成果

依據2.2.4 節所述，將薄雲影像資料庫中無雲 及薄雲兩時期影像對分別分割為 25×25 像元大小 的影像單元，並將每一影像分割轉為傅利葉頻譜 圖後，再以變遷分析之概念將兩時期頻譜圖相 除，即可觀察薄雲在傅利葉頻譜圖中所呈現的特 性。圖10 即為表 1 中薄雲影像資料庫之第一對影 像綠波段部分之成果，由結果可知，薄雲影響下 兩時期傅利葉頻譜相除後的影像會呈現一字形或 是十字形的特性。將表1 所有影像對所產生的 25

×25 像元大小之頻譜圖相除結果平均後，即產生薄 雲在傅利葉頻譜圖中之模型。最後經由公式(12)，

即可得到如圖11 至圖 13 所示的薄雲過濾器，並 可對表 2 中單時期之薄雲去除驗證影像進行薄雲 濾除，得到去雲影像。

(A) (B)

(C) (D)

圖 7 起伏小之非均質區厚雲霧處理之成果。(A) 案例一影像厚雲霧偵測去除之結果；(B)案 例一影像專家劃定之厚雲霧分佈圖；(C)案 例二影像厚雲霧偵測去除之結果；(D)案例 二影像專家劃定之厚雲霧分佈圖。

(A)

(C) (D)

(B)

圖8 起伏小之非均質區厚雲雲影偵測之成果。(A) 案例一原始影像；(B)案例一影像厚雲雲影偵 測結果 (綠色部分)；(C)案例二原始影像；(D) 案例二影像厚雲雲影偵測結果 (綠色部分)。

(A) 

(B) (C)

圖9 起伏大之均質及非均質混合區厚雲霧處理之 成果。(A)原始影像；(B)厚雲霧偵測去除之 結果；(C)厚雲雲影偵測結果 (綠色部分)。

(A)  (B) 

(C)  (D) 

圖 10 無雲及薄雲覆蓋影像綠波段部分之傅利葉 頻譜相除成果。(A)無雲影像；(B)薄雲影 像；(C)傅利葉頻譜相除成果；(D)(C)圖中 紅框部分之放大圖。其餘波段亦呈現相似 的特性。

圖11 綠波段之薄雲過濾器

圖12 紅波段之薄雲過濾器

圖13 近紅外波段之薄雲過濾器

本研究的薄雲去除驗證影像是以美國堪薩斯 州農業及聚落分佈區2005年4月24日之Landsat-5 TM影像為案例，比較前人發展之同態濾波 (賴格 英等，2003; Feng et al., 2004; Chen et al., 2009;

Zhong et al., 2010) 以及本研究發展之薄雲過濾器 的成果，並以時間相近的2005年5月10日Landsat-5 TM無雲影像比較兩種方法在去雲後是否能讓反 射率值接近無雲狀態。本研究在薄雲干擾區中取 12處共520個像元的屋頂、常綠植物及針葉林等區 域做為檢核點，繪製XY散佈圖，若趨勢線的斜率 愈接近1而截距愈接近0，以及迴歸的R²愈接近1，

則代表薄雲或去雲影像與無雲影像的光譜值愈接 近一致。

同態濾波之成果如圖14所示，由圖可知各波 段在斜率、截距以及R²方面的表現不一，去雲後 綠及近紅外波段之斜率較去雲前接近1，然而近紅 外波段的R²較去雲前低，且截距較去雲前大許 多。由整體來看，以同態濾波去雲的結果有低估 反射率值的現象。而本研究薄雲過濾器之成果如 圖所示，由圖可知薄雲去除後各波段趨勢線的斜

率都比去雲前的影像接近1，而截距亦皆較去雲前 的影像小，R²也都比去雲前影像大，然而整體來 看仍有低估的現象。與同態濾波的結果相比，本 研究發展去雲結果各個波段皆較去雲前得到改 善，雖然同態濾波在綠波段的改善程度較本研究 佳，但其它波段的表現則較不穩定，且低估的情 形較本研究之成果嚴重。因此可推論本研究薄雲 去除後的影像光譜值與完全無雲時的狀況仍較為 接近，表現較同態濾波佳。

由XY 散佈 圖雖 可 知本 研究 的 去雲 成果 較 佳，但不一定能直接證實本研究的去雲成果仍可 適用於影像分類以及生物量推估上。因此以下以 影像分類以及NDVI兩種方式來進行適用性的評 估。由圖16(A)和(C)可知，去雲後之影像在視覺上 薄雲霧影響確實有減少，可看清楚較多雲下的地 物，去雲處理前後影像分類之結果則如圖16(B)和 (D)所示。為了解薄雲過濾器對無雲區、雲影區之 影響，本研究請專家以人工方式將影像劃分為「雲 區」、「影區」及「無雲無影區」等三個區域，分 別計算薄雲處理前後的分類精度差異 (圖17)。由 表3可知，就全幅影像來探討，薄雲過濾器確實提 升了分類精度，並且McNemar檢定的z統計量也顯 示改善的程度具有統計上的高度顯著性 (z >>

1.96)。而就各分區來探討，過濾器對雲區的分類

精度提升最多，無雲無影區亦有少許提升，影區 的分類精度則反而下降，然而以McNemar的z統計 量來看，雖然薄雲過濾器顯著地降低影區的分類 精度，不過雲區和無雲無影區的精度提升則更具 顯著性。

而NDVI方面。由圖18可知，無雲與去雲處理 後NDVI影像相減後大於0.1的區域，較無雲與薄雲 覆蓋NDVI影像相減後大於0.1的區域顯著減少，

NDVI差值大於0.1的比例由去雲處理前佔全區 31.23%降為處理後的22.97%。再以方均根誤差 (root mean square error, RMSE) 來比較去雲前後 的NDVI差異，原始薄雲覆蓋影像與無雲影像 NDVI的RMSE值為0.21，而去除薄雲影像與無雲 影像NDVI的RMSE值則為0.15，兩者差異具有統 計顯著性 (p < 0.01)。

由上述結果可知，McNemar 檢定提供了除整 體精度和 kappa 值之外，另一種客觀的分類結果 差異解釋，證明了薄雲過濾器雖無法全面提升影 像各區之分類精度，然而其去雲的功效已有發 揮，對於 NDVI 分析而言，薄雲過濾器在統計上 也能有顯著的改善。不過值得注意的是，上述結 果也說明了除了雲影區之外，薄雲過濾器也可能 降低地形效應陰影區的分類精度或是 NDVI 的準 確性，此也是後續研究待解決的重點之一。

表 3 薄雲霧處理前後各區分類精度表

分區 薄雲霧處理前 薄雲過濾器處理後 z 統計量

全幅 77.74%*

0.297

81.69%

0.339

92.475

雲區 73.00%

0.183

79.18%

0.221

88.282

影區 79.52%

0.534

79.17%

0.521

-2.399

無雲無影區 81.74%

0.331

84.53%

0.365

49.673

*百分比代表整體精度，0~1 之間的小數代表 kappa 值。

(A) 

y = 0.4177x + 13.96  R² = 0.1277   

(C) 

y = 0.735x ‐ 0.172  R² = 0.2716   

(E) 

y = 0.4999x + 27.477  R² = 0.1914   

y = 0.6672x + 42.093  R² = 0.1289   

(F)  (D) 

y = 1.4128x ‐ 0.1412  R² = 0.3212   

(B) 

y = 0.8718x + 12.24  R² = 0.1630   

圖14 以同態濾波去雲前後各波段與無雲影像之反射率 XY 散佈圖。左排三張圖代表的 X 軸代表薄雲影 像的反射率值，Y 軸代表無雲影像的反射率值；右排三張圖代表的 X 軸代表去除薄雲後的影像反 射率值，Y 軸代表無雲影像的反射率值。

(A) 

y = 0.4177x + 13.96  R² = 0.1277   

(C) 

y = 0.735x ‐ 0.172  R² = 0.2716   

(D) 

y = 0.8552x ‐ 1.131  R² = 0.4803   

(E) 

y = 0.4999x + 27.477  R² = 0.1914   

(F) 

y = 0.5975x + 22.142  R² = 0.367   

(B) 

y = 0.5993x + 8.993  R² = 0.2468   

圖15 以本研究薄雲過濾器去雲前後各波段與無雲影像之反射率 XY 散佈圖。左排三張圖代表的 X 軸代 表薄雲影像的反射率值，Y 軸代表無雲影像的反射率值；右排三張圖代表的 X 軸代表去除薄雲後 的影像反射率值，Y 軸代表無雲影像的反射率值。

圖16 以影像分類法進行薄雲去除檢核之示意圖。(A)原始薄雲霧覆蓋影像；(B)薄雲霧影像之分類結果；

(C)去雲後之影像；(D)去雲影像之分類結果。

雲霧分佈區 雲影分佈區 圖17 薄雲霧處理前後分類精度差異探討之分區示意圖

(B) (C)

兩時期NDVI 差值 > 0.1 兩時期NDVI 差值 <= 0.1 (A)

圖18 以 NDVI 進行薄雲去除之檢核結果。(A)原始薄雲覆蓋影像；(B)無雲與薄雲覆蓋 NDVI 影像相減後 之結果；(C)無雲與去雲處理後 NDVI 影像相減後之結果。

4. 結論與建議

本研究最重要的貢獻，可歸納為以下二點。

(1) 影像利用度之提升：本研究僅以單時期以綠、

紅、近紅外波段等三波段影像進行去雲處理的 作業流程，可達到最小成本並符合目前產學界 多數研究的需求。利用本研究所提出的雲層偵 測和還原地物光譜資訊的方法，可應用於土地 利 用 判 釋 及 生 物 量 估 算 中 的 影 像 前 處 理 程 序，減少人工判釋和去除雲層的人力，增加衛 星影像的利用度。以往僅有少數研究探討去雲 過程對原本地物資訊的破壞，既使有也更少探 討去雲後影像能否用來進行自動化地物判釋 則欠缺探討，因此去雲影像對土地利用分類或 是生物量估算等相關研究的使用者來說用途 可能有限。而本研究則探討去雲後的影像「可」

或「不可」做何用途，以及是否在處理完後有 何影像資訊是被扭曲，並傳遞給後端使用者了 解的。

(2) 薄雲空間知識之歸納：相較於以往的研究方 法，如使用紅外光段的雲層偵測、傅利葉轉換 的高通濾波器，雖然皆是可用單時期影像進行 處理，但多需設定經驗閾值，其閾值可能只適 用於特定衛星影像或是特定地區。而本研究成 功歸納厚薄雲在衛星影像中的特性，以量化的 模型來更明確描述過去研究中所提及之雲的 知識。

而本研究之雲層處理在未來仍有一些課題可 再深入研究與改進。例如本研究目前實作的素材仍 是以單一影像中僅有單種雲層的案例，因此在多種 雲層同時出現在一張影像中時，本研究的各種雲霧 處理順序要如何決定，是未來將持續試作及探討的 部分。因為在現實的案例影像中，多為厚雲及薄雲 同時出現的狀況，如何同時精準描繪出無法進行還 原的厚雲區域以及還原受薄雲霧污染區域的地物 光譜值是各領域所關心的。另外，本研究的薄雲去 除影像若和去除前影像進行變遷分析，即可用以估 計雲層厚薄程度並進而估計薄雲干擾的潛在誤 差，但此方法所估計出來的誤差仍是以像元為基礎

(pixel-based) 的概念來呈現。由於高解析度影像資 源愈來愈豐富，物件導向式 (object-based) 分類法 的使用也愈來愈廣泛，因此未來可探討薄雲的影響 程度，在物件導向式分類法以及以像元為基礎的分 類法之間有顯著的差異。而雲影偵測方面，目前雲 影偵測雖已利用雲、雲影、太陽以及衛星的幾何關 係，由程式自動推算雲與雲影的相對位置，給予使 用者建議的閾值。但若是在較複雜的地形或是較大 範圍的雲層所產生的雲影，估算的閾值就需以半自 動之方式調整，偵測結果也可能會與地形效應造成 的 陰 影 產 生 混 淆 。 未 來 可 加 入 數 值 高 程 模 型 (digital elevation model, DEM)，或是判讀雲的類 別，來了解雲可能的高度，另外亦需深入探討地表 的特性、影像中的地物組成成份等影響因素。最後 是雪及沙漠區雲層處理的探討，因為雪及沙漠在可 見光及近紅外光之反射率與厚雲相近，易於混淆，

而本研究又是以綠、紅及近紅外光三波段進行厚雲 處理，因此後續需取高緯以及沙漠地區的影像進行 分析，來探討本研究在這些地區的適用性與限制 性。

參考文獻

內政部國土測繪中心 (2008) 國土利用調查成果 資 訊 網 http://lui.nlsc.gov.tw/LUWeb/. [2008, December 10]

內 政 部 營 建 署 (2008) 國 土 監 測 計 畫 http://www.cpami.gov.tw/chinese/index.php?opt ion=com_content&view=article&id=688&Itemi d=76. [2008, December 10]

徐 逸 祥 (2006) 遙 測影 像之 雲 層 偵 測及 干 擾 去 除，國立臺灣大學地理環境資源學研究所碩士 論文。

徐逸祥、朱子豪、劉英毓 (2006) 衛星影像的雲霧 偵測及干擾去除，國土資源遙感期刊，3：

23-28。

張立雨、陳繼藩、陳哲俊、林欣穎 (2007) 應用影 像 分 割 技 術 與 碎 形 理 論 於 福 衛 二 號 Quick-Look影像之雲覆蓋萃取，航測及遙測學

刊，12 (3)：273-281。

曾筱婷 (2005) 以多時段的 SPOT 衛星影像做雲 層自動去除，國立中央大學資訊工程研究所碩 士論文。

葉精國 (2007) 自動化衛星影像雲及雲陰影遮蔽 區偵測及修正之研究，國立中興大學土木工程 學研究所碩士論文。

蔡博閎 (2009) 衛星影像雲層遮蔽區域之移除與 填補演算法，國立成功大學測量及空間資訊學 系碩士論文。

賴 格 英 、 劉 春 燕 、 辜 曉 青 (2003) 基於 ERDAS IMAGINE 的遙感圖像去雲方法。

Ahmad, A. and Hashim, M. (2002) Determination of haze using NOAA-14 AVHRR satellite data,

MACRES bulletin: 15-27.

Bantges, R. J., Russell, J. E. and Haigh, J. D. (1999) Cirrus cloud top-of-atmosphere radiance spectra in the thermal infrared, Journal of Quantitative

Spectroscopy and Radiative Transfer, 63:

487-498.

Berendes, T. A., Berendes, D. A., Welch, R. M., Dutton, E. G., Uttal, T. and Clothiaux, E. E.

(2004). Cloud cover comparisons of the MODIS daytime cloud mask with surface instruments at the North Slope of Alaska ARM site, IEEE

Transactions on Geoscience and Remote Sensing,

42(11), 2584−2593.

Bsaibes, A., Courault, D., Baret, F., Weiss, M., Olioso, A., Jacob, F., Hagolle, O., Marloie, O., Bertrand, N., Desfond, V. and Kzemipour, F. (2009) Albedo and LAI estimates from FORMOSAT-2 data for crop monitoring, Remote Sensing of

Environment, 113: 716-729.

Chander, G. and Markham, B. (2003) Revised Landsat-5 TM radiometric calibration procedures and postcalibration dynamic ranges, IEEE

Transactions on Geoscience and Remote Sensing, , 41: 2674-2677.

Chen, W., Zhang, Z., Wang, Y. and Wen, X. (2009).

Atmospheric Correction of SPOT5 Land Surface Imagery, CISP '09. 2nd International Congress

on Image and Signal Processing, 1-5.

Dong, J., Kaufmann, R. K., Myneni, R. B., Tucker, C.

J., Kauppi, P. E., Liski, J., Buermann, W., Alexeyev, V. and Hughes, M. K. (2003) Remote sensing estimates of boreal and temperate forest woody biomass: carbon pools, sources, and sinks,

Remote Sensing of Environment, 84: 393-410.

Doraiswamy, P. C., Hatfield, J. L., Jackson, T. J., Akhmedov, B., Prueger, J. and Stern, A. (2004)

Crop condition and yield simulations using Landsat and MODIS, Remote Sensing of

Environment, 92: 548-559.

Du, Y., Guindon, B. and Gihlar J. (2002) Haze detection and removal in high resolution satellite image with wavelet analysis, IEEE Transactions

on Geoscience and Remote Sensing, 40 (1):

210-217.

Feng, C., Ma, J., Dai, Q. and Chen, X. (2004) An improved method for cloud removal in ASTER data change detection, 2004 IEEE International

Geoscience & Remote Sensing Symposium,

3387-3389.

Foody, G.M. (2004) Thematic map comparison:

evaluating the statistical significance of differences in classification accuracy,

Photogrammetric Engineering and Remote Sensing, 70, 627-633.

Hagolle, O., Huc, M., Pascual, D.V., & Dedieu, G.

(2010) A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENµS, LANDSAT and SENTINEL-2 images, Remote

Sensing of Environment, 114, 1747-1755.

Han, K. S., Champeaux, J. L. and Roujean, J. L.

(2004) A land cover classification product over France at 1 km resolution using SPOT4/VEGETATION data, Remote Sensing of

Environment, 92:52-66.

Irish, R.R., Barker, J.L., Goward, S.N., & Arvidson, T. (2006) Characterization of the Landsat-7 ETM_Automated Cloud-Cover Assessment (ACCA) Algorithm, Photogrammetric

Engineering & Remote Sensing, 72, 1179-1188.

Kaufman, Y.J. and Tanre, D. (1992) Atmospherically resistant vegetation index (ARVI) for EOS-MODIS, IEEE Transactions on Geoscience and Remote Sensing, 30, 261-270

Le Hégarat-Mascle, S. and André, C. (2009) Use of Markov random fields for automatic cloud/shadow detection on high resolution optical images,

ISPRS Journal of Photogrammetry and Remote Sensing, 64:

351-366.

Mannozzi, L., Giuseppe, F. D. and Rizzi, R. (1999) Cirrus cloud optical properties in far infrared,

Physics and Chemistry of the Earth (B), 24 (3):

269-273.

Melgani, F. (2006) Contextual reconstruction of cloud-contaminated multitemporal multispectral images, IEEE Transactions on Geoscience and

Remote Sensing, 44, 442-455.

Moro, G. D. and Halounova, L. (2007) Haze removal for high-resolution satellite data: a case study,

International Journal of Remote Sensing, 28:

2187-2205.

Otsu, N. (1979). A threshold selection method from gray-level histograms, IEEE Transactions on

Systems, Man, and Cybernetics, 9: 62-66.

Peterson, D., Whistler, J., Bishop, C., Egbert, S. and Martinko, E. (2009) The Kansas next-generation land use/land cover mapping initiative, ASPRS

2009 Annual Conference, Baltimore, Maryland.

Rakwatin, P., Takeuchi, W. and Yasuoka, Y. (2009) Restoration of Aqua MODIS Band 6 Using Histogram Matching and Local Least Squares Fitting, IEEE Transactions on Geoscience and

Remote Sensing, 47, 613-627.

Ratle, F., Camps-Valls, G. and Weston, J. (2010) Semisupervised Neural Networks for Efficient Hyperspectral Image Classification, IEEE

Transactions on Geoscience and Remote Sensing,

48, 2271-2282.

Richter, R. (1996) Atmospheric correction of satellite data with haze removal including a haze/clear transition region, Computers and Geosciences, 22 (6): 675-681.

Richter, R. (1997) Correction of atmospheric and topographic effects for high spatial resolution satellite imagery, International Journal of

Remote Sensing, 18 (5): 1099-1111.

Roy, D.P., Ju, J., Lewis, P., Schaaf, C., Gao, F., Hansen, M. and Lindquist, E. (2008) Multi-temporal MODIS–Landsat data fusion for relative radiometric normalization, gap filling, and prediction of Landsat data, Remote Sensing

of Environment, 112, 3112-3130.

Stowe, L. L. (1984). Evaluation of NIMBUS-7 THIR CLE and air-force 3-dimensional nephanalysis estimates of cloud amount, Journal of

Geophysical Research, 89: 5370−5380.

Tan, K., Piao, S., Peng, C. and Fang, J. (2007) Satellite-based estimation of biomass carbon stocks for northeast China's forests between 1982 and 1999, Forest Ecology and Management, 240, 114-121.

Wang, Y. P., Chang, K. W., Chen, R. K., Lo, J. C. and Shen, Y. (2010) Large-area rice yield forecasting using satellite imageries, International Journal of

Applied Earth Observation and Geoinformation,

12: 27-35.

Wang, Z., Jin, J., Liang, J., Yan, K. and Peng, Q.

(2005) A new cloud removal algorithm for multi-spectral images, Proceedings of SPIE, 6043: 60430W-1-60430W-11.

Yang, Y., Di Girolamo, L., and Mazzoni, D. (2007) Selection of the automated thresholding algorithm for the Multi-angle Imaging

SpectroRadiometer Radiometric Camera-by-Camera Cloud Mask over land,

Remote Sensing of Environment, 107: 159-171.

Zhang, Y., Guindon, B. and Cihlar, J. (2002) An image transform to characterize and compensate for spatial variations in thin cloud contamination of Landsat images, Remote Sensing of

Environment, 82:173-187.

Zhong, H., Shi, R., Liu, C., Gao, W. and Qu, P. (2010) Thin cloud removal of MODIS imagery based on AERONET data, Proceedings of SPIE - The

International Society for Optical Engineering,

7809.

1

Doctor, Department of Geography, National Taiwan University

Received Date: Feb. 21, 2012 2

Professor, Department of Geography, National Taiwan University

Revised Date: May. 02, 2012 2 Professor, Department of Information Management, Ling Tung University

Accepted Date: Jan. 21, 2013

*.

Corresponding Author, Phone: 886- 2-33665830 ext.12, E-mail：[email protected]

The Study on Cloud Processing in Optical Satellite Imagery

Yi-Shiang Shiu

¹

Tzu-How Chu

^2*

ABSTRACT

Cloud cover is an inevitable interference when mapping land use/cover with optical satellite imagery. In this study, we apply region growing processing to delineate unrecoverable thick cloud and use Fourier analysis to recover ground information from hazy areas with single temporal imagery. Several methodologies across literature successfully solve cloud problems, but most methods require additional cloud-free reference areas or imagery, which may be unavailable in the real world. Moreover, visual methods rather than quantitative methods are used for assessing results, which can be subjective and arbitrary. Most importantly, the feasibility of applying haze-off imagery to image classification is seldom discussed.

To overcome the existing limits, expert method is applied to assess the thick cloud delineation and image classification and normalized difference vegetation index (NDVI) is used to evaluate the recovery degree of ground information after the haze-off processing for quantitative verification of the results. This study revises the image enhancement and region growing algorithm to delineate unrecoverable thick cloud. Accuracy assessment shows the overall accuracy of delineation could be 90% above in each study area. For hazy areas, Fourier analysis is used to reduce haze interference and recover ground information. The proposed haze filter increases the overall accuracy of the whole scenes by about 4%. The overall accuracy of hazy areas in the imagery increases the most (by 6%), while that of shadow areas decreased slightly. The influence of haze on NDVI is also reduced with statistical significance (p < 0.01). Both thick cloud and hazy areas processing can be achieved with no cloud-free area or reference imagery required. Future applications include preprocessing of satellite imagery in land use/cover mapping, which can decrease the manpower to interpret and remove cloud areas and increase the usability of the satellite imagery.

Keywords:

Cloud removal, land features interpretation, region growing, Fourier analysis