航 航 測 測 及 及 遙 遙 測 測 學 學 刊 刊
Journal of Photogrammetry and Remote Sensing
發行人：陳良健
出版者：中華民國航空測量及遙感探測學會
地址：台北市文山區羅斯福路五段 113 號三樓 信箱：台北市郵政 93-158 號信箱
電話：886-2-8663-3468 886-2-8663-3469 傳真：886-2-2931-7225
電子信件：[email protected] 網址：http://www.csprs.org.tw
PUBLISHER: L. C. Chen
PUBLISHED BY: Chinese Society of Photogrammetry and Remote Sensing
Address: 3F, No.113, Sec.5, Roosevelt Road, Taipei, Taiwan Mail Address: P. O. Box. 93-158, Taipei, Taiwan
Tel: 886-2-8663-3468 886-2-8663-3469 Fax: 886-2-2931-7225
E-mail：[email protected] Web Site：http://www.csprs.org.tw 總編輯
曾義星
國立成功大學測量及空間資訊學系 電 話：886-6-275-7575 分機 63835 傳 真：886-6-237-5764
電子信件：[email protected]
EDITOR-IN-CHIEF Yi-Hsing Tseng
Department of Geomatics, National Cheng Kung University Tel: 886-6-275-7575 ext. 63835
Fax: 886-6-237-5764
E-Mail: [email protected]
編輯委員 EDITORIAL BOARD
農林
鄭祈全 (文化大學) 申雍 (中興大學)
Agriculture and Forestry
C. C. Cheng (Chinese Culture University) Y. Shen (National Chung Hsing University)
海洋與大氣 劉振榮 (中央大學) 何宗儒 (台灣海洋大學)
Oceanography and Atmospheric Science G. R. Liu (National Central University) C. R. Ho (National Taiwan Ocean University)
防災
劉進金 (達雲科技有限公司) 徐百輝 (台灣大學)
Disaster Management
J. K. Liu (LIDAR Technology Co., Ltd) P. H. Hsu (National Taiwan University)
土地利用與覆蓋 陳永寬 (銘傳大學) 詹進發 (政治大學)
Land Cover and Land Use
Y. K. Chen (Ming Chuan University) J. F. Jan (National ChengChi University)
地理資訊系統
何維信 (中華科技大學) 蔡榮得 (中興大學)
Geographical Information System
W. H. Ho (China University of Science and Technology) J. D. Tsai (National Chung Hsing University)
遙感探測
陳良健 (中央大學) 陳錕山 (中央大學)
Remote Sensing
L. C. Chen (National Central University) K. S. Chen (National Central University)
攝影測量與光達 廖揚清 (成功大學) 史天元 (交通大學) 趙鍵哲 (台灣大學)
Photogrammetry and LiDAR
Y. C. Liao (National Cheng Kung University) T. Y. Shih (National Chiao Tung University) J. J. Jaw (National Taiwan University)
封面圖片說明 About the Cover
本研究利用非監督分類遙測資料及監督學習估算遙測資料兩種方式，再分別加入當日少數定點水質採 樣資料作共克利金內插；另將高程(DEM)移至共克利金內插時作趨勢消除變量，藉以減少非監督分類時間(本 文後，非監督結合共克利金法加入 DEM 簡稱 UCKD。圖中可發現基隆河支流及左下方無測站支流，由於加 入了衛星光譜分類值作內插，因此可分類出符合該支流之污染顏色(數值愈高表示污染愈嚴重；由四個水質 參數項目 DO、BOD、SS、NH3-N 平均計算出河川污染指標(River Pollution Index，RPI))，可避免無測站支 流受主流或其他分流所影響，充分顯示以 UCKD 法建立河川水質污染面化圖之優點。
(封面圖片出處：比較監督及非監督模式與共克利金法估算河川污染指標面化圖，
第十七卷 第二期 第 77-93 頁)
Volume17, No2, August 2013, pp. 77-93
1國立雲林科技大學環境與安全衛生工程系 講師 收到日期:民國 100 年 03 月 28 日
2國立雲林科技大學環境與安全衛生工程所 碩士 修改日期:民國 100 年 09 月 24 日
3國立臺灣大學生物環境系統工程學系 博士後研究員 接受日期:民國 101 年 04 月 21 日
＊通訊作者, 電話: 05-5342601 ext.4489, E-mail: [email protected]
比較監督及非監督模式與共克利金法估算河川污染 指標面化圖
施明倫
1洪志豪
2*楊政儒
2摘 要
本研究擬結合SPOT 衛星遙測影像(Remote sensing, RS)及現場水質採樣分析資料，並以兩階段非監 督模式分類遙測圖結合共克利金法(Co-kriging, CK)內插建立台灣河川流域水質污染指標面化圖(River Pollution Index Map, RPIM)。本研究非監督模式採用自動最佳化分類結果，藉此分類衛星遙測影像涵蓋 河川水體之污染分佈情形，再加入高程，並以當日少數定點水質採樣資料為參考主要資料作共克利金空 間內插，兼顧遙測及現場實測兩種方法的優點建立河川水質污染面化圖，最後再與監督式共克利金法的 估算結果及傳統僅用監測站值空間內插結果作比較。綜合上述觀點非監督共克利金法具有簡單、快速及 改善監督共克利金法學習建模過程複雜的優點，且準確率亦未降低，顯示非監督共克利金法取代監督共 克利金法估算河川未設測站水質污染之優勢。
關鍵詞：遙測(RS)、非監督模式、共克利金(CK)、河川水質污染指標面化圖(RPIM)、
數值高程(DEM)
1. 前言
1.1 研究背景
目前遙測技術已廣泛應用於氣象、農漁、防 災、土地管理及環境監測，遙感探測技術具有快 速、簡便、同步、客觀、經濟等優勢，可以有效 的監測河川表面水質參數空間及時間上的變化。
有別於一般傳統的河川水質監測，採用人工現場 採樣實驗，再求得其河川水質污染指標，不但成 本高、耗時、同步性差，且往往只能獲得河川上 少數測站的水質狀況，並不能有效的掌握整條河 川水質之變化。因此本研究嘗試使用衛星遙測影 像，結合兩階段非監督模式及共克利金法應用於 河川水質污染分類的監測。以有效且半即時的掌 握河川水質污染指標之變化。本研究使用SPOT 衛星影像圖， 衛星影像之幾何扭曲校正，使得衛 星影像在土地利用監測及水體辨識之可行性，大
幅提升河川湖泊水質監測之潛勢。
1.2 文獻回顧
水質遙測之探討文獻，一般可分為指標污染 物參數項目的選擇、光譜變量的選擇及應用的模 式等。本研究文獻回顧共分為五部份，第一部份 提出應用遙測水質指標進行光譜相關性分析之研 究；第二部份探討衛星影像水體辨識採用之方法；
第三部份針對衍生波段的採用及變量篩選之分析；
第四部份則探討非監督模式應用在影像分類之文 獻；及第五部分多變量空間統計分析之應用。
1.2.1 衛星遙測水質指標可行性
近年來遙測技術已達成熟階段，也廣泛的應 用於土地分類及環境監測等，水質遙測方面亦有 相當多的研究，回顧文獻中可得知水質遙測技術 已應用於各種河川、湖泊及水庫污染，遙測的水
質包括優養化，其項目為懸浮微粒、葉綠素、溫 度及透視度(Ritchie et al., 1990；Ritchie et al., 1995；
Yang et al., 1999；Harma et al., 2001；Brivio et al., 2001)，或有機污染物 TOC、BOD、COD(Tripathi
et al., 1998；Wang et al., 2004)，皆具有高相關性
且不錯的估算結果，且由於河川本身為一綜合性 污染，各種污染物質彼此間相互干擾，而台灣常 用河川水質污染指標(River Pollution Index, RPI)為 河川標準，故本研究選用 RPI 作為研究項目，其 參數分別為DO、BOD、SS、NH3-N。1.2.2 衛星影像水體辨識
接著衛星影像選擇部份，由於 SPOT 衛星影 像解析度夠高，適合應用於水質遙測，其影像與 環保署監測站相對位置配合也較好，且以同等級 衛星影像作比較，SPOT 之衛星影像單價亦較低，
加 上 本 研 究 從 早 期 至 目 前 之 水 質 遙 測 皆 使 用 SPOT 衛星影像，為維持資料一致性，因此選用 SPOT 衛星影像作為研究遙測圖資，由於遙測影 像易受大氣及衛星本身系統影響，因此衛星影像 需先經輻射及大氣校正(雷祖強及鄭克聲，2000)，
而目前本研究僅作輻射校正，針對衛星感測器的 失真作修正，並未進行大氣散射校正，其主要藉 由太陽天頂角、校正係數及外太空輻射值的演算 方法，將影像的灰度值(Digital Number, DN)轉換 為反射率，再利用非監督式分類法中之反覆式自 我組織資料分析群集偵測法(ISODATA)(Nellis et
al., 1998)進行河道及非水體分類，接著以人工方
式將其河道擷取出來。1.2.3 衍生波段應用及變量篩選
為了提高非監督模式分類之準確率，本研究 嘗試經由變量篩選得到結果較佳之衍生波段，加 入模式進行分類，而本研究參照去年本實驗室以 STATISTICA 8.0 中 F 檢定進行變量篩選後的結果 (林珈禎，2009)，以 DEM、B3/B2、NDVI 變量組 合為最佳，其 R2達 0.6 以上，但當變量組合過多 時，反而導致非監督模式分類時間增長，故逐步 減少輸入變量以提高非監督分類效率。
1.2.4 非監督模式應用影像分類之 可行性
在河川分類法中，一般有監督及非監督式兩 種分類方法，大多數人採用監督式(如多變量迴歸、
類神經網路及支撐向量機)(Chen et al., 2004；張子 鵬，2006)，而前提必須建構一個具有多樣且代表 性的資料庫來提供學習，但在累積取得水質光譜 資料庫上，尚不足以代表所有可能水體污染，且 監督模式在參數最佳化選取過程非常複雜及繁瑣，
加上本研究樣本數不夠多，故本研究選用步驟較 為簡易之自行開發兩階段非監督方法(顏可翰，
2008)進行河川水體污染分類。
1.2.5 多變量空間統計分析應用
早期本實驗室使用的非監督模式(顏可翰，
2008)僅能最佳化分類河川水體，再經由人工主觀 判別其污染程度，因此往往會挑選有利於自己的 判別，並無法精準呈現出污染範圍值，藉此本研 究欲加入少量水質採樣實測值，並結合非監督分 類之衛星光譜值，綜合現場實測及遙測兩種方法 的優點以共克利金法(Co-kriging, CK)建立河川水 質污染面化圖，結合 CK 法後可改善早期非監督 模式需經由人工主觀判別之缺點，且自行精準估 算 河 川 污 染 範 圍 值 ； 由 於 相 關 研 究(Chappell (1998)；Rivero et al., 2007)顯示共克利金法(CK)中 之 ordinary kriging 估計誤差較小且穩定性高且於 (張永欣及李錫堤，2007)以 2004 年艾利颱風為例 進行雨量空間分佈特性之分析，以改善雨量站密 度不足的問題，故本研究最後以非監督模式結合 共克利金法建立河川污染空間分佈，並同步選用 監督模式結合共克利金法來作為比較模式及距離 平方反比內插法(Inverse Distance Weighted, IDW) 來作為比較模式。
2. 研究地區及基本資料
2.1 研究區域
本研究選用流域為淡水河、高屏溪，選擇條
件上盡量選取河道寬涵蓋流域面積大且水質監測 站較多之河川，而不選擇濁水溪原因為根據經濟 部水利署資料顯示，集集攔河堰設置後，濁水溪 水量減少一半，枯水期幾乎無水流排放，斷流嚴 重，加上監測站少，故僅擇取淡水河、高屏溪作 為本研究區域。
衛星拍攝日期則必須與環保署河川採樣為同 一天，由於河川採樣時間無法使數個測站配合衛 星拍照時間同步執行，故本研究假設河川水質於 數小時內為穩定狀態，並確認當日拍攝影像無任 何雲層遮蔽等干擾，故選擇十張衛星圖，分別為 淡水河七張(2004/11/03、2005/08/02、2007/07/03、
2008/03/04、2009/02/09、2009/05/08、2009/05/29)，
及 三 張 高 屏 溪 (2004/10/14 、 2009/01/03 、 2010/01/02)，共 95 筆測站資料，其中 2009/05/08、
2009/01/03、2010/01/02 為 SPOT 5 解析度較高(10 m)之衛星圖，接著利用 ERDAS IMAGINE 軟體進 行衛星影像前處理，以萃取流域各水質測站鄰近 水體的灰度值，研究位置如圖1 所示。
圖1 研究區域
2.2 水質資料
本研究採用行政院環境保護署監測的水質資 料，利用四個水質參數項目 DO、BOD、SS、
NH3-N 平均計算出河川污染指標(River Pollution Index，RPI)，數值愈高表示污染愈嚴重，即是將 四個水質項目之濃度依其點數進行平均計算，所
得到的點數再根據污染程度積分範圍轉換成未 (稍) 受污染、輕度污染、中度污染及重度污染。
本研究再配合研究區域同一時間之衛星影像作後 續光譜估算水質污染分類面化之研究。
3. 研究方法
本研究首先從行政院環保署取得河川水質參 數資料，做為本研究之輸入資料，再進行河道的 選取及光譜波段分析的評選，以萃取水質監測站 的水質參數及對應之像元波譜資料來判別分類，
並利用交叉驗證來估算其準確度。
3.1 水體像元選取
本研究因針對河川水體的辨識，河川水體寬 度一般較湖泊水庫或海洋範圍小，取樣上較為困 難，如果無法獲得解析度高之衛星影像，往往會 使河川水體取樣受到周圍雜訊之干擾，因此在遙 測影像取樣過程中，水體像元(pixels)的辨識為首 要的關鍵。
首先本研究先以ERDAS IMAGINE 軟體之非 監督式分類法(ISODATA)將水體及非水體部份分 類出來，利用水體光譜反射值比非水體較小之差 異來增加水體分辨度，以提升水體取樣的準確性，
最後再以人工圈選方式擷取河道之水體。
3.2 輻射校正
輻射修正最主要的目的是在校正大氣散射的 影響及影像的失真與雜訊，然又因日照的強弱、
大氣的條件、太陽的高度角、觀測角度與表面反 射率等因素的變化與系統性的儀器誤差，會造成 相同的地物類別具有不同的輻射量，而不同的地 物類別有可能會具有相近的輻射量，如此會造成 影像辨識上的困擾。因此遙測影像在研究之前需 要先經過輻射修正，而輻射修正最主要包括兩種，
衛 星 感 測 器 的 錯 誤 修 正(Correction for Sensor System Detector Error) 及 大 氣 散 射 (Atmosphere Scattering of Radiation) 的修正，而本研究使用模 式為非監督，目前僅使用衛星感測器的錯誤修正，
主要藉由太陽天頂角、校正係數及該衛星於外太 空輻射值的演算方法，將影像的灰度值轉換為反 射率。
輻射校正方法：考慮太陽光輻射角從地面反 射到衛星的反射率
反射率 =
E
ZG DN
cos
0
(1)
式中，DN (Digital Number)為衛星測得之灰度值，
G 是校正係數，θ
Z及E
0分別為太陽天頂角及外太 空的輻射值。3.3 波段篩選
本研究參考去年本實驗室(林珈禎，2009)所 使用之衍生波段組合，藉由STATISTICA 8.0 中之 變異數分析進行變量篩選，輸入各衍生變量及數 值高程變量與 RPI 點數，以變異數分析中之 F 檢 定進行變量篩
選，結果如表 1 所示，以此五個 變量之檢定結果較佳。
表 1 輸入變量與 RPI 點數之 F 檢定(林珈禎，
2009)
變量 F-value P-value DEM 23.14 0.000 B3/B2 9.76 0.000 NDVI 9.01 0.000 B3% 7.98 0.000 B2-B3 6.59 0.000
3.4 兩階段非監督式模糊及機 率群集演算法
本研究所使用非監督模式之演算法分為兩個 階段，第一階段決定群體的個數並選擇候選群中 心（candidate centers），第二階段利用第一階段 所選擇的候選群中心當作起始群中心並分類出結 果。圖 2 為兩階段非監督演算法流程圖。第一階 段可分為：『階層式減少候選群中心演算法』，
將可能的候選群中心挑選出來、『選擇候選群中 心演算法』，則是決定真正的候選群中心、『最
佳化模糊因數演算法』為決定模糊因數大小。第 二階段是利用『EM 群集演算法』分類出最後的 結果。
圖2 兩階段非監督演算法流程圖
3.4.1 階層式減少候選群中心演算法
階 層 式 減 少 候 選 群 中 心 演 算 法(HSC ： Hierarchical Subtractive Clustering algorithm)是參 考(Tao, 2002)的作法，將所有的資料點分割成若 干個子集合，針對每個子集合選出若干個候選群 中心，再將個別子集合所選出來的候選群中心集 合起來成一個新的集合，再由這個新的集合選出 若干個新的候選群中心交給之後的演算法挑選。
N 階層的目的是要減低演算法的運算量，所以可 以依據輸入的資料量來決定需要做幾層的減少，
如果輸入的資料不多的話，可以直接當成一個集 合選出這個階段的候選群中心，但是如果輸入的 資料量非常多的話，可以多增加一層以減少運算 的時間。
假設所有資料的集合為
C = {x
1, x
2,…,x
n}，其 階層式減少候選群中心演算選擇最佳化群數演算
最佳化模糊因子演算
否，m=m+0.1
選擇候選群中心
EM 群集演算法
是，利用距離最小的m 值找
最後群數及群中心位置
m≧mmax 輸入資料
分類結果
中
n 是輸入資料點的個數，而每個點是 p 維的向
量，xi= (x
i1, x
i2, …, x
ip)。在第一層的運算中，先 將集合 C 平均分成 NS 個子集合(在實際的程式裡，設定NS
n )，所以每個子集合的大小最多 只有一個的差距。假設每個子集合的大小為SS
(SS = n / NS)，所以可以定義子集合為NS λ x
x x
C
{
1,
2,...,
SS}, 1
(2) 接著定義『鄰域』(neighbor)，鄰域是指對一 個資料點來說，與它距離小於某個半徑內的點都 是它的鄰域，而這邊所定義的半徑（rmin）是對每 個點的每個維度算出其標準偏差值，取標準偏差 最小的值當作rmin，以數學式來說是min
pr
min {x
1p, x
2p, ..., x
np} )
(3) 另外需要定義密度函數（density function），一個資料點
x 的密度函數是指這個點的鄰域數，
也就是包含
x 且與 x 的距離小於 r
min點數，可用 下列式子表示
b
j
i j i
i x u r x x
D
1
min )
( )
(
(4)
The number of points in the neighborhood of
x
iς其 中
u 是 單位 步 階 函 數（ unit step function ）
0 otherwise.
0, if ) 1
(
u ，對一個子集合要選出候選群
中心，就是要先找到這個子集合中擁有最大密度 函數的點，並取這個點與其鄰域的質心當作候選 群中心，也就是說，假設點xc1
Cx 有最大的密 度函數（ Di
Dc1 maxi ），那麼此子集合所選出來 的第一個候選群中心可以表示為
c
1c
= mean(data points in theneighborhood of
x
c1 ) (5) 第二層HSC 的方法就像第一層的其中一個子 集合一樣，首先把第一層所有子集合所選出來的 候選群中心當成一個新的集合，並且利用找最大密度函數的方式選出若干個新的候選群中心，最 後得到的候選群中心就是這階段所要的結果。
3.4.2 選擇最佳群數及起始群中心 演算法
新目標函數的意義在於當多選到一個群中心 時，點對群中心的能量總和（Jfcm）將減少，而中 心對中心間的能量和（Jcc）將增加，所以當點對 群中心減少的能量和大於中心對中心增加的能量 和時，該群中心才會被選出來，可以藉著誤差的 大小來決定模糊因數的值，誤差越小表示重建回 來的集合越好。
濃縮了候選群中心的個數，再此利用 fuzzy 的概念從這些群中心裡面選出真正要的部份。利 用模糊因子演算法(FCM)的目標函數是
n1 i
n
1 j
m 2 ji cm
p
j i
f x -cc
J
(6)
如果在群數尚未確定的情況下，最小化這個 目標函數將無法減少任何群中心個數，也就是說，
假設原本有
n 個候選群中心，經由最小化 J
fcm之 後也會得到同樣n 個候選群中心，所以必須適當
的修改目標函數，在此增加群與群之間能量的參 數J
cc，2
1 1
i j i i j
cc cc cc
J
(7) 其中 β 是候選群中心的個數，所以新的目標函數
（Jnew）可以改寫成
1 1
n 2 1 i
n 1 j
m 2 ji p
i j i
j i j
i cc fcm new
cc cc -cc
x J J J
新目標函數的意義在於當多選到一個群中心 時，點對群中心的能量總和（Jfcm）將減少，而中 心對中心間的能量和（Jcc）將增加，所以當點對 群中心減少的能量和大於中心對中心增加的能量 和時，該群中心才會被選出來。
3.4.3 最佳化模糊因子演算法
群集的結果會因為模糊因數大小的不同而有 (8)
很大的差異，因此我們希望演算法能夠自動決定 一個好的模糊因數。(Chen & Lee, 2001)認為若是 針對遙測影像圖的話，模糊因數應該要設成 2.5 是最適當的。(Deer & Eklund, 2003)發現若是針對 純圖元（pure pixels）的影像，模糊因數應該設為 1.6，而如果是混合圖元（mixed pixels）的影像，
模糊因數應該設為 3。另外也有許多方法認為不 管針對任何資料的型態，模糊因數都應該設為 2。
因此在這裡採用(Okeke & Karnieli, 2006)的方法，
從一定範圍內中尋找一個最適當值。如以下簡述。
假設資料集
X = {x
1, x
2, …, x
k}且XRS，而X 是將 X 經過 FCM 之後所得到的群中心
* V* 和歸 屬 函 數 U* 的 乘 積 所 得 到 的 集 合 ，*
) ,..., , (
*1 *2 ** CS
c R
v v v
V
， 其 中 V*iRS ， 而cn
ik R
u
U* [*] 。假設藉著 FCM 執行的過程後所得 到 U 和 V，要重建或估算（predict）原本的資料 集 X，可以得到X*(U*,V*)。X 與
X 間一定會有誤差，
* 誤差的表示方式是：m X
X * , 0, (9) 所以可以藉著誤差的大小來決定模糊因數的值，
誤差越小表示重建回來的集合越好。
3.4.4 EM 群集演算法
EM(Expectation-Maximization) 利 用 第 一 階 段 的各個方法所得到的群數與群中心位置當作起始 條件，執行 E-Step 和 M-Step 兩步驟計算概似值 (likelihood)，當概似值的改變量小於某個門檻值 時，演算法就會停止執行。
兩個步驟如下式所表示：
E-Step (Expectation step)
K
k
k k i k
k k i k i
x f k P
x f k x P
k P
1
) ,
| ( ) (
) ,
| ( ) ) (
| (
(10)其中
K
k
k P k
P
1
1 ) ( 0 , 1 ) (
M-Step (Maximization step)
Ni
i N
i i i
k
x k P
x x k P
1 1
)
| (
)
| (
(11)
Ni
i N
i
k i i k
x k P
x x k P
1 1
2
)
| (
) )(
|
(
(12)
N
i P k xi
x N k P
1
)
| 1 (
)
|
( (13)
而概似值的表示方式如下：
N
i K k
i k
x f k P L
1 1
)
| ( ) (
log (14)
3.5 支撐向量機模式 (Support Vector Machine , SVM)
本研究所使用監督模式為支撐向量機模式 (SVM)，是機器學習的一種，也是一個新興的統 計學習系統，它屬於人工智慧的範疇，並且擁有 良好的學習與估算能力。
本研究所應用之監督模式並未自行學習建模，
而是採用林珈禎(2009)建立的學習資料庫來分類 及估算，其模式學習及驗證資料之採用拍攝時間 如表 2，共十五張衛星圖，分別為 98 筆學習樣本 數及50 筆驗證樣本數。
表2 學習及驗證數據組合(林珈禎，2009)
學習樣本數 驗證樣本數
98 筆 50 筆
拍攝時間 拍攝時間
1994/08/31、1995/01/06、
1996/12/05、1997/10/03、
1998/12/07、1999/08/10、
2000/03/15、2001/01/18、
2001/03/15、2002/12/02、
2003/05/07、2003/07/10、
2003/11/06、2004/01/06、
2004/09/07、
2004/10/14 2004/11/03 2005/08/02 2007/07/03 2008/03/04
3.6 空間統計分析
3.6.1 距離平方反比內插法 (Inverse Distance Weighted , IDW)
距離權重反比法係為地理資訊系統空間分析 外掛模組中提供之分析方法之一，為一傳統之空 間分析推估法，此方法針對每一個未知點的數值 推估，假設每一輸入點都有一個影響範圍，是利 用它鄰近的已知點之數值來進行加權內插運算，
而其中的未知點就是以與鄰近的已知點距離的遠 近來決定權重，隨著距離增加則影響減少，其概 念和物理學上的重力模式相似，公式如下：
N i
i N
i
i i
d w z d w y
x f
1 1
) ( )
( )
,
( (15)
其中，w(di
)是權重方程，z
i是第 i 個已知點的數
值，di是 i 點到未知點之間的距離。w(d
i) 的大
小是由 di的次方之反比來決定計算之，若採用
距離的高次方之反比為權重，則未知點受鄰近資 料點的影響增大。在這個內插過程中，關鍵在於 已知的樣本點來估算未知點的數值。3.6.2 共克利金內插法 (Co-kriging ,CK)
是將兩個或兩個以上具有空間相關性的區域 化隨機變數合併考慮進行空間資料推估，重點在 利用少數實測值輔助內插去推估未設測站位置資 料，來改善其推估誤差。其具有最佳線性不偏推 估(best linear unbiased edtimator)的特性。利用軟 體為Arc gis 9.2 的空間統計軟體。
線性：估計值為觀測值之最佳組合
mk k k
j n
j
j
u x v x
u
1 2 1 1
0
( ) ( )
(16) 不偏估計:估計值期望值與觀測值期望值相同)]
( [ ]
[u0 Eu xi
E (17)
將(16) 帶入(17)式得：
mk n k
j j
1 2
1
11 , 0
(18) (16)式最佳化：觀測值與預估值之差為最小] ) [(
]
var[ 0 0 0 0 2
2
0 u u E u u
將(16)式代入(18) 式，並利用 Langrange multiplier 及參數μ，則參數L 為：
n
i i i
k j m k
j n k
j
n i
m
j k i j
n i i
n
k i k i k
) ,u (u λ )
,v (v λ λ
) ,v (u λ
λ )
,u (u λ
λ L
0 1
2 2
2 1 2
1
n
1 j
k 2 1
2 0 0
cov 2
cov
cov 2
cov
-
1 2 ] ) u -
E[(u
mj λ2j
cov
(u0,vi)cov
(u0,u0)2
mk n k
j
j- μ λ
λ μ
1 2 1 1 2
1
1 2
2
將(19)式分別對 μ1、μ2及
1j、
k2取偏微分，並令其分式為零，得：
n
j k
n
j j
m
j j k j k
n
j j k j
m
j j k j k
n
j j k j
λ λ
) ,v (v μ
) ,v (v λ ) ,u (u λ
) ,v (u μ
) ,v (u λ ) ,u (u λ
1 1
1 2 1 0
1 1
1 2 1 0
1 1
0
; 1
cov cov
cov
cov cov
cov
以矩陣形式表現為
0
1
1 1 1 1 1 0 0 0
0 0 0 0 0 1 1 1
0 1
0 1
0 1
0 1
0 1
0 1
0 0 2
0 1
0 0 2
0 1
2 1 2
1
2 2 2 1 2 2 2 2 2 1
1 2 1 1 1 1 1 2 1 1
2 1 2
1
2 2 2 1 2 2 2 2 1 2
1 2 1 1 1 1 2 1 1 1
2 1 2 2 2 1 2 1 2 1 1 1
u v
u v
u v
u u
u u
u u
m n
v v v v v v v u v u v u
v v v v v v v u v u v u
v v v v v v v u v u v u
v u u u v u u u u u u u
v u v u v u u u u u u u
v u v u v u u u u u u u
m n
m m m
m m n m m
m n
m n
m n n
n n n n n
m n
m n
C C C C C C
C C
C C C C
C C
C C C C
C C
C C C C
C C
C C C C
C C
C C C C
C C
C C C C
得
1j、
1k後，代入(16 )式，得估計值Z0。即為 Ordinary Cokriging System 估計變異數為：(19)
(20)
(21)
n
j
j j
E Covu u Covu u
1
0 1
0 0
2 ( , ) ( , )
1 1
0
2 ( , )
m k
k
kCovv u
估計誤差為SE
E2 (23) 本 研 究 所 使 用 的 方 法 為 一 般 共 克 利 金 (ordinary cokriging) 。 以 環 保 署 水 質 採 樣 站 (sampling stations)實測計算得之河川水質污染指 標值，結合與水質採樣同一日期之非監督模式所 估算衛星圖空間水質污染面化值，利用空間統計 之 Co-kriging 進行內插，藉此方式估算衛星圖上 任一未設測站之水質空間污染情形。CK 法以抽 點之交叉驗證法進行模式評估，每次刪除一筆實 測資料，以剩餘所有資料進行內插，直至每一測 站皆被刪除驗證過，得以客觀評估空間估算之可 靠性。3.6.3 研究模式建立
本研究利用非監督分類遙測資料及監督學習 估算遙測資料兩種方式，再分別加入當日少數定 點水質採樣資料作共克利金內插(本文後，非監督 結合共克利金法簡稱 UCK；監督結合共克利金法 簡稱SCK)，但鑑於非監督分類衛星河道圖需耗費 相當長時間，將高程(DEM)移至共克利金內插時 作趨勢消除變量，藉以減少非監督分類時間(本文 後 ， 非 監 督 結 合 共 克 利 金 法 加 入 DEM 簡 稱 UCKD；監督結合共克利金法加入 DEM 簡稱 SCKD)，並比較 UCKD 及 SCKD 兩模式總運算時 間及精確度之優劣。 最後再與傳統僅用監測站值 作空間內插(本文後簡稱 IDW 法)結果作比較，茲 將各模式簡稱對照綜合列表3 如下。
表3 本研究使用模式之對照簡稱
模式 簡稱
非監督共克利金法 UCK 監督共克利金法 SCK 非監督共克利金加入高程 UCKD
監督共克利金加入高程 SCKD 共克利金法 CK 傳統內插法 IDW
4. 結果與討論
4.1 模式建立
本模式研究之主軸分為兩部份：第 1. 以非監 督分類遙測圖結合水質採樣點共克利金法(UCK) 來取代監督式估算遙測圖結合水質採樣點共克利 金法(SCK)。主要利用非監督模式自動最佳化分 類數的優點，接著以 CK 法利用少數實測值輔助 內插未設測站位置的估算值，每次抽取一個測站 資料，以其他測站內插對其進行交叉推估驗證，
直到每個測站皆被抽取完為止，最後比較非監督 共克利金法(UCK)與監督共克利金法(SCK)之估算 精確度。第 2. 為縮短非監督模式分類衛星河道圖 時間。嘗試減少變量藉以縮短非監督模式分類時 間，並將 DEM 變量移至共克利金內插時作趨勢 消除變量，以考慮水質因高度變化的趨勢。
本 研 究 參 照 林 珈 禎 (2009) 文 獻 中 以 STATISTICA 8.0 中之變異數分析 F 檢定之變量篩 選結果，選用 B3/B2、NDVI、DEM 來作為非監 督及監督模式之第一種輸入變量。然為了有效縮 短非監督模式分類衛星河道圖的時間，便將數值 高程(DEM)資料移至共克利金法作趨勢消除之內 插變量，因此第二、三種變量組合分別為 B3/B2、
NDVI 及 B3/B2，藉以比較未減少變量前及減少 變量後之估算精確率及分類時間。最後再加入傳 統內插法(IDW)，僅以環保署 RPI 實測值作平方 反比內插，並與 UCKD 法之估算結果做比較。本 次研究所使用電腦之硬體是 Intel 四核心 CPU Core 2 Quad Q6600 ， 時 脈 為 2.4G ， 記 憶 體 (RAM)2G。
4.2 估算結果
4.2.1 衍生變量組合
(B3/B2、NDVI、DEM)
本研究先選用 B3/B2、NDVI、DEM 作為非 監督模式之輸入變量，來分類十張河道圖，分為 七張淡水河 (2004/11/03、2005/08/02、2007/07/03、
(22)
2008/03/04、2009/02/09、2009/05/08、2009/05/29) 及 三 張 高 屏 溪 (2004/10/14 、 2009/01/03 、 2010/01/02)，共 95 筆測站資料。
本研究利用非監督及監督模式進行衛星圖分 類及估算後，再以共克利金內插法進行抽點交叉 驗證，藉此方式估算衛星遙測影像河川水體的污 染情形，表4 為 UCK 及 SCK 法綜合估算結果。
UCK 法之估算結果總平均 RMSE 為 1.0，MAPE 為 28%，R2為 0.76；SCK 法之估算結果總平均 RMSE 為 1.3，MAPE 為 33%，R2為0.64。
UCK 法中 2010/01/02 高屏溪流域估算結果之 RMSE 為 0.77 低於 SCK 法估算之 RMSE 值 0.84，
但 SCK 法 MAPE 卻僅高出 UCK 法 3%，顯示 UCK 法之估算誤差皆在上游監測站 RPI 值較小處，
因此少量誤差就會造成APE 較高之情形，而 SCK 法之估算誤差則集中在中下游 RPI 值較大之監測 站，如旗南、嶺口及里港橋，形成 MAPE 僅高出 3%之情況，綜合比較後可發現 UCK 法之十張河 道圖的估算結果 RMSE 及 MAPE 值皆比 SCK 法 稍佳。
總運算時間共分為衛星影像前處理、建立最 佳化模式、模式最佳化分類及估算、共克利金抽 點交叉驗證四部份探討，其結果顯示於表 5。 在
衛星影像前處理部分，兩模式皆須先將衛星圖之 河道擷取出來，而非監督模式僅需做單張衛星圖 感測器校正即可自行最佳化分類數，不需向監督 模式為了使學習樣本一致性，需同時進行多張的 校正，最後本研究僅進行感測器錯誤之修正，並 未進行大氣校正，而監督模式必需建立學習資料 庫，包含學習圖資之採樣及河道處理時間，因此 衛星影像前處理時間比非監督多了 1 小時；建立 最佳化模式部分，由於非監督模式可自行最佳化 分類，故不需花費任何時間，而監督模式學習前，
必需採用最佳化參數(C 值、γ 值、ε 值)之網格搜 尋法，來找尋最佳參數，其過程十分繁瑣且複雜 性高，故實際花費 3 小時建立最佳化模式時間，
因此衛星影像前處理及建立最佳化模式部分非監 督模式具有簡單、快速等優點；模式最佳化分類 及估算方面，非監督雖可自行最佳分類，但相對 花費的演算時間也較長，而優點是不用像監督模 式般需要大量且具代表性的學習樣本，只需要利 用少部分現場水質污染監測分類值輔以衛星遙測 圖之光譜分類值來當作最終水質污染判斷的參考；
另外兩模式之共克利金抽點交叉驗證所花費時間 是相同的，故影響總時間不大。
表4 UCK 法及 SCK 法之估算結果(變量：B3/B2、NDVI、DEM)
變量組合 (B3/B2、NDVI、DEM)
模式 UCK SCK
RMSE MAPE(%) R2 RMSE MAPE(%) R2
2004/10/14 1.2 69 0 1.3 62 0
2004/11/03 1.0 24 0.69 1.1 26 0.68
2005/08/02 1.6 30 0.60 1.9 38 0.45
2007/07/03 1.3 30 0.71 1.3 31 0.71
2008/03/04 0.8 18 0.81 0.8 19 0.83
2009/01/03 0.97 30 0.31 0.99 33 0.27
2009/02/09 0.67 16 0.89 0.7 17 0.88 2009/05/08 0.89 20 0.76 0.9 20 0.75
2009/05/29 0.65 17 0.88 0.66 17 0.88
2010/01/02 0.77 20 0.24 0.84 23 0.17
平均 1.0 28 0.76 1.3 33 0.64
表5 平均單張衛星圖之處理時間(三個變量組合) 衛星影像
前處理
建立最佳 化模式
模式最佳化分類及估算 (B3/B2、NDVI、DEM)
CK 法抽點 交叉驗證
總運算 時間
非監督 8 hr 0 hr 3 hr
50 min
約11 hr 50 min
監督 9 hr 3 hr 30 min 約13 hr
20 min * 監督模式之前處理時間包括建立學習資料筆數時間
表6 UCKD 及 SCKD 之估算結果(變量：B3/B2、NDVI)
變量組合 (B3/B2、NDVI)
模式 UCKD SCKD
RMSE MAPE(%) R2 RMSE MAPE(%) R2
2004/10/14 1.2 67 0 1.3 69 0
2004/11/03 1.1 24 0.69 1.0 24 0.7
2005/08/02 1.6 30 0.60 1.8 37 0.47
2007/07/03 1.2 29 0.74 1.3 31 0.71
2008/03/04 0.8 18 0.84 0.83 19 0.83
2009/01/03 0.97 34 0.31 0.98 33 0.31
2009/02/09 0.67 18 0.89 0.69 18 0.88
2009/05/08 0.87 20 0.76 0.92 21 0.74
2009/05/29 0.64 18 0.89 0.70 19 0.86
2010/01/02 0.75 19 0.24 0.8 22 0.22
平均 1.0 28 0.77 1.1 30 0.75
* 高程(DEM)移至共克利金內插時作趨勢消除變量 綜合上述以 B3/B2、NDVI、DEM 作為輸入 變量之模式分類時間來看，非監督模式分類單張 衛星圖大約需要3 小時，總運算時間達 12 小時，
而監督模式大約需要 30 分，總運算時間約為 13 小時 20 分，約比非監督多了 1 小時 30 分，但非 監督模式最佳化分類時間仍稍嫌過長，因此為了 有效縮短非監督模式的分類時間，本研究試圖將 高程(DEM)移至共克利金內插時作趨勢消除變量，
共同進行空間資料推估，使非監督模式之輸入變 量減少為 B3/B2、NDVI，且監督模式也輸入相同 變量，再比較UCKD 及 SCKD 兩模式之估算精確 率及估算時間。
4.2.2 衍生變量組合(B3/B2、NDVI)
本研究嘗試減少非監督模式之輸入變量，將 高程(DEM)移至共克利金內插時作趨勢消除變量，
考慮水質因高度變化的趨勢，共同進行空間資料 內插，試圖減少非監督模式的分類時間，以提高 模式效率。本節模式之輸入變量為 B3/B2、NDVI，
分類 10 張衛星圖，其綜合結果如表 6，總運算時 間為表7。UCKD 法估算結果總平均 RMSE 為 1.0，
MAPE 為 28%，R2為0.77；SCKD 法估算結果總 平均RMSE 為 1.3，MAPE 為 36%，R2為0.75。
結果顯示UCK 法未將高程(DEM)移至共克利 金法前之估算結果總平均 RMSE 為 1.0，MAPE 為 28%，R2為 0.76，而將高程(DEM)移至共克利 金法作趨勢消除之內插變量後並未造成估算誤差 提升，比較後結果差異不大，另外 SCKD 法減少 高程(DEM)後其誤差比 UCKD 法來得稍差，顯示 減少變量後UCKD 法估算準確率比 SCKD 法稍佳。
綜合上述比較，兩模式減少變量後仍具有相當高 的準確率，但非監督模式的分類時間則大幅降低，
原三變量組合(B3/B2、NDVI、DEM)的平均估算 時間要 3 小時，而減少後的兩變量組合(B3/B2、
NDVI)只需要 1 小時 30 分，整整縮短 1 個半小時 的模式分類時間，但將高程(DEM)移至 CK 法作 趨勢消除變量，雖造成抽點交叉驗證時間小幅增 加 20 分鐘，但非監督總運算時間只需約 10 小時 40 分，已比監督模式少了 3 小時，最後顯示將 DEM 移至共克利金法所增加時間對整體影響不大，
且減少變量後非監督模式分類衛星圖之估算效率 更佳。
4.2.3 衍生變量組合(B3/B2)
本研究試圖再縮短非監督模式分類衛星河道 圖的時間，將 NDVI 捨去僅輸入 B3/B2 之單一變 量，分類 10 張衛星圖(2004/10/14、2004/11/03、
2005/08/02、2007/07/03、2008/03/04、2009/02/09、
2009/05/08、2009/05/29、2009/01/03、2010/01/02)，
同樣將高程(DEM)移至 CK 法作趨勢消除變量，
探討其估算精確率及模式估算時間，其綜合結果 如表8。UCKD 法估算結果總平均 RMSE 為 1.0，
MAPE 為 28%，R2為0.77；SCKD 法估算結果總 平均RMSE 為 1.1，MAPE 為 30%，R2為0.74。
兩模式估算誤差方面，單一變量(B3/B2)之 UCKD 法平均估算結果 RMSE 與兩個變量組合 (B3/B2、NDVI)之平均估算結果相同，且 SCKD 法不論單一或兩變量組合其平均估算結果亦皆相 同，RMSE 值為 1.1，MAPE 值則同樣 30%，而仔 細比較兩模式後可發現單一變量與兩變量組合之 估算結果差異不大，皆具有相當穩定的準確率。
SCKD 法之估算誤差 RMSE 值僅大於 UCKD 法 0.1，MAPE 值則大為 2%，總結來說 UCKD 法估 算結果僅比SCKD 法稍佳。
表 9 結果顯示僅輸入單一變量(B3/B2)，非監 督模式分類時間更為省時，平均模擬單張衛星圖 時間只要30 分鐘即可分類完成，比起兩個變量組 合(B3/B2、NDVI) 1 小時 30 分的分類時間，約縮 短 3 倍，表示非監督模式分類衛星圖資所需時間 已與監督模式相同，加上非監督模式衛星圖前處
理及建立最佳化模式簡單、快速，無監督模式學 習建模及參數網格搜尋之耗時及複雜性，顯示非 監督模式分類程序簡潔之效率。
4.2.4 距離平方反比內插法 (Inverse Distance Weighted, IDW)
由於目前環保署所公佈的河川流域空間面化 圖皆是以單純的 RPI 值進行內插產生，未考慮其 他變因資料，雖然過程快速，但產生的判斷誤差 也相對較大。有鑑於此，本研究另將10 張衛星圖 僅以環保署監測站值之河川水質污染指標(RPI)，
以距離平方反比權重法(IDW)進行內插，藉此方 式估算各張衛星圖之水質空間污染情形，並與 UCKD 法(變量：B3/B2、NDVI 及 B3/B2)進行比 較，其綜合結果如表10。
表10 中 IDW 法總平均 RMSE 為 1.2，MAPE 為32%，R2 為 0.72；而 UCKD 法(二個變量至單 一變量)之總平均 RMSE 則位於 1.0，MAPE 為 28%，R2 為 0.76-0.77，在 MAPE 值的部分，兩 者所得知結果大致上差不多，而 RMSE 值的結果 IDW 法卻高出 UCKD 法許多，這表示 IDW 法估 算 誤 差 皆 集 中 在 RPI 範 圍 值 較 大 的 部 分 ， 而 UCKD 法估算誤差集中在 RPI 範圍值較小的部分，
因此 UCKD 少量的誤差變會造成 MAPE 較高的 情況，總結來說 IDW 法之估算誤差仍比 UCKD 法來得高。以 IDW 法內插估算淡水河(2009/02/09、
2009/05/08) 之 RMSE 值 0.85 及 1.03，而高屏溪 (2009/01/03、2010/01/02) 之 RMSE 值為 1.01 及 0.84，此四張衛星圖估算誤差皆較高，其中較大 誤差集中在淡水河之秀朗、華江及中正橋與高屏 溪之嶺口及里港橋等中上游測站，主要應該是傳 統內插IDW 法僅考慮未知測站與鄰近已知測站間 平面距離為內插依據，故容易造成估算偏差。
表7 平均單張衛星圖之處理時間(兩個變量組合) 衛星影像
前處理
建立最佳 化模式
模式最佳化分類及估算 (B3/B2、NDVI)
CK 法抽點 交叉驗證
總運算 時間 非監督 8 hr 0 hr 1 hr 30 min
1 hr 10 min
約10 hr 40 min
監督 9 hr 3 hr 30 min 約13 hr
40 min * 監督模式之前處理時間包括建立學習資料筆數時間
表8 UCKD 及 SCKD 之估算結果(變量：B3/B2)
變量組合 (B3/B2)
模式 UCKD SCKD
RMSE MAPE(%) R2 RMSE MAPE(%) R2
2004/10/14 1.2 67 0 1.2 68 0
2004/11/03 1.0 23 0.70 1.1 24 0.69
2005/08/02 1.5 30 0.61 1.8 38 0.46
2007/07/03 1.3 31 0.71 1.3 31 0.70
2008/03/04 0.9 17 0.79 0.8 19 0.81
2009/01/03 0.98 35 0.31 0.97 33 0.28
2009/02/09 0.7 18 0.88 0.70 19 0.88
2009/05/08 0.87 21 0.77 0.94 23 0.72
2009/05/29 0.66 18 0.88 0.72 19 0.86
2010/01/02 0.78 21 0.22 0.84 23 0.14
平均 1.0 28 0.77 1.1 30 0.74
* 高程(DEM)移至共克利金內插時作趨勢消除變量
4.3 河川流域水質污染面化圖 (RPIM)
因單一變量(B3/B2)較為省時且模式效能評估 結果較佳，故本研究利用非監督及監督模式輸入 單一變量進行估算，皆以 CK 法並加入傳統內插 IDW 法建立 RPIM，最後選用 2009/05/08 三張不 同模式之面化圖作比較，如圖 3 至 5 所示，進一 步比較其成效，探討此兩方法實際運用之可行性。
圖3 為 UCKD 法之面化結果，圖 4 為 SCKD 法之面化結果，圖5 為 IDW 法之面化結果，其中 UCKD 法及 SCKD 法中淡水河之華中橋實測值為 7，但面化的分類結果介於 6.1-7 之間，有些許的 誤 差 ， 另 外 上 游 的 部 份 SCKD 法 分 類 藍 色 較 UCKD 法稍多，但總結而言，UCKD 法及 SCKD 法兩模式皆有加入衛星光譜值作調整，因此兩模
式建立之河川流域水質污染面化圖無大太差異。
接著以傳統內插IDW 法(圖 5)與 UCKD 法作 比較，可清楚發現在有監測站附近的內插值基本 上差異不大，但在右上角基隆河支流，靠近關渡 橋的污染流段，其判別為黃色(RPI 值介於 4.1-5.0)，
另外於左下方支流則判別為紅色(RPI 值介於 4.1- 5.0)，這主要原因為 IDW 法僅依靠鄰近監測站資 訊作平方內插，而在無監測站支流亦僅能仰賴最 近測站作水體污染分類，但這可能存在有嚴重的 錯誤，因不同支流污染物也不盡相同，且支流不 應該受到主流所影響，在圖3 以 UCKD 法面化圖 中，可發現基隆河支流及左下方無測站支流，由 於加入了衛星光譜分類值作內插，因此可分類出 符合該支流之污染顏色，可避免無測站支流受主 流或其他分流所影響，充分顯示以 UCKD 法建立 河川水質污染面化圖之優點。
表9 平均單張衛星圖之處理時間(單一變量組合) 衛星影像
前處理
建立最佳 化模式
模式最佳化分類及估算 (B3/B2)
CK 法抽點 交叉驗證
總運算 時間 非監督 8 hr 0 hr 30 min
1 hr 10 min
約9 hr 40 min
監督 9 hr 3 hr 30 min 約13 hr
40 min * 監督模式之前處理時間包括建立學習資料筆數時間
表10 IDW 法與 UCKD 法之總估算結果
變量組合 無 B3/B2、NDVI B3/B2
模式 IDW UCKD UCKD
RMSE MAPE(%) R2 RMSE MAPE(%) R2 RMSE MAPE(%) R2 2004/10/14
1.5 77 0 1.2 67 0 1.2 67 0 2004/11/03
1.1 20 0.69 1.1 24 0.69 1.0 23 0.70 2005/08/02
1.6 36 0.52 1.6 30 0.60 1.5 30 0.61 2007/07/03
1.5 39 0.59 1.2 29 0.74 1.3 31 0.71 2008/03/04
0.9 22 0.8 0.8 18 0.84 0.9 17 0.79 2009/01/03
1.01 35 0.14 0.97 34 0.31 0.98 35 0.31 2009/02/09
0.85 21 0.84 0.67 18 0.89 0.7 18 0.88 2009/05/08
1.03 22 0.69 0.87 20 0.76 0.87 21 0.77 2009/05/29
0.79 18 0.84 0.64 18 0.89 0.66 18 0.88 2010/01/02
0.84 23 0.24 0.75 19 0.24 0.78 21 0.22
平均 1.2 32 0.72 1.0 28 0.77 1.0 28 0.77
圖3 淡水河流域水質污染面化圖—UCKD 法(2009/05/08) 關渡橋4.25
重陽橋5
忠孝橋6.25 華江橋6.25
華中橋7
中正橋5.5
秀朗橋2.75
直潭橋1 屈尺堰1
■：0.0 – 1.0
■：1.1 – 2.0
■：2.1 – 3.0
■：3.1 – 4.0
■：4.1 – 5.0
■：5.1 – 6.0
■：6.1 – 7.0
■：7.1 – 8.0
■： ＞ 8.1
