行政院國家科學委員會專題研究計畫成果報告:地理資訊系統SuperGIS 之地理統計演算元件平台技術開發

行政院國家科學委員會專題研究計畫 成果報告

地理資訊系統 SuperGIS 之地理統計演算元件平台技術開發

計畫類別： 個別型計畫 計畫編號： NSC93-2415-H-002-071-CC3 執行期間： 93 年 11 月 01 日至 94 年 10 月 31 日 執行單位： 國立臺灣大學生物環境系統工程學系暨研究所 計畫主持人： 林裕彬 共同主持人： 鄭克聲，張尊國 計畫參與人員： 吳振發,史瓊雯 報告類型： 精簡報告 處理方式： 本計畫為提升產業技術及人才培育研究計畫，不提供公開查詢

中 華 民 國 95 年 1 月 24 日

行政院國家科學委員會補助產學合作研究計畫成果完整報告

地理資訊系統 SuperGIS 之地理統計

演算元件平台技術開發

計畫類別：產學合作研究計畫

計畫編號：NSC 93－2415－H－034－008－CC3

執行期間：93 年 8 月 1 日至 94 年 7 月 31 日

計畫主持人：林裕彬

、王能超

共同主持人：

張尊國、鄭克聲

計畫參與人員：

吳振發、周黛君、林啟陽、史瓊雯、洪萬程

處理方式：完整報告內容因涉及專利、技術移轉案或其他智慧財產

權，不予公開。

執行單位：

中國文化大學景觀系、崧旭資訊股份有限公司

中 華

民 國 九 十

四 年 十 月

三 十 一 日

中文摘要

地理統計目前已廣泛應用於各領域中，包括社會與自然科學，其是以區域化 變數理論為基礎，利用變異圖進行一區域化隨機變量之資料結構分析，而得到區 域化變數空間相依性，以克利金等方法推估未採樣點數值之最佳線性無偏估計； 若能再將提供空間資料儲存、展示、處理及運算的地理資訊系統與地理統計方法 相互結合，不但能兼具地理資訊系統探討空間資料與屬性資料的查詢、展示與基 本統計能力外，更能增進地理資訊系統之空間資訊分析與估算功能，了解空間資 料的空間變異性與相依性，並進行空間推估與模擬，進而提昇地理資訊系統輔助 空間資料決策分析的能力。因此，地理統計演算元件對於地理資訊系統的相關應 用甚為重要，且不可或缺的有效工具。 本研究擬首先回顧與整理地理統計之相關理論與計算方法，並篩選出最適合 且普遍應用之方法以做為構建地理統計演算元件之基本方法；此外，亦設計各方 法之測試例題進行試驗，以評估開發之地理統計元件演算法之求解績效。本研究 計畫之具體成果除可建置一套屬於國人自行研發的架構在 GIS 上的地理統計演 算元件平台技術，提升相關產品之功能與競爭力外，並可為產學雙方提供一個很 好的合作模式，增進學術理論之實務應用。 關鍵詞：地理統計、地理資訊系統、空間結構、地理統計運算方法、地理統計運 算元件

Abstract

Geostatistics, based on regionalized variable theorem incorporating the spatial or temporal characteristics of actual data into statistical estimation processes to do linear unbiased estimation or simulation, is widely applied in many fields included social and natural study. Combining geostatistical techniques and geographic information system may not only provide the basic functions of geographic information system, but also improve the functions of spatial information analysis and estimation to realize spatial variations, dependences, estimation and simulation of spatial data, especially to improve the ability of decision making support of geographic information system. Therefore, in this project, we first review the geostatistical methods and algorithms to select suitable, useful and applicable geostatistical methods which are the bases of the geostatistical components of the geographic information system SuperGIS. A number of test examples also set up to evaluate the performances of the components comparing to most popular geostatistical software. SuperGeo Technologies Inc. performed the coding the algorithm components, user interface and functional modules. The performances of these codes reported and feedback to designer in the SuperGeo Inc.. Final, this team produces high performance geostaistical components for SuperGoe which is a native Geographic information system in Taiwan.

Keywords: Geostatistics, Geographic Information System, Spatial Structure, Geostatistical methods, Geostatistical algorithm components

目 錄

中文摘要………I Abstract………..II 目錄………..III 一、前言………..1 二、研究目的………..1 三、研究方法………10 四、結果與討論………31 五、結論與建議………39 六、參考文獻………....40 計畫成果自評………..45

一、前言

地理資訊系統（Geographic Information Systems）為一套整合各項相關地理 資訊作業，架構於一完整豐富的地理資料庫之上，並具有資料擷取、編修、更新、 儲存、查詢、處理、分析及展示等功能的軟體工具。目前已普遍運用於自然科學 與社會科學中。隨著地理資訊系統的重要性與普及性增加，國內各專家業者對於 地理資訊系統的需求也日益增多，當前國內市場佔有率最高的 GIS 軟體為國外 所開發之 Arcview 與 Mapinfo 兩套軟體。但由於國外軟體開發者的時空背景與開 發原意不同，因而在國內時常發生使用上的困擾與尋求技術支援不足的問題。因 此，崧旭資訊秉持多年執行專案的經驗中，有感於地理資訊系統在國內所面臨的 問題，是必須國內業者自行從軟體核心開始解決的問題，也希望國內 GIS 專業 精英人才與技術能根留台灣，因此於 2002 年開始積極專研於國內本土 GIS 市場 的需求調查與 GIS 軟體的研發工作，於 2003 年開始推出一系列的 GIS 軟體，於 2004 年正式推出 SuperGIS 桌上型地理資訊系統。隨著桌上型地理資訊系統的產 出，各領域之 GIS 專才開始強列希望與期許崧旭資訊在資料的分析上能有更強 大與實用的分析功能，也希望能將當前具有高度空間分析與推估能力的地理統計 方法加入地理資訊系統中，使國人自行研發產出的地理資訊系統能更加的強大， 與國外軟體並駕齊驅，提昇產業競爭力。而相關專才參與，讓本土的地理資訊系 統的發展能更加不受侷限與多元化，遂積極尋求國內大專院系及研究機構合作。 二、研究目的 本研究案之目的即由中國文化大學景觀學系與崧旭資訊進行產學合作計 畫，共同合作開發 SuperGIS 之地理統計演算元件平台技術開發，並由台灣大學 生物環境系統工程學系提供相關理論資詢，以提昇產品的分析功能與協助相關領 域決策之運用性，加強產業的競爭力。 地理統計目前已廣泛應用於各領域中，包括社會與自然科學，其是以區域化 變數理論為基礎，利用變異圖進行一區域化隨機變量之資料結構分析，而得到區 域化變數空間相依性，以克利金等方法推估未採樣點數值之最佳線性無偏估計 （Journel 與 Huijbregts，1978；Cressie，1993；Isaaks 與 Journel，2002）；若能 再將提供空間資料儲存、展示、處理及運算的地理資訊系統與地理統計方法相互 結合，不但能兼具地理資訊系統探討空間資料與屬性資料的查詢、展示與基本統 計能力外，更能增進地理資訊系統之空間資訊分析與估算功能，了解空間資料的 空間變異性與相依性，並進行空間推估與模擬，進而提昇地理資訊系統輔助空間 資料決策分析的能力。因此，地理統計演算元件對於地理資訊系統的相關應用甚 為重要，且不可或缺的有效工具。 當前地裡統計的發展已漸趨成熟，並已運用於相當多的層面，包含有土壤(鄭

森源與萬鑫森，1994；Chang 等人，1999、2000；Lin 與 Chang,，2000；Lin 等 人, 1995, 2000, 2001；林裕彬等人，1998、1999、2001、2002；Lin，2002；林裕 彬與張尊國，1999；徐貴新等人，1999；張尊國等人，2001；莊凱瑋等人，1996、 1999)、地質（Jaquent 與 Carniel，2001, 2003）、水文（Lin, 1995；詹益璋等人， 2001、2002；Tung 等人 2003；Lin 等人，2004；Lin，2004）、農業(Diodato 與 Ceccarelli, 2004；林裕彬與譚義績，1999；林允斌等人，2004)、氣象(Dalezios 等人，2002；Nanus 等人，2003；Rubel 與 Hantel，2001；Lin 與 Rouhani，2001)、 生態(Franklin 與 Mills, 2003; Wylie 等人, 2002；Lin 等人，2002；林裕彬等人， 2001；)、森林(Nanos 與 Montero, 2002；Clark 等人，2004；Hudak 等人，2002)、 環境規劃與工程（Cheng 等人，2004；Uphoff 等人，2004；Ozturk 與 Nasuf，2002； Lin 等人，2004）等。搜尋目前國內各大學最普遍電子期刊全文查詢系統 SDOS (www.sdos.ejoural.ascc.net) 與 SwetsWise (http://www.swetswise.com/eAccess /searchArticles.do) 系統，分別以 kriging 與 geostatistic 為查詢關鍵字，查詢的範 圍分別為文章名稱(title)、摘要(abstract)、所有欄位(all fields)，已初步瞭解目前 地理統計學目前於國外期刊發展情形（詳見表 1）。於 SwetsWise 系統中所查詢 到篇數較多，以 kriging 為關鍵字查詢 1995 至 2004 年所有欄位，計有 388 篇， 若以 geostatistics 為關鍵字，查詢到 421 篇，顯示地理統計學相關研究蓬勃發展。 表 1 地理統計學於電子全文查詢系統查詢統計表 全文檢索 資料庫 資料年期 Kriging (查詢 title) Kriging (查詢 abstract) Kriging (查詢 all field) geostatistics (查詢 title) geostatistics (查詢 abstract) geostatistics (查詢 all fields) 1995-2004 65 323 369 47 106 324 SDOS 2001-2004 24 125 146 13 37 139 1995-2004 151 - 388 114 - 421 SwetsWise 2001-2004 103 - 278 63 - 276 資料來源：本研究整理。 除了數量上的分析外，在應用領域的分析上，國外部分，以 Kriging 為關鍵 字，查詢 SDOS 系統 abstract 欄位之結果，計有 125 篇，逐一分析結果發現應用 於土壤、水文相關領域的研究數量最多，其次為氣候及方法論之研究，並且拓展 至植物、動物、生態、地質、環境條件、環境規劃、農業等，甚至是機械工程、 電腦科學與醫學等（詳見表 2）。國內部分，同樣以 kriging 為關鍵字分別查詢 國家圖書館中所有國內期刊與博碩士論文，並分析應用於各領域狀況，結果發現 國家圖書館收錄發表於國內期刊者計有 52 篇，博碩士論文有 56 篇，主要集中於

的情況較少，詳細統計內容參見表 2。 表 2 國內外地理統計應用領域分析表 SDOS (2001~2004) 國內期刊 (1994~2004) 國內博碩士論文 (1992~2004) 領域別 篇數 篇數 篇數 土壤 _{31 18 19} 水文 24 18 27 方法論、作業研究 12 2 0 氣候 12 6 2 地理 8 0 0 植物與森林 7 2 0 地質 5 0 6 動物 4 0 0 生態 4 0 0 環境品質 4 3 0 環境規劃 4 1 0 農業 3 0 0 機械工程 _{2 0 0} 電腦科學 2 0 0 醫學 2 0 0 其他 1 2 2 合計 125 52 56 資料來源：本研究整理。

而地理統計軟體亦相當多，最普遍的有 GS+（Deutsh 與 Journel, 1992）、GSLIB

（Gamma Design Software，1995）、VARIOWIN（Pannatier，1996）等，以及外 掛於 Arcview（ESRI，1998）或 ArcGIS 的地理統計模組，及其他相關軟體，彼 此間優劣差異比較如下：

表 3 地理統計相關軟體分析表

主要功能 軟體名稱 研發期程 作者 系統平台 原始碼

基本統計 估計 模擬 其他

價格

3DField 2.4.5 1999-2003 Vladimir Galouchko MS Windows

95/98/2000/XP

未公佈 ⊙ ⊙ ⊙ 免費

3plot98 4.60 1997-1999 M. Kanevski, S. Chernov, V.

Demyanov

MS Windows95

or Windows 3.x

with Win32s

未公佈 ⊙ ⊙ ⊙ 免費

AGROMET 1.0. 1995 By P. Bogaert, P. Mahau and

F. Beckers

MS-DOS, UNIX C++ ⊙ 免費

BMElib 2001~ George Christakos, Marc

Serre, Patrick Bogaert

Windows &

UNIX

MATLA

B

⊙ ⊙ ⊙ ⊙ BMElib 免費

COSIM 1986 MIT DOS FORTRA

N

⊙ 免費

DACE 2.5 2002-2003 Hans Bruun Nielsen,

Lophaven and Jacob

MATLAB Matlab ⊙ ⊙ 免費

EasyKrig 2.1. 2001 Dezhang Chu MATLAB Matlab ⊙ ⊙ 免費

ECOSSE 2001 Isobel Clark, Noel Cressie,

William Harper

Windows

95/98/2000/NT

未公布 ⊙ ⊙ ⊙ ⊙ $US1,000

ExploStat 2000-2003 Lodewijk Hazelhoff,

Almere

MS Windows 未公布 ⊙ ⊙ 免費

GeoPack 1990 S.R. Yates(1) and M.V. Yates

(2)

DOS 未公布 ⊙ ⊙ 200 US$

GeoR 1.3.11 1997-2003 Paulo Justiniano Ribeiro Jr.,

Peter J. Diggle

Linux, Windows,

UNIX

R ⊙ ⊙ ⊙ ⊙ 免費

GSII 1997 Geostat Systems International

Inc. Windows (95, NT) 未公布 ⊙ 免費 Geostatistical Toolbox 1.30.

1990 Roland Froidevaux DOS 未公布 ⊙ 免費

Geostatistics

Package from Cyze

and Associates

2002 Cyze & Associates, Ltd. Windows 未公布 提供與 GSLIB 及

SAGE

2001 搭 免費

Geostatistics with

SAS

1999-2000 SAS/STAT OS/2, UNIX,

OpenVMS Alpha,

Windows

未公布 ⊙ 未知

Geostatistis

Analyzer

2003 4MationGeo, Inc Windows，需要

Statistica 6.0 或

MapInfo

Professional 7.0/

未公布 ⊙ 495 US$.

Geostokos Toolkit 1997～ Isobel Clark Windows

95/98/2000/NT

未公布 ⊙ ⊙ ⊙ ⊙ 收費

GMS 3.0 1997 Engineering Computer

Graphics Laboratory of

Brigham Young University

UNIX, Windows 未公布 ⊙ ⊙ ⊙ US$ 5,095 (Unix) or US$

3,395

(Windows)

GRIDSTAT 1995-2001 Applied Computer

Engineering

Windows (3.1.,

95, NT), UNIX

未公布 ⊙ ⊙ ⊙ 495 US$

GS+ 7 1989 - 2004 Gamma design Software Windows

XP/2000//NT

未公布 ⊙ ⊙ ⊙ ⊙ 1000 US$、199 US$（學生用）

gstat 2.3.8 1996-2002 Edzer J. Pebesma (and

others)

Win32 (NT/9x)、

Linux

未公布 ⊙ ⊙ ⊙ 免費

GsTL 2001-2002 Nicolas H. Remy 支援 ANSI C++的 編譯程序

C++ ⊙ ⊙ 免費

GViz 1997 Correlations Company Windows (95,

NT)

未公布 ⊙ ⊙ US $ 15,000

ISATIS 4.1 1997-2002 Geovariances et Ecole des

Mines de Paris

UNIX,

WinNT/2000/XP

未公布 ⊙ ⊙ ⊙ ⊙ 未知

Isomap 1999 GeoSoft di ing. Windows 95,98,

NT

未公布 ⊙ ⊙ ⊙ 未知

Kriging C code 1995 Chao-yi Lang IBM RS/6000、

HP-700

未公布 ⊙ 免費

LYNX Geosystems

Inc 4.8.

1997-1998 Lynx Geosystems Inc Windows 95, NT 未公布 ⊙ ⊙ ⊙ ⊙ US$ 4,950

MacGeoPak 不祥 Distributed by RockWare, Inc. Mac 未公布 ⊙ ⊙ ⊙ ⊙ 495 US$ . MATLAB Kriging Toolbox 4.0. 1996 - Yves Gratton Centre Eau Windows 95 Matlab ⊙ ⊙ ⊙ 免費

Software

RC2 _{1999 Reservoir Characterization,}

Research and Consulting

(also known as (RC)2₎

UNIX, some PC 未公布 ⊙ ⊙ ⊙ 未知

SADA 2.0 2000 - 2002 Robert Stewart、Brad Lyon、

Tom Purucker

Windows 32bits 未公布 ⊙ ⊙ 免費

SaGA 1.0. 1996 Kirill K. Pankratov MATLAB Matlab ⊙ ⊙ 40 US$. SAGA GIS 1.0. 2003-2004 O. Conrad、 A.Ringeler、 J.

B. Windows 未公布 ⊙ ⊙ ⊙ ⊙ 免費 SAGE2001 1.04 SAGE95 (1997) -SAGE2001( 1999)

Edward H. Isaaks Windows 95, 98,

NT

未公布 ⊙ ⊙ US $515

SAND 1997 Norwegian Computing

Center Applied Research and

Development Norsk

Regnesentral

Unix 未公布 ⊙ ⊙ ⊙ 未知

SGS 2.3 1992-1999 Frieder Tonn. UNIX (SunOs,

AIX, Linux, ...)

Standard-C

⊙ ⊙ 免費

Sigma View 不祥 Landmark Graphics Corporation

Unix 未公布 ⊙ ⊙ 未知

Snowden VISOR 1997 - 2003 Snowden Group Windows 未公布 ⊙ ⊙ ⊙ 未知

Spatial Functions

for kriging and

point pattern

analysis for R and

SPlus 1997 Venables、 Ripley、W. N. Venables、B. D. Ripley Windows 或 UNIX SPlus 和 R ⊙ ⊙ 免費

SPATSTAT 2003 Adrian Baddeley、Rolf

Turner Windows 或 UNIX SPlus 和 R ⊙ ⊙ 免費

Spherekit 1.0 1997 Cort Willmott、Rob Raskin、

Chris Funk、Scott Webber、

Mike Goodchild UNIX (DEC Alpha、SGI Challenger) Fortran 和 C ⊙ 免費

Splus for Arcview

& S+Spatial Stats

1.5

1996-2000 Insightful MS Windows,

UNIX

version

1998)

SURFACE III 1995 Kansas Geological Survey

University of Kansas

MAC 未公布 ⊙ ⊙ 30 US$

SURFER 8 1997-2002 Golden Software, Inc. Windows 未公布 ⊙ ⊙ ⊙ 600 US$ Uncert 1.20b 1994-1997 William L. Wingle、Eileen P.

Poeter、Sean A. McKenna

UNIX C ⊙ ⊙ ⊙ 免費

Vario96 1996 John G. Davis DOS 未公布 ⊙ ⊙ ⊙ 未知 Variowin 2.21 1996-1999 Yvan Pannatier Shell

International Exploration and

Production

Win 3.x.、

Win95、WinNT

未公布 ⊙ ⊙ 免費

VESPER 1.0 1999 Minasny, B.、 McBratney,

A.B.、Whelan、B.M.

WIN95/98 & NT

(Service Pack 3)

未公布 ⊙ ⊙ 免費

於國內外成功將地理統計技術與地理資訊系統相互結合的軟體為美國 ESRI 的 ArcGIS Geostatistical Analyst 模組，以及 ArcObjects 中相關地理統計的元件。

以 ArcGIS Geostatistical Analyst 地理統計模組而言，其內容包含有以下部分：  數學統計分析部分包含有

（1） Inverse Distance Weighting,IDW （2） Global Polynomial Interpolation （3） Local Polynomial Interpolation （4） Radial Basis Functions,RBFs  地理統計部分包含有 （一）克利金的型態 （1）一般克利金（Odinary Kriging） 1. 一般克利金（Odinary Kriging） 2. 簡單克利金（Simple Kriging） 3. 泛克利金（Universal Kriging） （2）指標克利金（Indicator Kriging） （3）聯合克利金（Co- Kriging）

（二）克利金分析模式工具 （1）資料轉換方法 1. 對數轉換（Logarithmic） 2. Arcsine 3. 資料標準化（Normal score） 4. 平方根（Square root） （2）分析資料的趨勢 1. 整體性多項式分析 2. 區域性多項式分析 （3）自相關分析（Autocorrelation） 1. 變異圖分析（Variogram） 2. 交叉變異圖分析（Cross- Variogram） （4）變異量分析（Variography） 1. 選擇套用模式 a.球形模式（Spherical） b.指數模式（Exponential） c.高斯模式（Gaussian） 2. 半變異圖/共變異數（Semivariogram/Covariance）分析 3. 等方向或非等方性變異分析（Anisotropy） 4. 小採樣間隔內的變異或觀測值得誤差分析（nugget） 5. 共變異數函數(cross covariance)和共變異數圖分析 6. 最小加權平方數推估 7. 雲叢分布展示方式 a.網格方式 b.多邊形方式 8. 檢核資料二元分布情形 （5）驗證方法

1. 共變異量分析套配模式的優劣 2. 確認推估的品質 3. 交叉比對套用模式 4. 利用滑鼠顯示推估的數值 （三）推估結果包含 （1）推估結果的標準差分布圖 （2）推估結果分布圖 （3）指標誤差

（4）Quantile map (over- and underpredicted values)

以 ArcObjects 中相關地理統計的元件而言，其包含於空間分析模組中的 物件群中，主要相關於地理統計模組演算元件的內容包含有以下部分：

 IInterpolationOp

 IDW( geoDataSet as IGeoDataSet,power as Double,

radius as IRasterRadius,[Barrier as variant] ) as IgeoDataSet  Krige( geoDataSet as IGeoDataSet,

semiVariogramType as esriGeoAnalysisSemiVariogramEnum, radius as IRasterRadius,outSemiVariance as Boolean,

[Barrier as variant] ) as IgeoDataSet

參數包含有

 Spline( geoDataSet as IGeoDataSet,

[numPoints as variant] ) as IgeoDataSet 參數包含有  IRasterRadius  IRaster  IRasterLayer  IRasterAnalysisEnvironment

本研究計畫除針對上述 ArcGIS Geostatistical Analyst 所提供的地理統計 分析功能與 ArcObjects 所提供的地理統計演算元件進行本土化建置外，亦將 加入其他尚未開發或較不周全的地理統計演算方法，並透過建立操作範例的 方式來加強使用者操作界面之親和性與易懂性，使崧旭資訊公司 SuperGIS 產品能夠創造力與產業競爭力，成為亞洲最佳的地理資訊系統解決方案。 三、研究方法 （一）研究流程與分工 本計畫主要分為三大部分，第一部份為中國文化大學透過相關理論整理 與分析，然後篩選出適合做為構建地理統計演算元件之方法。第二部分為崧 旭資訊專業工程師利用 Visual C++程式語言撰寫相關地理統計元件與模組 功能介面。第三部份為中國文化大學透過實驗設計的方法，設計出各類問題 之測試例題，並與其他地理統計專業軟體比較，以評估演算法之計算績效。 整體而言，中國文化大學與崧旭公司之合作模式如下圖所示。

圖 1 研究流程與分工圖 於實驗設計部分，本研究計畫主要將程式設計所撰寫出來的演算法結 果，設計各類問題之測試例題，透過與當前可獲得地理統計相關專業軟體， 例如 GS+、GSLIB、GEOEAS、VARIOWIN 等，所使用之方法與其推估出 來的結果相互比較，了解其差異性或相同性，以評估演算法之求解的績效。 （二）方法 1、地理統計 （1）起源 地理統計（geostatistics）又稱為地質統計，於 1950 年代初期開始成形， 1960 年代法國統計學家 G. Matheron 進行許多相關基礎理論研究，形成一門 新的統計學分支。由於他首先是在地理學領域(如採礦學、地質學)中發展與 應用，因此被稱為地理統計學。發展的簡史整理如下： 1940 年代，地質學家進行礦場資源定量估計與評價時，發現傳統統計 學由於沒有考慮到地質樣本或區塊體積的大小，沒有考慮區位的空間變化特 徵，亦沒有考慮到空間相關性，因此在地質估計與評價上存在許多不能克服 的問題。後來發現，地質變量不是純粹的隨機變亮，因此不能用簡單的統計

中國文化大學

崧旭資訊

理論支援 演算法分析 元件功能分析 元件規格設定 整合測試 情境例題測試 解題績效分析 文件編寫 訓練教材講義 使用手冊 效率評估 評估機制建立 測試實驗設計 情境例題產生 程式設計 元件程式撰寫 操作界面設計 功能模組組合 回饋修正

方法進行估計與評價；同時數理統計中提出自相關和變異函數概念，成為空 間或時間連續性或相關的變量數據處理的基本方法。 1950 年代初期，南非的礦山工程師 D. G. Krige 與 H. S. Sichel 以南非 金礦儲量計算的問題出發，提出用樣本空間位置和相關程度來估計區塊位置 及存量，使其估計誤差最小的儲量計算方法，為 kriging 方法的雛形。G. Matheron 於 1962 年提出區域化變量概念，並出版”應用地理統計學導論”。 1970 至 1980 年代大量地理統計研究與應用書籍相繼出版，例如： M. David (1977)出版”礦產儲量地理統計學評價”，A. G. Journel(1978)出版”採礦地理統

計學”，I. Clark (1979)出版”實用地理統計學”， B.D. Ripley (1979)的”空間統

計學”等。此外， 1989 年 E. H. Issaks 與 R. M. Srivostava 合著” 應用地理統 計學導論”，主要介紹地理統計學應用的方法；1991 年 N. Cressie 出版”空間 數據統計學”，著重於理論證明與推導。近幾年來此理論已被廣氾應用於生 態、水文、環境污染、遙測、人文等空間資料分析、估計。 (2)地理統計(Geostatistics)定義、假設、條件 地理統計學主要是以區域化變量(regionalized variable)理論為基礎，研究 那些分佈於空間上並顯示出一定結構性與隨機性的自然現象。當一個變量呈 現空間分佈時，稱為”區域化”，這種變量常常反應某種空間現象的特徵，用 區域化變量描述的現象稱之為”區域化現象”。 地 理 統 計 (Geostatistics) 是 以 區 域 化 變 數 理 論 為 基 礎 ， 利 用 變 異 圖 (Variogram) 進 行 一 區 域 化 隨 機 變量 Z(x) 之 資 料 結 構 分析 (Data Structure

Analysis)而得到區域化變數空間相依(Spatial dependence)性，以克利金法

(Kriging Method) 推 估 未 採 樣 點 數 值 之 最 佳 線 性 無 偏 估 計 (Best Linear

Unbiased Estimate，簡稱 B.L.U.E.) 。由於地理統計學中，空間抽樣只能得 到一對樣本的數值 Z(x)與 Z(x+h)，不可能在空間上同一點取得重複，這就 在統計推論上出現了問題，為克服這個問題，必須對區劃變量 Z(x)做一些假 設：

若一隨機變量域之平均值、變異數及共變異數滿足下列條件，稱為二階 定常性。 (A)在整個研究區內，區域化變量 Z(x)的數學期望值對任意 x 存在且等 於常數。 平均值：E[Z(x)]=m ，m 為常數 (1) (B)在整個研究區內，區域化變量的共變異數函數對任意 x 與 h 存在且 平穩。 變異數：Var[Z(x)]=2 _，2_為常數 (2) 共變異數：Cov[Z(x), Z(x+h)] = E{Z(x)Z(x+h)}-E[Z(x)]E[Z(x+h)] = E[X(x)Z(x+h)]-m2=C(h) (3) 式中 h 為點 x+h 及 x 之距離，表示空間中任意兩點之共變異數與絕對位 置無關，只與相對距離有關。 B.內在假設(Intrinsic Hypothesis)： 若區域化變量 Z(x)滿足下面兩個條件，稱該區域化變量為內在假設或弱 二階定常性。 (A)在整個研究區域內，區域化變量 Z(x)的增量[Z(x)-Z(x+h)]的數學期望 值對任意 x 與 h 存在且等於零。 E[Z(x) - Z(x+h)] = 0 (4) (B)在整個研究區域內，區域化變量的增量[Z(x)-Z(x+h)]的共變異函數對 任意 x 和 h 存在且平穩。 Var[Z(x)-Z(x+h)]=E[Z(x)-Z(x+h)]2=2(x,h)=2(h) (5) 式中，(h)為變異圖(Variogram)或半變異圖(Semivariogram)。 C.正定條件（positive definite condition）

設 Z（x）為一平穩性隨機函數，其數學期望為 m，共變異數為 C（h），

 



 n i i iZ x Y 1  (6) 則 Y 也為一隨機變量，其變異數 Var[Y] 0≧ ，並可表示為：



( ) 0 1 1   



  n i n j j i j i C x x Y Var  (7) 當 C(Xi-Xj) 必須滿足上述各變異數 0≧ ，則稱 C(h)為正定。 因為



h =C(0)-C(h)，則：



( ) 0 1 1   



  n i n j j i j i x x Y Var  (8) 在(xi xj)_滿足條件



  n i i 1 0  的條件下，Var{Y} 0≧ ，稱 為「條件」



h 正定的。正定條件在分析共變異數函數和變異函數特性時非常有用。地理統 計學中研究的共變異數 C（h）必須是個非負定的函數，即由C(x_i x_j)構成 的共變數函數矩陣，必須是個非負定的矩陣。而變異函數必須是條件非負定 函數，即由-(x_i x_j)構成的矩陣必須是條件非負定矩陣，它的條件是



  n i i 1 0  成立。 （3）變異圖分析(Variography) 架設 Z(x)為區域化隨機變量，並滿足二階定常性與內在假設，h 為兩樣 本空間分隔距離，Z(xi)和 Z(Xi+h)分別是區域化變量 Z(x)在空間位置 xi 和 xi+h 上的觀測值(i=1,2,….,N(h))，其算式如下：





_    _{ } 



 2 ) ( ) ( ) ( 2 1 ) ( ) ( 1 x z h x z h N h h N i  (9) 以 h 為橫軸，γ(h)為縱座標，畫出變異函數曲線，用以表現區域化變量 Z(x)的空間變異特性，稱為實驗變異圖（如圖 2），為空間分析與結構分析 的有效工具。

圖 2 實驗變異圖示意圖

（4）一般典型變異圖重要參數

一般典型變異圖，具有三個重要之參數分別為：(1)Sill；(2) 影響範圍 (Range)；(3) 碎塊效應(Nugget Effect)；(h)最後會逼近一定值，此值稱為 Sill，此時之距離 h 稱為影響範圍，表示在此範圍外資料之間無相關性。另 一特徵是圖形在半變異數軸的截距，稱為碎塊效應，產生的原因可能是小距 離發生在小採樣間隔內的變異或是觀測值之實驗誤差。 圖 3 一般典型變異圖 （5）等向性與非等向性變異 地理統計利用變異圖表現空間變異，且變異圖模式為距離之函數。當空 間變數之變異圖僅為距離之函數，且不因空間方向不同而改變時，此稱為等 向性變異。若變異因空間方向不同而變化，則此現象稱之為非等向性變異。 換言之，非等向性變異之變異圖模式為距離與方向之函數。因此，可改寫為：





_





    ( , ) 1 2 ) , ( ) ( ) , ( 2 / 1 ) , (      N h i i i z x h x z h N h (10) ： 沿 點 xi 與 xi+h 之角度； N h( , )_{：沿 點 xi 與 xi+h 之角度樣本點間距為 h 之組數。} 若有二區域化變數 z 及 y，其聯合變異圖可由下式求得。



       ( ) 1 )] ( ) ( [ )] ( ) ( [ ) ( 2 1 ) ( ˆ h N i i i i i ZY z x z x h y x y x h h N h  (11) 一般而言，等向性變異與非等向性變異之特性判定之步驟相同，皆是由 沿各方向進行變異圖模式之套配，進而得到沿各方向之變異圖。換言之，可 得到沿各方向之 Sill 值、影響範圍與碎塊效應，最後即可判定該空間資料為 等向性或非等向性。若沿各方向之變異圖模式之 Sill 值、影響範圍與碎塊效 應皆相同，則為等向性。若不相等則為非等向性。非等向性變異依其結構不 同又可分為 幾何型 非 等向(Goemetric Anisotropy) 、區 域型非 等向 (Zonal Anisotropy)與混合型非等向性(Mixed Anisotropy)。 A.幾何型非等向變異 幾何型非等向之變異圖，其 Sill 值僅依距離變化而改變，不因方向而 變。換言之，於空間中於任何角度，當樣本點間距離相同時其 Sill 值亦相同， 但是影響範圍因角度不同而異。且影響範圍最大之方向為變異圖方向，此範 圍稱為主要影響範圍(Major Range)。此時之變異圖稱為該角度之變異圖 (Directional Variogram)。影響範圍最小者稱為最小影響範圍。若將各角度之 空間影響範圍繪於玫瑰圖上則呈現橢圓形。

B.區域型非等向性變異 區域型非等向之變異圖，其 Sill 值不僅依距離變化而改變，且因方向而 變。換言之，於空間中任一角度，當距離相同時其 Sill 值不相同，且影響範 圍因角度不同而異。區域型非等向之變異圖模式由兩個以上之非等向變異圖 組成。一般之自然現象較少出現此類之非等向性。 圖 5 區域型非等向變異圖 C.混合型非等向性變異 混合型非等向性變異為幾何與區域非等向性兩者之混合，其 Sill 值與影 響範圍皆因方向不同而變。自然界常出現此一空間變異。例如地質之變化常 於水平方向較垂直方向之影響範圍大，且垂直方向之變化較水平方向大。 圖 6 混合型非等向性變異 （6）變異圖理論模式 實際上變異函數模型γ(h)是未知的，往往要以有效的的空間取樣數據進 行估計，對各種不的 h 值可計算出一系列的 γ(h)值。截至目前為止，統計學 中最常用的模型包括：

A.球形模式(spherical model) 球形模型的一般公式為： a h a h C C a h a h C C h h                ,0 , ) 2 1 2 3 ( 0 , 0 ) ( 0 3 3 0  (12) 式中 C0 為碎塊常數，C+C0 為 Sill 值，C 為結構變異數，a 為影響範圍。 當 C0=0，C=1 時，稱為標準球形模型。在原點處的斜率為 3C/2a，與 Sill 值焦點的橫座標為 2a/3。球形模型是地理統計學應用最廣的理論模型。 圖 7 球形模式變異圖 B.指數模式(exponential model) 指數模型的一般公式為： (13) 當 h=3a 時，1ea 1e3 0.951 h ，即當 h=3a 時，(h)C0 C， 因此指數模型的影響範圍為 3a。當 C0=0，C=1 時，稱為標準指數函數模型。           _ 0 ), 1 ( 0 , 0 ) ( 0 C e h C h h a h 

圖8 指數模型變異圖 C.高斯模式(Gaussian model) 指數模型的一般公式為：           _ 0 ), 1 ( 0 , 0 ) ( 2 2 0 C e h C h h a h  (14) 當 h=γ(h)時 ， 1 ₂ 1 3 0.95 1.0 2        e a h e ， 即 當 h= 3a 時 ， C C h) ₀  (  ，因此高斯模型的影響範圍為 3 。當 C0=0，C=1 時，稱為a 標準高斯函數模型。 圖9 高斯模式變異圖

比較上述三個模型發現，球形模型的影響範圍最小，指數模型的影響範 圍最大，而高斯模型的影響範圍位於中間。球形模型與指數模型通過原點都 有切線存在，高斯模型切線不存在。 圖 10 球形、指數、高斯模型比較圖 D.線性 sill 模型 線性 sill 模型的一般公式為：            a h C C a h Ah h C h , 0 , 0 , ) ( 0 0  (15) 式中 A 為常數，表示直線的斜率。當 h=0 時，γ(h)=C0，當 0<h a≦ ，γ(h)=Ah 為一直線，當 h>a 時，γ(h)＝C0+C，因此這種模型的影響範圍為 a，sill 值 為 C0+C。 E.純碎塊效應模型 純碎塊效應模型的一般公式為：       0 , 0 , 0 ) ( 0 h C h h  (16) 式中 Co>0 為先驗變異數。這種模型相當於區域化變量為隨機分佈，樣 點間的共變異數函數對於所有距離均等於 0，也就是說變量的空間相關性不

圖 11 線性 sill 模型 圖 12 純碎塊效應模型 （7）空間估計 克利金法（kriging）又稱為空間局部估計或空間局部插值，是地理統計 學兩大主要內容之一。該方法是建立在變異函數理論及結構分析理論基礎 上，在有限區域內對區域化變量的取值進行無偏最佳化估計的一種方法。 克利金法與普通的估計不同點，在於其最大限度地利用空間取樣所提供 的各種信息。在估計未知樣點數值時，它不僅考慮落在該樣點的數據，而且 考慮鄰近樣點的數值；不僅考慮待估樣點與鄰近已知樣點的空間位置，並考 慮各鄰近樣點彼此間的位置關係，並利用已知觀測值空間分佈的結構特徵。 使得其估計比其他傳統的估計方法更精確，更符合實際環境，避免系統誤差 的出現。 常用的克利金估計方法有一般克利金 (Ordinary Kriging)、指標克利金 (Indicator Kriging)、聯合克利金 (Co-Kriging)等三類，以下分別詳述之。

A.普通克利金 在克利金法區域化變數 Z(x)的數學期望值為已知或未知，分別有兩種估 計方法，當E

 

Z(x) m，m 為已知常數時，此時的估計法稱為簡單克利金 法。當 m 為未知常數時，稱為普通克利金法。 (A)簡單克利金法 在簡單克利金法中，區域化變數 Z(x)的數學期望值為已知常數，即

 

Z x m E ( )  _，若 m x Z x Y( ) ( ) ，則E

  

Y(x) E Z(x)m



0，其共變異函數存在， 且為E



Y(x),Z(x)



C(x,y)。假設一待推估區塊為 V，其實際值為 Y(V)，中 心點為 x，則實際值可表示為：







    vY x dx _v vZ x dx m Z V m v V Y( ) 1 ( ) 1 ( ) ( ) (17) 若待估計區塊 V 的範圍內有一組 n 個已知樣本點x_i(i1,2,...,n)，其觀 測值為 Z(xi)，則Y(xi)Z(xi)m，Y(V)的估計值 * v Y _{是由這 n 個已知樣點觀} 測值 Y(xi)的線性組合，表示為：



 n i i i v Y x Y 1 * ) (  (18) 在滿足無偏性與最佳化的條件下， * v Y 成為 Y(V)的線性、無偏、最佳化 估計量，由此可得出簡單克利金方程組：



  n j i j i jC x x C x V 1 ) , ( ) , (  (19) 上式為一個 n 階方程組，可解出j(j1,2,...,n)，經整理後，可得簡單 克利金估計變異數 2 E  ：



    n i i i E C V V C x V 1 2 ) ( ) , (   (20)



 n i i i K Y x Y 1 * ) (  (21) 因為前面令 Y(x)=Z(x)-m，並非區域化變數 Z(x)於區塊 V 內的真實值 Z(V)。因此，Z(V)的克利金估計量應為：









         n i n i i i i i K K m Y m Y x m Z x m Z 1 1 * * ) ( ) (  









         n i n i n i i i i n i i iZ x m Z x m m 1 1 1 1 ) 1 ( ) ( ) (     因此





    n i n i i i i K Z x m Z 1 1 * ) 1 ( ) (   (22) (B)普通克利金法 當區域化變數 Z(x)的數學期望值E

 

Z(x) m，m 為未知常數時，採用 普通克利金法。假設待估計區塊為 V，中心為 x，其平均值為 Zv，則 dx x Z V Z V V 



( ) 1 (23)

 



_

 

 V V E Z x dx m V Z E 1 ( ) (24) 在待估計區塊 V 的區域範圍內，有一組已知 n 個樣點xi(i1,2,...,n)， 其觀測值為 Z(xi)，其數學期望值也為 m，即E

 

Z(x) m。令 * V Z 是 Zv 的線 性估計量，由 n 個已知的樣本觀測值 Z(xi)構成的線性組合，即：



 n i i i V Z x Z 1 * ) (  (25) 同樣在滿足無偏性與最佳化的條件下， * V Z 成為 Zv 的線性、無偏、最佳 化估計量。經整理後可得普通克利金估計變異數 2 K  ，即



    n i i i K C V V C x V 1 2 ) , ( ) , (    (26) 在變異函數存在條件下，也可用變異數函數表示普通克利金線性方程組 與普通克利金估計變異數，即：          





  n i i n j i j i i x x x V 1 1 1 ) , ( ) , (      (27)      



   ) , ( ) , ( 1 2 V V V x_i n i i K (28) 將解出的i(i1,2,...,n)帶入公式(9)中得出普通克利金估計量 * K Z ，即



 n i i i K Z x Z 1 # ) (  (29) 待估計區塊 V 區域範圍內以 xi 為中心的區塊 vi 時，普通克利金線性方 程組和普通克利金估計變異數為：          





  n i i n i i j i jC v v C v V 1 1 1 ) , ( ) , (    (30)

 

_

    n i i i K C V V C v V 1 2 ) , ( ,    (31) 也可以用變異函數表示為：          





  n i i n j i j i j v v v V 1 1 1 ) , ( ) , (      (32)



   n i i i K v V V V 1 2 ) , ( ) , (      B.聯合克利金 聯合克利金法是利用兩個變數間的相互關性，以其中易於觀測的變數對 另一變數進行局部估計的方法。聯合克利金法可顯著提昇估計精度及採樣效 率。但在實際應用中，聯合克利金法要求要有一個已知的互相關函數，這就 需要在很多地點同時採樣，測定兩個函數間的互相關性。聯合克利金法是建 立在聯合區域化變數理論基礎上，藉由建立共變異數函數(cross covariance) 和共變異數圖(cross variogram)模型，然後以聯合克利金法對未抽樣區欲進行 估計。 (A)共變異數函數理論 研究聯合區域化現象與研究單一變數現象方法相類似。假設 K 個聯合 區域化變數 Z1(x), Z2(x), …,Zk(x)，組成一組 K 為區域化變數的向量



Z₁(x),Z₂(x),...,Z_k(x)



。在觀測前，他是一個 K 維區域化向量，觀測後可視 為 K 維空間向量。在二階定常性假設條件下，聯合區域化變量為： (a)每一個Z_K(x)(K 1,2,...,k)的數學期望存在平穩，即

 

Zk x mk E ( )  (33) (b)每個區域化隨機變數 Zk(x)，Zk'(x)的共變異函數為：



'



' '( ) k( ) k( ) k k kk h E Z x Z x h m m C     (34) 式中K,K’=1,2,…,K，當 h≠0時，K,K’的順序不能任意顛倒。當 h=0 時，



'



' '(0) k( ) k ( ) k k kk E Z x Z x m m C    (35) 上式為 Zk(x)與 Zk’(x)的變異數，或是普通統計學中兩個變數之間的共 變異數。 (c)滿足內在假設條件下，Zk(x)的增量



Z_k(x)Z_k(xh)



的數學期望為

零，則每一對區域化變數 Zk(x)與 Zk’(x)的共變異函數存在，且為：









( ) ( ) ( ) ( )



2 1 ) ( _{' '} ' h E Zk x Zk x h Zk x Zk x h kk       (36) ' kk  與C 僅與 h 有關，與 x 無關。kk' (B)共變異數與共變異圖計算公式 在點 x 與點 x+h 處，分別測得兩個變量的觀測值 Zk(x),Zk’(x),Zk(x+h), Zk’(x+h)，則共變數計算公式為：











      ( ) 1 1 ' ' '( ) ( ) ( ) ( ) ) ( 1 ) ( ' h N j i k k j k j k k i k x Z x h Z x Z x h m m Z h N h k Ck (37) 式中







  n i i k k Z x h N m 1 ' ' ( ) ) ( 1 ，







  ( ) 1 ) ( ) ( 1 N h j j k k Z x h N m ，N(h)為樣本對數。 共變異圖的計算公式為：











       ( ) 1 1 ' ' ' ( ) ( ) ( ) ( ) ) ( 2 1 ) ( h N j i j k j k i k i k k k z x Z x h Z x Z x h h N h  (38) (C)聯合克利金線性估計量 聯合克利金法在理論上與普通克利金法本質相同，可以用推倒普通克利 金法的過程推導聯合克利金。假設研究區域上有 K 個變數構成聯合區域化 變數Z_K(x)(K 1,2,...,k)，滿足二階定常性與內在假設條件，即共變數與共 變異圖存在，即

















( ) ( ) ( ) ( )



2 1 ) ( ) ( ) ( ) ( ) ( ' ' ' ' ' ' h x Z x Z h x Z x Z E h m m x Z x Z E h C m x Z E k k k k k k k k k k k k k k          (39) 估計中心點為 x0，範圍為 Vk0(x0)的平均值，即



1 Z (x)dx Z

在範圍 Vk0(x0)附近有 n 個已知觀測樣本Zak(ak 1,2,...,nk)，其小區域 範圍為 vak, 則 Zak 的平均值，為



 zak k ak ak Z x dx v Z 1 ( ) (41) Zvk0 的估計量為Z_V*_K0，是 K 個聯合區域化隨機變數全部有效觀測值的 線性組合，聯合克利金線性估計量：



  k k n a ak k V k k Z Z 1 1 0 *  (42) (D)聯合克利金估計方程組與聯合克利金估計變異數 若在點 x0 處的某變數平均值為 u0，在 x0 附近有兩個已觀測的聯合區 域隨機變數u i(i=1,2,…,n)和vj(j=1,2,…,m)。則平均值 u0 的估計值 * 0 _構成 的聯合克利金線性估計量為：





   n i n j j i i iu bv a 1 1 * 0  (43) 其中 ai (i=1,2,…,n)，bj(j=1,2,…,m)為聯合克利金權重係數，為使 * 0 為 u0 的最佳無偏線性估計，則需滿足無偏性與最佳化條件，經整理後，兩個 變數的聯合克利金線性方程組之一般式為：                     













      m i i n i i n i m i j j i i j i i n i m i j j i i j i i b a v u C u v v C b v u C a u u C u u v C b u u C a 1 1 1 1 0 2 1 1 0 1 0 1 ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( (44) 上述為 n+m+2 階線性方程組，解該線性方程組可得聯合克利金權重係 數 ai 與 bj，然後帶入公式(43)，可得聯合克利金線性無偏最佳估計量 * 0 ，

) ( ] ) ( [ ] ) ( [ 0 ] ) ( [ 1 )] , ( [ z F z x z prob z x z prob z x z prob z x I E          此時的聯合克利金估計變異數為：





      n i m j j j i i COK C u u aC u u b C v u 1 1 0 0 1 0 0 2 ) , ( ) , ( ) , (   (45) C.指標克利金 指標克利金是一種非參數的地理統計學，其不必去除重要而實際存在的 高實際值數據，進行各項分析，並可求得在一定風險條件下，未知量 Z(x) 的估計量及空間分佈。 (A)指標函數 假設在空間範圍 D 中，於 x 位置取得樣本值隨機變數 Z(x)，選定一臨 界值 Z，則在 x 上的每一樣本點xZ_{上定義 Z 的指標函數為：} i(x,Z)=    , 0 , 1 在 D 上的任一區域AD_{內，低於 Z 值的隨機變數 Z(x)所佔區域 A 的} 比例表示為：

 



  Ai x Z dx A Z A I( , ) 1 ( , ) 0,1 (46) 高於 Z 值的隨機變數 Z(x)所佔區域 A 的比例表示為： ) , ( 1 ) , (A Z I A Z I 



Z x Z x A



P    1 ( ) , ₍₄₇₎ 在臨界值 Z 的條件下，隨機函數 I(x,Z)為二項式分配，其期望值為：



I x z



prob



Z x Z



E ( , )  ( ) I(x ,z)與 I2(x ,z)的期望期可以表示成下式 (48) 當 x 點上樣本值 x(x)≦Z 當 x 點上樣本值 x(x)＞Z

) ( ] ) ( [ ] ) ( [ 0 ] ) ( [ 1 )] , ( [ 2 2 2 z F z x z prob z x z prob z x z prob z x I E          ) ( )] , ( ) , ( [ )] , ( [ )] , ( ) , ( [ )] , ( [ )] , ( [ ]] , [ ) , ( [ ) , ( 2 2 z F z x I z h x I E z x I E z x I z h x I E z x I E z h x I E z x I z h x I E z h CI               ) ( ) ( )] , ( [ )] , ( [ )] , ( [ 2 2 2 z F z F z x I E z x I E z x I Var     )] , ( ) , ( [ ) ( )] , ( ) , ( [ )] , ( [ )] , ( ) , ( [ 2 / 1 ) , ( 2 2 1 z x I z h x I E z F z x I z h x I E z x I E z x I z h x I E z h             (49) 距離 h，兩位置的指標函數 I(x, z)與 I(x+h, z)其共變數、變異數和半變 異數分述如下： (50) (51) (52) (B)指標克利金方程組 假設空間範圍 D 中，有 N 個有效樣本點，在 D 範圍內之一區域AD_內 有 n 個有效樣本點



Z(x_),x_A,1,2,...,n



，在設定臨界值之後，可得樣 本的指標函數



Z(x_),x_A,1,2,...,n



，則I(A,Z)的估計值可以表示為：



 n Z i x Z Z A I 1 * ) , ( ) ( ) , (     (53) 若給定一系列臨界值 Zl



Z_l,l1,2,...,L



時，改寫為：



 n l l l Z i x Z Z A I 1 * ) , ( ) ( ) , (     (54) 為求算I*(A,Z)，在滿足無偏與估計變異數最小化的條件下，計算權重 係數_(Z)(1,2,...,n)，須解以下指標克利金方程組：

         





  n n i i Z Z A x Z x x Z 1 1 1 ) ( ) , , ( ) , , ( ) (             (55) n ,..., 2 , 1   指標克利金的共變異數 2 1 K  表示為：



   n i i K Z x A Z A A Z 1 2 1 ( ) ( , , ) ( , , )         (56) 以共變異函數可表示為：          





  n n i i Z Z A x C Z x x C Z 1 1 1 ) ( ) , , ( ) , , ( ) (           (57) n ,..., 2 , 1  



    n K C A A Z Z C x A Z 1 2 1 ( , , ) ( ) ( , , )       (58)

四、結果與討論

本研究應用 Visual C++程式開發語言撰寫空間統計分析模組，以精靈導引 的方式，進行相關克利金的變異圖分析與推估，並以擴充模組的方式，外掛於 SuperGIS 桌上型地理資訊系統上。

（一）開啟外掛 SuperGIS SpatialStatistic Analyst 模組

本研究成果透過 Visual C++程式開發語言，將空間統計分析模組的成 果架構在 SuperGIS 桌上型地理資訊系統中，以空間統計分析工具列呈現於 系統中，如圖 13 所示。 圖 13 空間統計分析模組工具列 空間統計分析模組主要提供三種克利金分析方式，分別為一般克利 金、聯合克利金與指標克利金三種方式。於一般克利金與指標克利金分析 中，主要提供單一圖層中的單一屬性變數的變異分析，即只能選擇單一資 料來源；聯合克利金分析方式中，則使同用一圖層中的兩組屬性變數進行 分析，即選擇同一圖層的兩組屬性資料來源，如圖 14 所示。 圖 14 匯入分析資料來源與選取分析方式

（二）等向性與非等向性變異 地理統計利用變異圖表現空間變異，且變異圖模式為距離之函數。 當空間變數之變異圖僅為距離之函數，且不因空間方向不同而改變時，此 稱為等向性變異。若變異因空間方向不同而變化，則此現象稱之為非等向 性變異。等向性與非等向性的變異結果呈現，如圖 15 與圖 16 所示，使用 者可以依據手動與鍵入數值方式，了解該筆數值在不同的角度內，其變異 數值的分布與方向性；此外，透過調整容忍角度來決定鄰近點在 Bandwidth 範圍內是否被列入搜尋目標。透過 Bandwidth 決定哪些配對點要繪製於半 變異圖時所需的搜尋寬度。擷取各方向性的資料進行變異數分析，使分析 結果更接近於實際，目前 Bandwidth 的寬度預設值為 6。 透過 Surface 該區塊中所呈現的圖示結果，提供使用者更了解當改變 搜尋資料的方向、角度與容忍寬度時，被截取到的顯示資料，依據所搜尋 到的兩兩變數之間的變異數值與方向值而有所不同，顯示於圖面上的資料 結構亦隨之改變，如圖 15 與圖 16 所示。由圖面上所呈現之結果，提供使 用者更容易了解哪種方向與角度之資料結構較具有相關性。 圖 15 未設定方向性角度 圖 16 設定方向性角度 系統於等向性與非等向性上提供三種參數，供使用者自行調整資料 搜尋計算的方向結構，如圖 17 所示，分述如下： （1）Angle Direction

用者依據資料結構之不同，透過方向的界定找出同一結構組之資料。 （2）Angle Tolerance 設定資料結構方向性的容忍範圍，即θ角度方向上，上下容忍某一 α 角範圍。在此範圍內之角度皆視為同一方向性之資料結構。 （3）Bandwidth 即搜尋資料隨搜尋角度無限擴大時，設定某一寬度大小時即收斂不 再擴大搜尋範圍。 圖 17 方向性設定參數示意圖 （二）變異圖分析 （1）半變異圖分析 針對兩兩資料進行變異函數γ(h)計算，依據不同的 h，將所有資料加 以匯整入統計圖表中加以呈現，橫座標顯示兩兩變數之間的距離，縱座標 則顯示兩兩變數之間所計算出來的變異函數γ(h)值，及統計圖上每一點之 值，為一組變數計算出之變異函數γ(h)值，實際成果展示如圖 18 所示。 圖 18 半變異圖分析圖

（2）實驗半變異圖分析 所謂實驗半變異圖，即將某一區間之所有變異函數值，加以平均獲 得一個平均的變異函數值，本研究主要將所擷取範圍內之資料分為 10 個 區間，每一區間獲得一個平均的變異函數值，加以相互連接成一曲線，以 利套配模式，如圖 19 所示。 圖 19 實驗半變異圖 （三）模式套配 本研究成果主要提供三種模式套配上，分別為球型模式、高斯模 式與指數模式。使用者依據不同的使用需求，自行選用較適宜的套配模 式。系統會依據你所選擇的套配模式，自行先給予一個最佳的套配結果 （Best Fit），使用者亦可以透過手動的方式，自行調整相關參數，以 獲得更好的套配模式解。可調整參數分別為： （1）碎塊效應（Nugget Effect） 圖形在半變異數軸的截距，稱為碎塊效應，如圖?所示，產生 的原因可能是小距離發生在小採樣間隔內的變異或是觀測值之實 驗誤差。系統會自行以最佳線性無偏估計方法先給予一組最佳的 參數值，使用者依據自行對資料結構的了解，手動調整該參數值 的數值，如圖 20 所示。 （2）影響範圍（Major Range） (h)最後會逼近一定值，此值稱為 Sill，此時之距離 h 稱為影 響範圍，表示在此範圍外資料之間無相關性，如圖 20 所示。透過

（3）Partial Sill (h)最後會逼近一定值，此值稱為 Sill，該參數結構示意如圖 所 示，如圖 20 所示。 圖 20 參數示意圖 （4）區間間隔大小（Lag Size） 區間間隔的大小，依據已知樣本點的最遠相關性計算範圍而 定，即距離（h）＝區間間隔大小 × 區間各數（Distance（h）＝

Lag Size（h） × Number of Lag）。系統會預先依據計算已知兩兩

樣本點間最大的範圍（Distance），預設的區間間隔數（Number of Lag）為 10，因此，系統自行給予的區間間隔大小（Lag Size）參 數為 Lag Size（h）＝Distance（h）/ Number of Lag，如圖 21 所示。 （5）區間間隔數（Number of Lag） 區間間隔數的數量隨使用者的需求自行設定，預設個數為 10，主要影響實驗半變異圖所呈現的結果與模式套配的結果，如 圖 21 所示。 （6）鄰近計算涵蓋樣本數（Neighbors to Include） “Neighbors to Include”該參數，主要針對克利金推估時，每一 已知樣本點，在某一範圍內，只少必須搜尋到已知點的各數，如 圖 21 所示。

Sill =C C₀= Total Variability

C= Structured Variability Range 0 C =Nugget Distance （h） γ（h）

圖 21 變異圖參數設定介面 （四）克利金推估

本研究提供之克利金估計方法有一般克利金 (Ordinary Kriging)、指標 克利金 (Indicator Kriging)、聯合克利金 (Co-Kriging)等三類。透過上述之樣 本分析、模式套用與校正驗證等流程後，如圖 22 所示，依據不同的推估方 法而獲得最後的分析結果。 A B C D E F

最後推估之結果以網格式圖層呈現，透過與地理資訊系統的結合，將該網 格圖層展示，並可套疊其他向量資料進行分析研究，如圖 23。 圖 23 一般、指標與聯合克利金推估結果 （五）統計分析 透過 Cross Validation 介面進行推估結果驗證，將已知實際值與推 估值相互比較，獲得誤差值，以了解推估結果之正確性，如圖 24 所示。 圖 24 推估結果正確性驗證

（六）實驗設計

將相同的空間資料分別帶入 SuperGIS SpatialStatistical Analyst（以下簡 稱 SS）、GS+（以下簡稱 GS）與 ArcGIS GeoStatistical Analyst（以下簡稱

AG）三套軟體中進行不同的模式套配，如圖 25、圖 26 與圖 27 所示。結果 呈現，SS 與 GS 在介面上以實驗半變異圖呈現，整體趨勢大致相同。 （a） （b） （c） 圖 25 球型模式套配（a）SS（b）GS（c）AG （a） （b） （c） 圖 26 指數模式套配（a）SS（b）GS（c）AG （a） （b） （c） 圖 27 高斯模式套配（a）SS（b）GS（c）AG 透過球型模式、指數模式與高斯模式三種不同的模式套配結果與變異 函數計算，所獲得的 Nugget、Sill 與 Range 三個參數帶入 GS+統計驗證中， 獲得的 r2_{與 RSS 值進行比較，如表 4、表 5 與表 6，結果顯示，相同資料}

結果為最佳，其次為 SuperGIS SpatialStatistical Analyst 模組所獲得的套配結 果。 表 4 球型檢驗統計表 球型模式 r2 RSS SS 0.799 42.4 GS 0.803 41.4 AG 0.797 44.5

註：SS:SuperGIS SpatialStatistical Analyst; GS:GS+; AG:ArcGIS GeoStatistical Analyst

表 5 指數檢驗統計表

指數模式 r2 RSS

SS 0.8 46.7

GS 0.816 38.7

AG 0.795 36.1

註：SS:SuperGIS SpatialStatistical Analyst; GS:GS+; AG:ArcGIS GeoStatistical Analyst

表 6 高斯檢驗統計表

高斯模式 r2 RSS

SS 0.748 57.7

GS 0.785 45.2

AG 0.769 119

註：SS:SuperGIS SpatialStatistical Analyst; GS:GS+; AG:ArcGIS GeoStatistical Analyst

五、結論與建議

國人結合學術單位與產業界自行研發的地理統計分析模組（SuperGIS Spatial Statistical Analyst）已成功上市，目前已使用於環境工程與觀光兩方面， 亦有許多相關單位正在測試軟體，同時獲得各方的肯定。此軟體的執行效能與 分析結果，亦能於國外 GIS 軟體並駕齊驅，由於空間統計分析模組的加入， 使得國人自行研發的桌上型地理資訊系統 SuperGIS，其對於空間的分析能力 更能與國外所開發知 GIS 軟體相互競爭。

於環境工程領域方面，目前主要應用於地下水水質監測，透過地理統計 分析模組的推估分析功能，提供管理者在有限的調查地下水井中，簡便、有效 率且準確的找出該研究區所有未採樣點的推估數值，大大降低龐大外業調查的 困擾與增強全面性研究區水質狀態的評估依據，提升水質監測的時效性。 於觀光方面，目前應用於生態旅遊景點規劃結果，群眾對生態的影響效 應。透過有限點位的動物、植物或昆蟲等的實際調查，了解該研究區內的生態 分部或遷徙情形，進而了解人類的活動對該區生態的實質影響。提供管理者更 完善且全面性的評估決策依據，及時改善不良的規劃設計，永續生態與環境。 本計畫成功的建立起產學合作模式，將學術理論及國外經驗轉換為國人 自行研發之空間統計模組，有助於未來空間統計應用於國內環境調查、規劃與 監測之可行性，同時大幅提昇學術研究的能力及軟體研發能力，未來可依不同 的研究主題與特性，設計適合的軟體功能與介面，不再受限於國外軟體制式的 架構與功能。 六、參考文獻

1、 Cheng, K. S., Wei, C. and Chang, S. C. (2004). Locating landslides using multi-temporal satellite images. Advances in Space Research, 33:296-301. 2、 Clark, M. L., Clark, D. B. and Roberts, D. A. (2004). Small-footprint lidar

estimation of sub-canopy elevation and tree height in a tropical rain forest landscape. Remote Sensing of Environment, 91:68–89

3、 Cressie N (1993) Statistics for spatial data. John Wiley & Sons, New York 4、 Dalezios, N. R., Loukas, A. and Bampzelis, D. (2002). Spatial variability of

reference evapotranspiration in Greece. Physics and Chemistry of the Earth, 27:1031-1038

5、 Deutsh C. and A. G. Journel, GSLIB, Geostatistical Software Library and User’sGuide,Oxford University Press.,New York,1992.

6、 Diodato, N. and Ceccarelli, M. (2004). Multivariate indicator Kriging approach using a GIS to classify soil degradation for Mediterranean agricultural lands. Ecological Indicators, 4:177-187

7、 ESRI (Environmental systems Research Institute) (1998) Arcview 3.0; ESRI, Redlands, California

8、 Franklin, R. B. and Mills, A. L. (2003). Multi-scale variation in spatial heterogeneity for microbial community structure in an eastern Virginia agricultural field. FEMS Microbiology Ecology, 44:335-346

Version 2.3. Gamma Design Software, Plainwell, MI, 1995.

10、 Hudak, A. T., Lefsky, M. A., Cohen, W. B. and Berterretche, M. (2002). Integration of lidar and Landsat ETM+ data for estimating and mapping forest canopy height. Remote Sensing of Environment, 82:397-416.

11、 Isaaks E. H. and A. G. Journel, Applied Geostatistics, Oxford University Press Inc.

12、 Jaquet, O. and Carniel, R. (2001). Stochastic modelling at Stromboli: a volcano with remarkable memory. Journal of Volcanology and Geothermal Research, 105: 249-262.

13、 Jaquet, O. and Carniel, R. (2003). Multivariate stochastic modelling: towards forecasts of paroxysmal phases at Stromboli. Journal of Volcanology and Geothermal Research, 128:261-271.

14、 Journel, A. G. and C. J. Huijbregts, 1978, Mining Geostatistics, Academic Press. .

15、 Lin Y.P., T.K. Chang and T.P. Teng, 2001, Characterization soil lead by comparing of sequential Gaussian simulation, simulated annealing simulation and kriging methods, Environmental Geology, 41(1/2), 189-199. (SCI, EI) NSC 87-2621-P-002-012

16、 Lin, Y.-B., Y.-C. Tan, Y.-P. Lin*, C.-W. Liu , 2004, Geostatistical method to delineate anomalies of multi-scale spatial variation in Hydrogeological changes due to the ChiChi earthquake in the ChouShui River alluvial fan in Taiwan, Environmental Geology . (Accepted)

17、 Lin, Y.P. , 2002, Multivariate geostatistical methods to identify and map spatial variations of soil heavy metals, Environmental Geology 42, 1-10. (SCI, EI) NSC 87-2621-P-002-012 and NSC 89-2621-B-002-004.

18、 Lin, Y.P. and S. Rouhani, 2001, Multiple-Point Variance Analysis for Optimal Adjustment of A Monitoring Network, Environmental Monitoring and Assessment, 69(3), pp. 239-266.

19、 Lin, Y.P. and T. K. Chang, 2000, Simulated annealing and kriging method for identifying the spatial patterns and variability of Soil heavy metal, Journal of Environmental Science and Health, Part A - Toxic/Hazardous Substance & Environmental Engineering, A35(7), 1089-1115.

20、 Lin, Y.P. and T.K. Chang, 2000, Geostatistical Simulation and Estimation of the Spatial Variability of Soil Zinc, Journal of Environmental Science and Health, Part A - Toxic/Hazardous Substance & Environmental Engineering, A35(3), 327-347.

21、 Lin, Y.P., 2004, Assessing Impacts of Spatial Variability on Geostatistical Optimally Adjusting an Environmental Monitoring Network. (Submitted to

Environmental Monitoring and Assessment

22、 Lin, Y.P., C. C. Li, and Y. C. Tan, 2000, Geostatistical Approach for Identification of Transmissivity Structure at Dulliu area in Taiwan, Environmental Geology, 40 (1/2), 111-120.

23、 Lin, Y.P., T.K. Chang, C.W. Shi and C.H. Tsen, 2002, Factorial and indicator kriging method using geographic information system to delineate spatial variations and pollution sources of soil heavy metals. Environmental Geology, 48(2), 900-909.

24、 Lin, Y.P., T.P. Teng and T.K. Chang, 2002, Multivariate Analysis of the soil heavy metal pollution and Landscape patterns analysis in Changhua County in Taiwan, Landscape and Urban Planning, 62(1), 19-35.

25、 Lin, Y.P., Y.C. Tan and S. Rouhani, 2001, Identifying Spatial Characteristics of Transmissivity using Simulated Annealing and Kriging Methods, Environmental Geology, 41(1/2), 200-208.

26、 Lin, Y.P.,, D. G. Chian, S.H. Lin and T.K. Chang, 2004, Assessing impacts of ChiChi earthquake and typhoons on landscape patterns of the ChenYuland watershed in Central Taiwan. (To be submitted to Landscape and Urban Planning)

27、 Nanos, N. and Montero, G. (2002). Spatial prediction of diameter distribution models. Forest Ecology and Management, 161:147-158.

28、 Nanus, L., Campbell, D. H., Ingersoll, G. P., Clow, D. W. and Mast M.A. (2003). Atmospheric deposition maps for the Rocky Mountains. Atmospheric Environment, 37:4881-4892.

29、 Ozturk, C. A. and Nasuf, E. (2002). Geostatistical assessment of rock zones for tunneling. Tunnelling and Underground Space Technology, 17:275–285.

30、 Pannatier Y (1996) VARIOWIN Software for spatial data analysis in 2D. Springer-Verlag, New York.

31、 Rubel, F. and Hantel, M. (2001). BALTEX Precipitation Analysis: Results from the BRIDGE Preparation Phase. Phys. Chem. Earth (B), 26(5-6):397-401. 32、 Tsun-Kuo Chang, Guey-Shin Shyu, Yu-Pin Lin, and Nan-Chang Chang, 1999,

Geostatistical Analysis of Soil Arsenic Content in Taiwan, Journal of Environmental Science and Health, Part A-Toxic/Hazardous Substance & Environmental Engineering, A34(7), 1485-1501.

33、 Tsun-Kuo Chang, Nan-Chang Chang, Guey-Shin Shyu, and Yu-Pin Lin, 2000, Mapping Soil Mercury in Taiwan Using Geostatistics and Geographic Information Systems, J. of the Chinese Institute of Environmental Engineering, 10(2), 77-85.