應用影像分割技術與碎形理論於福衛二號 Quick-Look 影像之雲覆蓋萃取

(1)

Volume 12, No.3, September 2007, pp. 273-281

應用影像分割技術與碎形理論於福衛二號 Quick-Look 影像之雲覆蓋萃取

張立雨

¹

陳繼藩

²

陳哲俊

³

林欣穎

⁴

摘要

對於應用光學衛星影像進行地面各項土地資源調查而言，影響影像可用度之最大因素即為影像之含雲量。因此對於所接收到之影像若能快速的提供相關雲覆蓋量之評估，不但對於影像使用者在選擇影像上有更良好之參考，也能在影像接收上提供適當資訊以作為未來接收規劃之依據。目前，福衛二號衛星影像在雲量評估上主要是利用系統所產生之 QUICK-LOOK 影像進行人工判釋，而 QUICK-LOOK 影像為衛星影像接收系統產生作為影像檢索用途之影像。而應用 QUICK-LOOK 影像進行雲覆蓋判釋時，主要是先將原始 QUICK-LOOK 影像分割為八等分，並且在每等分中以人工方式給定五種等級之雲覆蓋指標來標示雲含量之多寡。此項指標除必須花費較多時間進行人工判釋外，並且其結果常會受到操作者之主觀認定而受影響。同時，僅使用五種等級所產生之雲覆蓋資訊對於實際應用上亦略顯不足。因此本研究提出利用影像分割並配合碎形理論之分析方法對 QUICK-LOOK 影像進行定量且自動化之雲覆蓋評估。方法上主要先利用區塊成長進行影像分割後，再利用非監督性分類法 K-MEAN 將分割後之區塊進行非監督性分類來產生訓練雲區區塊。最後利用所產生之訓練雲區區塊計算雲區之碎形特性，並以此碎形特性進行所有雲區區塊之萃取。研究成果顯示，將所產生之結果轉換為相當人工判釋所得之雲量指標後再與人工判釋之結果比較，其結果顯示兩者具有一致性。

關鍵詞：QUICK-LOOK 影像、雲區覆蓋量評估、碎形理論

1. 簡介

目前資源衛星接收站提供了國內外各機構與政府單位所需的衛星影像資料以進行廣泛的科學研究與監測應用。而在這些相關研究計畫與應用系統的發展中，就進行地面各項土地資源調查的資料供給面而言，最關鍵的一項的瓶頸即為影像所能提供的週期。而由於台灣地處多雲雨的亞熱帶地區，

所以就應用光學影像於地表面之資訊取得而言，雲霧為影響可用影像之獲取週期的重要因素，若以過去統計數字而言，台灣地區接收到一張無雲影像的

機會僅有 10%。同時，依照過去影像使用之經驗與統計，以 SPOT 系列資源衛星提供涵蓋全台灣地區地表面之衛星影像，其合理之週期約在二至三個月左右。因此若可用影像提供的週期能夠有效的縮短，則各項應用成果在質與量上必然可以有顯著的提升，而其中關鍵性的因素則為衛星影像檢索系統的發展。一個有效率的檢索系統與流程，不但可以提升影像資訊的可用度，並可同時輔助衛星影像的接收規劃與運作，可說是在珍貴的資源衛星資料獲取上，同時達到節流與開源之目的。

對於影像查詢應用而言，影像品質分析以及評

1國立中央大學太空及遙測研究中心助理研究工程師

2國立中央大學太空及遙測研究中心與土木工程學系副教授

3國立中央大學太空及遙測研究中心教授

4國立中央大學太空及遙測研究中心助理研究工程師

收到日期:民國 96 年 02 月 27 日修改日期:民國 96 年 10 月 30 日接受日期:民國 96 年 10 月 30 日

(2)

估系統的建立是相當重要的一環。對於衛星影像來說，影像本身之雲量即為一個重要參考指標。目前現行的系統大多以人工視覺之方式進行雲量評估並給定相對之品質參數。此一過程不但耗費時間人力，並且可能因判釋之操作人員不同而可能有標準不一致的情形發生。在過去的研究中，已有不少專家學者針對如何自動化進行影像中之雲區之分類與分析提出解決的方法，例如應用類神經網路系統進行衛星影像之雲區區塊之分類(Bankert and Aha, 1996)，或是使用模糊邏輯方法針對衛星影像上不同種類之雲區區塊建立模糊歸屬函數後再進行分類(Baum et al., 1997)，其他亦有利用影像質地分析 (Image Texture Analysis)進行衛星影像雲區區塊分類與分析(Christodoulou et al., 2003)。同時，利用影像處理技術進行雲量評估並給定相對之品質參數，不但可有效降低人力之需求，並可提高評估品質之穩定性。

在接收站的建置中，QUICK-LOOK 影像處理系統主要的功能是提供在影像查詢階段的瀏覽資訊，其所產生之 QUICK-LOOK 影像主要以 JPEG 格式進行影像壓縮。目前資源衛星接收站對於福衛二號衛星影像之影像雲量評估是以人工方式進行，其方式主要將原始 QUICK-LOOK 影像分成八等份，

之後再將每等份分別給定由 A~E 五種雲量等級來界定每等份之含雲量(EADS, 2003)，其等級與含雲量之百分比如表 1 所示，而 QUICK-LOOK 影像之八等份分割與雲量評估之範例則如圖 1 所示。

表 1、雲量等級與含雲量之百分比之對照表等級含雲量之百分比 (CP)

A CP = 0%

B 0% < CP≦10%

C 10% < CP≦25%

D 25% < CP≦75%

E 75% < CP≦100%

(a) (b) (c)

圖 1、QUICK-LOOK 影像之八等份分割與雲量評估之範例. (a)一張八等份分割的 QUICK-LOOK 影像 (b) 由人工對每一等份影像進行雲量評估後之結果 (c) 每一等份之雲量最後以等級來表示並儲存於資料庫中。

本研究為了有效降低影像雲量評估之人力需求以及提高評估品質，因此提出兩階段之方法來自動化的針對 QUICK-LOOK 影像萃取雲區區塊。首先在第一階段中利用區塊成長(Adams and Bischof, 1994)進行影像分割後，再以非監督性之分類法 K-MEAN(Tou and Gonzalez, 1974)將分割後之區塊

進行非監督性分類來產生訓練雲區區塊。之後於第二階段中再利用前步驟所產生之雲區區塊進行雲區之碎形特性分析，並以所得之雲區碎形特性進行所有雲區區塊之萃取。為了進一步說明區塊成長以及 K-MEAN 分類法的演算過程，在此敘述如下：

(1) 區塊成長法

(3)

區塊成長法(region growing)主要是要將影像上光譜值相似的像元群聚成區塊，主要可分為三個階段：

A. 選擇起始種子點

首先，需設定區塊分割的條件作為分割合併之依據，如光譜距離門檻。接著，選擇某個像元當作起始種子點(seed pixel)，種子點之選擇可由影像的任何位置出發，一般而言種子點均為亂數產生，並給予種子點像元一區塊編號(label)。選定種子點後便以此種子點之光譜值作為參考光譜值，並將種子點周圍鄰近四方向或八方向之像元光譜值與參考光譜值進行光譜距離之計算，若所得之光譜距離小於門檻，則此像元則指定為與種子點相同之區塊編號，並同時以相同區塊編號之像元光譜平均值更新參考光譜值。

B. 區塊成長

將起始種子點周圍四方向或八方向偵測得到相同區塊編號之像元，當作新的種子點出發，如同上一步驟中對於新種子點周圍四方向或八方向之像元進行偵測工作。若新種子點周圍仍有光譜距離小於門檻之像元產生，則將新產生之像元再當做新的種子點繼續向外偵測。此一向外迭代、增長之過程稱為區塊增長。區塊不斷的增長直到周圍無法找到光譜距離小於門檻的像元，則此一區塊建置完成。

C. 區塊分割完成

當一個區塊建置完成後，我們對於影像上未分割區域之像元，再選擇一個種子點，作為下一個區塊增長之核心，並以此進行下一個區塊之增長工作。接著，當影像上所有的像元都有其歸屬之區塊後，還需將面積過小的區塊像元與其鄰近之區塊進行合併之工作，而合併之方式則可由光譜差異判斷。最後，影像區塊分割即為完成。

(2) K-Mean 分類法為一非監督式之分類法，其分類過程如下：

A. 在影像光譜空間中選取與數目與分類類別數相等之數個初始中心作為類別中心。

B. 計算影像中所有像元與數個類別中心之光

譜距離，並將每個像元之類別指定到光譜距離最短的類別上。

C. 按照每個像元在前步驟所指定之類別，以所有同類別之像元光譜平均值重作為新的類別中心。

D. 重複步驟 B 進行迭代，直到所有或是預設百分比的像元其類別不再改變，亦或是迭代次數超過上限時則停止分類。

2. 研究方法

在本研究中所發展之模式主要可分為兩個階段：訓練階段以及分析階段，其流程如圖 2 所示。

(1) 訓練階段

在訓練階段中首先以區塊成長之方法對輸入 QUICK-LOOK 影像進行影像分割。一般而言，雲區區塊在全色態或多光譜影像中均呈現相當高之光譜值，亦即其在光譜上和地面之類別在統計上有顯著不同之分布(Gotoh and Fujii, 1998)。圖 3 為一典型多光譜自然色波段並含有雲區區塊影像之 Scatter Plot，雲區大多產生在圖中所示灰度較為頂端之位置。因此在訓練階段中採用非監督性之 K-MEAN 分類法針對分割後之區塊光譜平均值進行分類，唯在分類初始時給定之初始類別中心分別於 Scatter Plot 上(0,0,0)以及(255,255,255)此二位置。

由於雲區類別和地面類別在統計上有顯著不同之分布，因此在 K-MEAN 迭代之過程中，(255,255,255) 此一中心理論上應會逐漸收斂至雲區之光譜中心，而(0,0,0)此一中心則會收斂至其他地面類別之光譜中心，最後即可得到兩類別之輸出。

(2)分析階段

在第二階段進行雲區碎形特性分析前，必需先行判斷前階段所產生之雲區品質後才可視其為訓練雲區：

z 雲區類別與其他地面類別是否有任一類別之像元總數為零

YES：影像全為像元總數不為零之類別

NO：兩類別中心距離是否有足夠之分離度

(4)

z YES：雲區類別為可靠之雲區

取雲區類別區塊中，區塊各波段光譜平均值較接近(255,255,255)之前 1/2 個區塊作為訓練雲區並以碎形理論分析其他所有之雲區區塊

z NO：雲區類別不是可靠之雲區

影像各波段平均值之總平均＞ Threshold

影像均為雲區類別

影像各波段平均值之總平均 ≦ Threshold

影像均為其他地面類別

在判斷前階段所產生之雲區品質中，為了取得最具有代表性之雲區區塊進行後續碎形理論分析，因此僅選取雲區類別區塊中，各波段光譜平均值較接近(255,255,255)之前「1/2」個區塊作為訓練樣本。在經過上述之條件判斷後即可進入分析階段，其中關於類別分離度以及 Threshold 在本研究中是以經驗值設定之。分析階段主要應用訓練階段所得之雲區區塊進行碎形特性分析後，再以此來產生所有之雲區區塊。根據雲區所具有之自我相似 (Self-Similarity)特性，雲區區塊之碎形特性應為近似，過去文獻在關於應用碎形理論進行雲區之辨識與萃取亦有相當之討論(Carvalho and Silva Dias, 1998; Gotoh and Fujii, 1998)。在本研究中主要是以 Box Counting Fractal Dimension(Peitgen et al., 2004)的觀念來分析雲區區塊碎形關係。理論上，具有自我相似特徵之物件間在其度量基準之尺度Δs 上與其對物件進行度量所得之總量 N(Δs)應呈正比之關係，如式(1)所示：

s D

c s

N(Δ )≈ (1/Δ ) (1) 其中c 為常數。而式(1)可利用對數之運算關係改寫為式(2)：

[

( )

]

log(1/ ) '

logN Δs ≈D Δs +c (2) 式(2)說明了尺度Δs與總量N(Δs)在對數平面上為一線性關係，其中c'=log(c)。此外，由於福衛二號影像之單一像幅涵蓋範圍為24公里x24公里，此範圍就雲區特性變化之空間尺度而言相對較

小，因此在本研究中假設在單一像幅中其雲區之特性相同，換言之即為在單一像幅中之雲區碎形特性應相同。另一方面，大的雲區區塊在影像上應可呈現出較小的雲區區塊有更多的變化，因此雖然大的雲區區塊與小的雲區區塊在同一幅影像上其解析度是相同的，但若將大的雲區區塊等比例縮小至小的雲區區塊的大小時，應可視大的雲區區塊是以

「較小度量基準之尺度Δs」所度量到之結果，也就是說假設將所有不同體積之雲區都放大縮小至相同的比例時，此時因所有的雲區將會有近似的碎形特性，因此對每一個區塊可以度量到的N(Δs)值即可由(3)式表示之

∑

=

= Δ ⁿ

i i

h R h s N

1 max

)

( (3)

此時n為區塊總像元數。在度量高度方面則以區塊中每個像元之各波段光譜值之平均值計算，因此每個像元的各波段光譜值之平均值為hⁱ，而各波段光譜值之平均最大值為h^max，此值在BYTE影像上應為255。此外由於在做比例的放大縮小時，在高度方面由於不論大小雲區其光譜值變化範圍均相同，因此就較大的雲區在高度上(指使用各波段光譜值之平均值時)應該要增加的分割就以乘上比例調整值R來調整之，並且本研究中之比例調整值R 則以雲區底面面積之平方根來近似之，也就是n^1/2。換言之，亦即所有的區塊都要放大縮小至相同的比 例後再做比較。而(3)式可改寫為(4)式，其中 h 表示 區塊之灰度平均值。

max 2 3

max 1

max

1 max

) 1

( h

n h h Rn h n h

h R n h R h s N

n

i i n

i

i = × = =

=

Δ

∑ ∑

=

(4)

此時之度量尺度Δs 在本研究中則以前述之比例調整值來定義，如(5)式所示

2 /

/ 1

1 = ⁻

=

Δs R n (5) 將式(4)及式(5)以對數表示可得式(6)及式(7)：

[ ]

n

h s h

N log

2 log 3

) ( log

max

+

=

Δ (6)

n s log

2 ) 1 / 1

log( Δ = (7)

(5)

此時我們可將前一階段所產生之訓練雲區區塊利用式(6)及式(7)進行式(2)之回歸計算，此時將可得到可以代表雲區區塊之碎形特徵直線之參數 D 及 c＇。

如圖 4 所示，理論上影像分割所得之其他區塊若其為雲區時，其區塊應有之尺度Δs 與其體積總量 N(Δs)在前述所定義之對數平面上(如圖 4 中之綠色點所示)，應該相當靠近代表雲區區塊之碎形

特徵直線(如圖 4 中之黃色直線)。因此在本研究中以代表雲區區塊碎形特性之直線為中心，取相互間距離為Δ且斜率與碎形特性直線相同的兩條直線作為限制 (如圖 4 中之橘色虛線)，最後將落於此限制範圍內之所有區塊作為判釋所得之雲區區塊。此外，門檻值Δ是在 log[N(Δs)]方向上計算，並且是以訓練雲區回歸時所產生之 log[N(Δs)]之值的最大殘差絕對值再向上及向下平移 1.5 倍來計算。

圖 2、本研究之雲區區塊萃取流程

圖 3、典型自然色波段含有雲區區塊之 Scatter Plot，雲區大多在圖中所示之頂端位置(Gotoh and Fujii, 1998)。

(6)

圖 4、雲區區塊之碎形特徵直線與其他區塊在對數平面上之關係

(a) (b)

(c) (d)

圖 5、應用福衛二號 QUICK-LOOK 影像進行測試之成果。圖(a)及圖(b)分別為獲取時間 05/13/2005，位置 Path/Row 為 587/330 的 QUICK-LOOK 原始影像以及分析所得的雲區套疊成果，圖(c)及圖(d)分別為獲取時間 09/03/2005，位置 Path/Row 為 584/342 的 QUICK-LOOK 原始影像以及分析所得的雲區套疊成果。其中黃

色區塊為第一階段所得雲區區塊，紅色區塊為經第二階段碎形分析後所得雲區區塊。

(7)

3. 研究結果

圖 5 為應用福衛二號 QUICK-LOOK 影像進行測試之成果。由其結果可看出利用本研究所提方法偵測出之雲區區塊在定性上與視覺的觀察具有一致性。

為了更進一步說明本方法對於一般 QUICK-LOOK 影像之雲區區塊之偵測具有一般性與代表性，因此利用 2005 全年度福衛二號多光譜

衛星 QUICK-LOOK 影像作為測試，並且針對現有系統所產生之雲量指標與本研究所開發之模式加以比較。此一部份總共測試至 2855 幅影像，共 22840 個雲量指標值(每幅影像有 8 個雲量指標值)。比較時須先將本研究所開發之模式所產生之網格雲區判釋結果轉換為相當原系統判釋所得之雲量指標，之後再與原系統之結果比較。圖 6 為經此一致性分析後統計所得之直方圖，表 2 為在一致性分析後，雲量指標誤差在數量上之累積百分比。

0 2000 4000 6000 8000 10000 12000 14000 16000 18000

-4 -3 -2 -1 0 1 2 3 4

Grade Difference

Count

圖 6、經一致性分析後統計所得之直方圖

表 2、一致性分析後，雲量指標誤差在數量上之累積正確百分比

指標等級差異資料筆數累積正確百分比

0 15667 68.6%

±1 20949 91.7%

±2 22108 96.8%

±3 22663 99.2%

±4 22840 100.0%

由圖 6 可觀察到，相較於原人工判釋之結果，

模式所得之雲量有微幅之低估，分析其原因為由於雲區大多屬於高亮度影像區塊，此類區塊對於人工

判釋上常會受視覺之影響而有範圍擴大之錯覺，因此會造成人工估算上之過度估計所致，因此若視雲量指標誤差 ±1 等級為可接受之誤差範圍時，由表

(8)

2 可得知本模式之精確度可達 91.7%。

4. 結論

由實驗結果顯示，本研究所發展模式應用於福衛二號多光譜 QUICK-LOOK 影像之雲區區塊萃取時，所產生之雲區區塊與實際之雲區具有足夠之代表性。此外，將所產生之結果轉換為相當人工判釋所得之雲量指標後再與人工判釋之結果比較，其結果亦顯示兩者具有相當之一致性。若考慮本模式之結果與人工判釋之雲量指標相差 ±1 個指標視為可接受之範圍時，其正確率可達 91.7%。對於少部分雲區覆蓋區域無法正確偵測之問題，未來將會在方法上作進一步之修正。

參考文獻

Adams, R. and L. Bischof, 1994, Seeded region growing, IEEE Transactions on Pattern Analysis and Machine Intelligence, 16(6), pp. 641-647.

Bankert, R. L. and D. W. Aha, 1996, Improvement to A Neural Network Cloud Classifier, Journal of Applied Meteorology, 35, pp. 2036-2039.

Baum, B. A., V. Tovinkere, J. Titlow, Hnd R. M. Welch, 1997, Automated Cloud Classification of Global AVHRR Data Uusing A Fuzzy Logic Approach, Journal of Applied Meteorology, 36, pp.

1519-1540.

Carvalho, L. V. and M. F. Silva Dias, 1998, An Application of Fractal Box Dimension to the Recognition of Mesoscale Cloud Patterns in Infrared Satellite Images, Journal of Applied Meteorology, 37(10), pp. 1265-1282.

Christodoulou, C. I., S. C. Michaelides and C. S.

Pattichis, 2003, Multifeature Texture Analysis for the Classification of Clouds in Satellite Imagery, IEEE Transactions on Geoscience and Remote Sensing , 41(11), pp. 2662-2668.

EADS, 2003, IAPS Operator manual, EADS ROCSAT-2 IAPS Terminal technical document:

R2T-MO-1-00020-SDE, Paris.

Gotoh, Kazuo and Y. Fujii, 1998, A Fractal Dimensional Analysis on the Cloud Shape Parameters of Cumulus over Land. Journal of Applied Meteorology, 37(10), pp. 1283-1292.

Peitgen, Heinz-Otto, H. Jürgens and D. Saupe, 2004, Chaos and Fractals - New Frontiers of Science, 2nd Ed. Springer Verlag, New York.

Tou, J. T. and R. C. Gonzalez, 1974, Pattern Recognition Principles, Addison-Wesley, Massachusetts.

(9)

Using Segmentation Techniques and Fractal Theory to Extract Cloud Coverage for Formosat-2

QUICK-LOOK Images

L. Y. Chang

¹

C. F. Chen

²

A. J. Chen

³

H. Y. Lin

⁴

ABSTRACT

The percentage of cloud coverage is generally considered the most important factor to affect the availability of optical satellite images for land resource applications. Consequently, an effective cloud assessment process can not only offer general users a guideline for image browsing, but also provide valuable information for orbit planning of ground receiving station. At present, the cloud coverage of Formosat-2 data is visually examined and manually estimated by the operators on the QUICK-LOOK images. The QUICK-LOOK images basically are the reduced images created by satellite data receiving system for browsing purpose. As to the current estimation of cloud coverage, every QUICK-LOOK image is divided into 8 portions and the cloud coverage for each portion is indicated with a cloud index with 5 grades. Apparently, the visual and manual cloud assessment method is labor-intensive, time-consuming, and dependent on the operator’s experience. Furthermore, the existing 5-grade cloud index is inadequate to describe the accurate percentage of cloud coverage. Thus, there is an urgent need for developing an automatic approach for cloud assessment. In this study, a two-stage scheme based on image segmentation and fractal theory is used to automatically extract cloud pixels from QUICK-LOOK images. The experimental results indicate that the cloud areas extracted by proposed scheme are very similar to the manual interpretation. In addition, when the automatically extracted cloud coverage is transferred to the existing cloud index, a good agreement between our automatic approach and current manual estimation is achieved.

Key Words: QUICK-LOOK image, Cloud assessment, Fractal theory.

1 Assistant research engineer, Center for Space and Remote Sensing Research, National Central University.

2 Associate professor, Center for Space and Remote Sensing Research and Department of Civil Engineering, National Central University.

3 Professor, Center for Space and Remote Sensing Research, National Central University.

4 Assistant research engineer, Center for Space and Remote Sensing Research, National Central University.

Received Date: Feb. 27, 2007 Revised Date: Oct. 30, 2007 Accepted Date: Oct. 30, 2007