第八講空間插值與地統計一階方法倒距離權重(IDW)插值趨勢面分析二階方法區域化變量協方差函數,變異函數,變異函數模型交叉驗證地統計：克里金(Kriging)方法概述：理解不同的克里金方法幾種不同的克里金方法

普通克里金(OK)、簡單克里金(SK)、泛克里金(UK)

指示克里金（IK）、協克里金（Ck）整理ppt1第八講空間插值與地統計一階方法整理ppt1一、一階方法倒距離權重（IDW）插值趨勢面分析整理ppt2一、一階方法倒距離權重（IDW）插值整理ppt2IDW插值方法假定每個輸入點都有著局部影響，這種影響隨著距離的增加而減弱。步驟:計算未知點到所有點的距離；計算每個點的權重：

權重是距離倒數的函數。計算結果：

1.倒距離權重(IDW)插值整理pptzhaoy@

zhaoy@IDW插值方法假定每個輸入點都有著局部影響，這種影響隨∑

=14

最大6

最小示例：IDW插值（求圖中0點的值）整理pptzhaoy@

zhaoy@∑=1示例：IDW插值（求圖中0點的值）整理ppt整理pptzhaoy@

zhaoy@整理pptzhaoy@zhIDW插值的一般模型

對所有或選定的i進行計算典型地，α

的取值是1或2除以距離權重的和，保證了權重加起來等于1整理pptzhaoy@

zhaoy@IDW插值的一般模型整理pptzhaoy@henu.e432100160200例：倒距離加權(IDW)插值結果整理pptzhaoy@

zhaoy@432100160200例：倒距離加權(IDW)插值結果整理IDW插值的缺點IDW不能得到大于樣本最大值或小于樣本最小值的估計。對于高程表面，這‘抹平’了峰和谷。（除非它們的高點和低點是樣本的一部分。因為估計值為均值，得到的表面將不通過樣本點整理ppt8IDW插值的缺點IDW不能得到大于樣本最大值或小于樣本最小值IDW插值在ArcGIS中的實現整理ppt9IDW插值在ArcGIS中的實現整理ppt9IDW插值在ArcGIS中的實現整理ppt10IDW插值在ArcGIS中的實現整理ppt10Idaho州降雨量等直線圖IDW插值在ArcGIS中的實現整理ppt11Idaho州降雨量等直線圖IDW插值在ArcGIS中的實現整2.趨勢面分析趨勢面分析(trendsurfaceanalysis)：用數學模型來模擬（或擬合）地理數據的空間分布及其區域性變化趨勢的方法。整理ppt122.趨勢面分析趨勢面分析(trendsurfacea確定性插值函數擬合最普通的形式：多項式e.gy=ax2+bx+c局部（local）：分段地全局（global）local1storderglobal1storder整理ppt13確定性插值函數擬合local1storderglobal趨勢面分析是一種整體內插法，該方法假設一般趨勢與局部變化無關，并利用曲面方程來模擬待定點附近地形表面的一般趨勢。通常使用的是1次、2次、3次趨勢面，過高次的趨勢面不利于反映空間趨勢，并可能存在趨勢面的“畸變”。其中，2次趨勢面可用待定點附近的6個數據點來計算方程式系數。整理ppt14趨勢面分析是一種整體內插法，該方法假設一般趨勢與局部變化無關DeterministicSolutionsFirstOrderPolynomialInterpolationPredictedModelMeasuredSecondOrder(third,fourth,etc.)

PolynomialInterpolationLocalPolynomialInterpolationRadialBasisFunction(Spline)Interpolation整理ppt15DeterministicSolutionsFirstO趨勢面的性質與特點是一種光滑的數學曲面，能集中地代表地理數據在大范圍內的空間分布變化趨勢。與實際地理曲面不同，它只是實際曲面的一種近似。實際曲面包括趨勢面和剩余曲面兩部分，即：實際曲面＝趨勢面+剩余曲面整理ppt16趨勢面的性質與特點整理ppt16設Zi(xi，yi)表示某一地理特征值在空間上的分布。其中(xi，yi)為平面上點的坐標。任一觀測點Zi可分解為兩個部分，即：趨勢面剩余面整理ppt17設Zi(xi，yi)表示某一地理特征值在空間上的分布。其中(一階趨勢面FirstOrderTrendSurface一階趨勢面殘差ResidualsfromFirstOrderTrendSurface整理ppt18一階趨勢面一階趨勢面殘差整理ppt18趨勢面參數的確定（最小二乘法）使每一個觀測值與趨勢值的殘差平方和最小，即按建立多元線性方程的方法，使Q對系數b0，b1，…，bn求偏導，并令這些偏導數等于零，得趨勢面的正規方程組，解正規方程組，即可求出系數，從而得到趨勢面方程。整理ppt19趨勢面參數的確定（最小二乘法）整理ppt19因為任何函數在一定范圍內總可以用多項式來逼近，并可調整多項式的次數來滿足趨勢面分析的需要，一般來說，多項式的次數越高則趨勢值越接近于觀測值，而剩余值越小。整理ppt20因為任何函數在一定范圍內總可以用多項式來逼近，并可調整多項式多項式趨勢面的數學模型整理ppt21多項式趨勢面的數學模型整理ppt21整理ppt22整理ppt22二階趨勢面SecondOrderTrendSurface(1,x,y,x2,y2,xy)T三階趨勢面ThirdOrderTrendSurface

(1,x,y,x2,y2,xy,x3,y3,x2y,xy2)T整理ppt23二階趨勢面三階趨勢面整理ppt23趨勢面的具體計算方法與步驟：原始數據列表；等間隔選取縱橫坐標網，將原始數據點入坐標；按多元線性回歸分析方法求出趨勢面的正規方程組，解出參數；從趨勢值等值線圖中，獲得地理要素的區域性變化規律；用F分布對趨勢面進行擬合程度檢驗。整理ppt24趨勢面的具體計算方法與步驟：原始數據列表；整理ppt24即用雙三次多項式擬合趨勢面。雙三次多項式（樣條函數）插值整理ppt25即用雙三次多項式擬合趨勢面。雙三次多項式（樣條函數）插值整理該曲面模型有16個待定系數（Cij,i,j=0,1,2,3;）。通常用4個數據點（規則格網的4個頂點）的4個函數值組成的4×4方程組求解（如圖）。這4個函數值是高程Z、x方向斜率R、y方向斜率S，以及扭矩T：56781234916151011121314整理ppt26該曲面模型有16個待定系數（Cij,i,j=0,1,2,3其中Z保證曲面通過格網的4個數據點，R、S、T保證曲面在這4個數據點處光滑連續。雙三次多項式（樣條函數）內插法是規則格網插值的常用方法之一。這種方法通過一系列曲面片段來拼接地形表面，最終得到一個1階、2階連續的表面。該方法屬于局部插值，計算負擔中等；對于平滑表面擬合效果最好，對于起伏的表面擬合效果最差。整理ppt27其中Z保證曲面通過格網的4個數據點，R、S、T保證曲面在這4趨勢面分析在ArcGIS中的實現整理ppt28趨勢面分析在ArcGIS中的實現整理ppt28趨勢面分析在ArcGIS中的實現整理ppt29趨勢面分析在ArcGIS中的實現整理ppt29趨勢面分析在ArcGIS中的實現Idaho州降雨量等直線圖整理ppt30趨勢面分析在ArcGIS中的實現Idaho州降雨量等直線圖整二、二階方法IDW插值和趨勢面方法的缺陷：IDW：距離權重函數的選擇和“鄰居”的定義是其“致命缺陷”（Achilles’heel）。趨勢面分析：Inasense,trendsurfaceanalysisletsthedataspeakforthemselves,whereasIDWinterpolationforcesasetstructureontothem.整理ppt31二、二階方法IDW插值和趨勢面方法的缺陷：整理ppt31二、二階方法區域化變量協方差函數&

半變異函數整理ppt32二、二階方法區域化變量整理ppt32地統計學是以區域化變量理論為基礎，以變異函數為主要工具，研究那些在空間分布上既有隨機性又有結構性，或空間相關和依賴性現象的學科。協方差函數和變異函數是以區域化變量理論為基礎建立起來的地統計學的兩個最基本函數。地統計學的主要方法之一，克里金方法(Kriging)就是建立在變異函數理論和結構分析基礎之上的。整理ppt33地統計學是以區域化變量理論為基礎，以變異函數為主要工具，研究當一個變量的取值與其空間位置有關時，就稱為區域化變量（regionalizedvariable）。區域化變量常常反映某種空間現象的特征，用它來描述的現象稱之為區域化現象。區域化變量，亦稱區域化隨機變量，Matheron(1963)將它定義為以空間點x的三個直角坐標為自變量的隨機場區域化變量具有兩個最顯著，也是最重要的特征：隨機性和結構性。整理ppt34當一個變量的取值與其空間位置有關時，就稱為區域化變量（reg隨機變量隨機函數隨機過程隨機場區域化變量與時間有關的隨機函數帶有多個(2個以上)自變量的隨機函數以空間點的三個直角坐標為自變量整理ppt35隨機變量隨機函數隨機過程隨機場區域化變量與時間有關的隨機函數http://cg.ensmp.fr/Presentation/Matheron/Matheron_en.shtml

ProfessorGeorgesMatheron(1930-2000.8.7)法國數學家和地質學家

整理ppt36http://cg.ensmp.fr/Presentatio區域化變量的功能：

由于區域化變量是一種隨機函數，因而能同時反映空間變量的結構性和隨機性。一方面，當空間點x固定后，Z(x)就是一個隨機變量，這體現了其隨機性。另一方面，在空間兩個不同點x與x+h處的區域化變量值具有某種程度的相關性，這體現了其結構性。整理ppt37區域化變量的功能： 由于區域化變量是一種隨機函數，因而能同時區域化變量的組成部分

數據點結構性可以用均值和常數趨勢表示空間相關數據通常呈現正空間相關性

隨機性測量誤差，其他誤差

整理ppt38區域化變量的組成部分數據點結構性整理ppt38distanceelevation結構性隨機性實際值整理ppt39distanceelevation結構性隨機性實際值整協方差函數與變異函數數學期望、方差和協方差數學期望：一階原點矩方差：二階中心矩協方差：二階混合中心矩整理ppt40協方差函數與變異函數數學期望、方差和協方差整理ppt40協方差函數 類似地，當Z(x)是區域化變量時，對于任意兩點si和sj，空間隨機過程的協方差函數為： 相關系數和方差分別定義為：整理ppt41協方差函數整理ppt41整理ppt42整理ppt42若，則過程是二階平穩的，即均值與方差獨立于空間位置并在研究區域上是常數。于是有：

稱為協方差圖或過程的協方差函數，稱為相關圖或相關函數。顯然，協方差函數僅依賴于向量差h，當h=0時，整理ppt43若若獨立性僅是距離的函數，與方向無關，則空間過程是各向同性(isotropy)的。協方差函數就只依賴于距離向量h：

若在給定距離和方向上，不同位置數值差異的均值和方差為常數，則下式成立：半變異函數，也稱半方差函數。(semi-variogram)整理ppt44若獨立性僅是距離的函數，與方向無關，則空間過程是各向同性(i半變異函數(semi-variogram)—區域化變量的基本研究工具

—半變異函數就是區域化變量增量平方的數學期望之半區域化變量在i、i+h點的值步長為h的樣品對數步長(h)：在一定方向上，距離為h的矢量整理ppt45半變異函數(semi-variogram)區域化變量在i、i方差（Variance）變異函數（Variogram）h對比：方差&變異函數整理ppt46方差（Variance）h對比：方差&變異函數整理ppthaVetordistanceh建立經驗半變異函數(Semivariogram)測度空間變異。對于間隔距離為h的每一對Z(x)和Z(x+h)，測度它們之間差的平方。??h??h??h??h??h??h整理ppt47haVetordistanceh建立經驗半變異函數(Se半變異云圖半變異云圖：數據集中所有點對的“半變異函數—距離圖”。

如果空間相關性存在：

互相靠近的已知點趨向于具有較小的半方差

互相遠離的已知點趨向于具有較大的半方差整理ppt48半變異云圖半變異云圖：數據集中所有點對的“半變異函數—距離圖半變異云圖（semivariogramcloud）整理ppt49半變異云圖（semivariogramcloud）整理pp分組(Binning)如果數據很大，樣點對的數目將迅速增加并且變得難以操作。

因此，可以將樣點對分組，也就是步長(h)分組。如，可將距離在0～1之間的作為第一組；h在1～2之間的作為第二組；依此類推。Binningisaprocessthataveragessemivariancedatabydistanceanddirection(hey,itsweighted!!!)整理ppt50分組(Binning)如果數據很大，樣點對的數目將迅速增加并Binnedsemivariogram整理ppt51Binnedsemivariogram整理ppt51模型擬合

半變異函數必須擬合一個數學函數或模型，以用于估計任何給定距離的半方差。整理ppt52模型擬合半變異函數必須擬合一個數學函數或模型，以用于估計任理論變異函數模型實踐中，常用的是變異函數圖：偏基臺值:C(partialsill)塊金值:C0(nugget)變程:a(range)h基臺值(sill)notrelatedanymore變程范圍內才有結構性變化（有規律的變化）反映隨機性大小：主要來源于區域化變量Z(x)在小于抽樣尺度h時所具有的內部變異；另外還有抽樣分析誤差。變異函數是一個單調不減函數。當h超過某一個范圍，例如變程，變異函數不再增大，而是趨于一個極限值，即為基臺值。實際上等于區域化變量的先驗方差。即，即基臺值與塊金值之差，表示數據中存在空間相關性引起的方差變化范圍。整理ppt53理論變異函數模型實踐中，常用的是變異函數圖：偏基臺值:C塊金Thesillisthevalueatwhichthesemivariogramlevelsoff(itsasymptoticvalue)Therangeisthedistanceatwhichthesemivariogramlevelsoff(thespatialextentofstructureinthedata)Thenuggetisthesemivarianceatadistance0.0,(they–intercept)Asemivariogramisaplotofthestructurefunctionthat,likeautocorrelation,describestherelationshipbetweenmeasurementstakensomedistanceapart.Semivariogramsdefinetherangeordistanceoverwhichspatialdependenceexists.整理ppt54Thesillisthevalueatwhich整理ppt55整理ppt55ArcGIS：GeostatisticalAnalystprovidesthefollowingfunctionstochoosefromtomodeltheempiricalsemivariogram:CircularSphericalTetrasphericalPentasphericalExponentialGaussianRationalQuadraticHoleEffectK-BesselJ-BesselStable理論半方差模型的類型整理ppt56ArcGIS：GeostatisticalAnalyst例：假設某地區降水量Z(x)（單位：mm）是二維區域化隨機變量，滿足二階平穩假設，其觀測值的空間正方形網格數據如圖1所示（點與點之間的距離為h=1km）。試計算其南北方向及西北和東南方向的變異函數。整理ppt57例：整理ppt57圖1

空間正方形網格數據（點間距h=100m）

從圖1可以看出，空間上有些點，由于某種原因沒有采集到。如果沒有缺失值，可直接對正方形網格數據結構計算變異函數；在有缺失值的情況下，也可以計算變異函數。只要“跳過”缺失點位置即可（圖2）。整理ppt58圖1空間正方形網格數據（點間距h=100m）從圖1可以首先計算南北方向上的變異函數值，由變異函數的計算公式可得

=385/72=5.35圖2缺失值情況下樣本數對的組成和計算過程

☉為缺失值

整理ppt59首先計算南北方向上的變異函數值，由變異函數的計算公式可得整理ppt60整理ppt60

同樣計算出最后，得到南北方向和西北—東南方向上的變異函數（下表）。同樣可以計算東西方向上的變異函數。

方向：南北

方向：西北—東南

h1h2h3h4h5hh1.41h2.82h4.24h5.65h7.07hN(h)

362721135N(h)

322113825.359.2617.5525.6922.907.0612.9530.8558.1350.00整理ppt61同樣計算出方向：南北方向：西北—東南h1h2h3h4hsemivariogram南北向semivariogram西北—東南向整理ppt62semivariogramsemivariogram整理pp3.交叉驗證（CrossValidation）對每種插值方法重復下面的步驟，實現對不同插值方法的比較：從數據集中除去一個已知點的測量值；用剩余的點估計除去點的值；比較原始值和估計值，計算出估計值的預測誤差。針對每個已知點，進行上述步驟，然后評價不同插值方法的精確度。常用的評價指標是均方根（RMS）：整理ppt633.交叉驗證（CrossValidation）對每種插交叉驗證計算每個點的RMS計算某種插值方法的平均RMS?????123選擇某種插值方法45整理ppt64交叉驗證計算每個點的RMS計算某種插值方法的平均RMS???Kriging內插和外推背后的理論由法國數學家

GeorgesMatheron基于DanielG.Krige的碩士論文而發展，Krige是距離加權平均的開創者，用于當時南非的金礦勘探。TheEnglishverbistokrigeandthemostcommonnouniskriging.三、克里金插值方法克里金（kriging）插值：根據隨機場中待測值鄰居位置的觀測值對待測值進行插值的一組地統計技術。(e.g.,把

高程z作為地理位置的一個函數。)整理ppt65Kriging內插和外推背后的理論由法國數學家George克里金（Kriging）插值法，又稱空間局部估計或空間局部插值法，是地統計學的主要內容之一。克里金法建立在變異函數理論及結構分析基礎之上，它是在有限區域內對區域化變量進行無偏最優估計的一種方法。其實質是利用區域化變量的原始數據和變異函數的結構特點，對未采樣點的區域化變量值進行線性無偏最優估計。克里金法的適用條件：如果變異函數和相關分析的結果表明區域化變量存在空間相關性。整理ppt66克里金（Kriging）插值法，又稱空間局部估計或空間局部插假設： 所有的隨機誤差都具有二階平穩性，即，隨機誤差的均值為零，并且任意兩個隨機誤差之間的協方差只與兩者之間的距離和方向有關，而與它們的具體位置無關。克里金(Kriging)方法整理ppt67假設：克里金(Kriging)方法整理ppt67克里金方法的基本形式確定趨勢項空間相關隨機誤差項理解不同的克里金模型對誤差項的假設：期望值為0，并且和之間的自相關不取決于s點的位置，而取決于位移量h。為確保自相關方差有解，必須允許某兩點間的自相關可以相等。如，下面有箭頭相連的兩對位置點假設具有相同的自相關性。整理ppt68克里金方法的基本形式確定趨勢項空間相關隨機誤差項理解不同的克趨勢值可以被簡單地賦予一個常量，即，在任何位置處。理解不同的克里金模型如果在任何時候趨勢已知，無論趨勢是否是常量，都形成簡單克里金模型。趨勢也可以表示為：若趨勢中的系數未知，就是泛克里金模型。如果未知，就是普通克里金模型。整理ppt69趨勢值可以被簡單地賦予一個常量，即，在任何理解不同的克里金模型可以對Z(s)進行一下變換。例如，可以把它換為一個指示變量，若Z(s)低于一定的閾值（如空氣中的溴氧濃度值0.12ppm），則變為0；若Z(s)高于一定的閾值，則變為1。然后對高于閾值的情況進行預測。便形成了指示克里金模型。現在看一下方程的左邊。整理ppt70理解不同的克里金模型可以對Z(s)進行一下變換。例如，可以把普通克里金（OrdinaryKriging，OK）簡單克里金（SimpleKriging，SK）泛克里金（UniversalKriging，UK）指示克里金（IndicatorKriging，IK）整理ppt71普通克里金（OrdinaryKriging，OK）簡單克里協克里金（Cokriging，CK）假設和是未知常數。感興趣的變量是Z1,同時利用

Z1的自相關以及與Z2之間的交叉相關。整理ppt72協克里金（Cokriging，CK）假設和TypesofKrigingSimplekrigingisoptimalestimationofarandomfield,e.g.T(x),withaknownmean,m(x),andaknowncovariancePTT(x,x’).Ordinarykrigingisoptimalestimationofarandomfield,e.g.T(x),withanunknownconstantorlinearlytrendingmean,butaknownsemivariogramgTT(x,x’).Universalkrigingisoptimalestimationofarandomfield,e.g.T(x),withanunknownpolynomialtrendingmean,butaknownsemivariogramgTT(x,x’).整理ppt73TypesofKrigingSimplekriging運用克里金(Kriging)方法進行插值的過程可以分為兩步：第一步：變異函數分析。即，樣點的空間結構量化分析，是指對樣點數據擬合一個空間獨立模型。第二步：對未知點值進行預測。利用第一步擬合的變異函數、樣點數據的空間分布以及樣點數據值對未知點進行預測。整理ppt74運用克里金(Kriging)方法進行插值的過程可以分為兩步：普通克里格(OK)模型為：

其中，s=(x,y)為空間位置。如，觀測點(1,5)，Z(1,5)=100。未知均值整理ppt75普通克里格(OK)模型為： 其中，s=(x,y)為空間位置首先假設區域化變量滿足二階平穩假設，其數學期望為m，協方差函數及變異函數存在。即：假設在待估計點(x)的鄰域內共有n個實測點，即x1，x2，…，xn，其樣本值為Z(xi)。那么，普通克里格(OK)的插值公式為：第i個觀測點的權重預測點觀測點整理ppt76首先假設區域化變量滿足二階平穩假設，其數學期望為這個簡單的模型類似于倒距離權重法(IDW)，不同的是IDW的權重只與距預測點的距離有關，而普通克里格(OK)方法的權重則取決于半變異圖、距預測點的距離和預測點周圍觀測值的空間關系。x1x2x3x4x0InverseDistanceWeightsKrigingLocalMeans整理ppt77這個簡單的模型類似于倒距離權重法(IDW)，不同的是IDW的最優(Best)：要求盡可能小，從而預測值將盡可能的接近未知值線性(Linear)

:weightedlinearcombinationsofthedata無偏性(Unbiased)：一些點的預測值比真值高，一些比真值低，平均起來預測值和真值的差為0。要求權重系數值和為1。Estimation

普通克里金（OK）的目標：BLUE整理ppt78最優(Best)：要求在無偏估計的限制條件下，通過最小化可得到克里金方程：包含所有樣點對的半變異值基于點i和j之間距離的半變異值包含每一個觀測點與預測點之間的半變異值基于第i個觀測點與預測點之間距離的半變異值待定整理ppt79在無偏估計的限制條件下，通過最小化可得到克里金方程：包含所有普通克里金（OK）的計算過程計算樣點對之間的距離&方差，理論半變異值&平均理論半變異值擬合模型計算權系數與預測克里金方差整理ppt80普通克里金（OK）的計算過程計算樣點對之間的距離&方差，理論計算樣點對之間的距離和方差半變異值=0.5×方差示例：普通克里金方法半變異值計算結果位置距離方差(Zi-Zj)2半變異值距離計算整理ppt81計算樣點對之間的距離和方差半變異值=0.5×方差示例：普如果數據很大，樣點對的數目將迅速增加并且變得難以操作。因此，可以將樣點對分組，也就是將步長(h)分組。在該例中，首先將距離在1～2之間的作為第一組；h在2～3之間的作為第二組；依此類推（下表）。步長間距樣點對間距平均距離半變異值平均半變異值半變異值的步長分組表整理ppt82如果數據很大，樣點對的數目將迅速增加并且變得難以操作。因此，擬合模型 為了預測，需要用理論半方差模型擬合經驗半方差，以描述其變程、基臺值和塊金(range,sill,&nugget)。四個常用模型：整理ppt83擬合模型 為了預測，需要用理論半方差模型擬合經驗半方差，以描Gaussian:Linear:Spherical:Exponential:Where:c0=nuggetb=regressionslopea=rangec0+c=sillAssumesnosillorrange整理ppt84Gaussian:Linear:Spherical:Expo擬合模型 用平均半變異值為縱坐標，步長（抽樣間距）為橫坐標，做理論半變異圖。半變異值=斜率×距離

=13.5×h生成伽瑪矩陣：如，樣點(1,5)與樣點(3,4)的半變異值為：13.5×2.236=30.19整理ppt85擬合模型半變異值=斜率×距離生成伽瑪矩陣：如，樣點(1,伽瑪矩陣：最后一行的1和0根據無偏估計的限制條件求得。伽瑪矩陣的逆矩陣：整理ppt86伽瑪矩陣：最后一行的1和0根據無偏估計的限制條件求得。伽瑪矩計算權系數與預測 普通克里金的權系數矩陣為：向量g由未知點生成。如用點(1,4)來計算，計算該點與每個觀測點的距離，如(1,4)與(1,5)、(3,4)、(1,3)、(4,5)、(5,1)的距離。預測點(1,4)的g向量計算結果如下表。半變異值

=13.5×h整理ppt87計算權系數與預測向量g由未知點生成。如用點(1,4)來計算，計算權系數：權系數觀測值乘積克里金預測結果正如預期，權系數隨距離的增加而減小，但由于將各點的空間分布也考慮進去，其結果比直接的距離權重要好。整理ppt88計算權系數：權系數觀測值乘積克里金預測結果正如預期，權系數隨權系數為負值的處理DeutschClaytonV.Correctingfornegativeweightsinordinarykriging[J].Computers&Geosciences,1996,22(7):765-773.整理ppt89權系數為負值的處理DeutschClaytonV.Co克里格方差為了估計預測結果的不確定性，將權系數向量的每一行和g向量的每一行相乘，然后將結果相加，得到預測克里格方差，其平方根就是克里格標準差。g向量權系數乘積克里格方差克里格標準差如果假設誤差是正態分布的，則95%的置信區間為： 克里格預測值±1.96×克里格標準差

=[95.49,109.75]整理ppt90克里格方差為了估計預測結果的不確定性，將權系數向量的每一行和這僅是一個小例子，但從中可以看出幾個重要特點：Kriging是計算密集型的方法。需要一個合適的軟件。雖然一些GIS軟件提供半方差估計、建模和Kriging，該鄰域大多嚴肅的工作者使用特殊軟件如GSLIB、Variowin，或GS+。所有的結果依賴于為從樣點數據估計的半方差所擬合的模型，以及相關假設。其中包含一些主觀決定（多少個距離組合？擬合的模型是什么？基臺值和塊金值取多少？）。除了所討論的普通克里金，還有其他不同類型的克里金方法，比如，當存在均值漂移（drift）時采用UniversalKriging；CoKriging擴展到同時考慮兩個或更多變量。

整理ppt91這僅是一個小例子，但從中可以看出幾個重要特點：整理ppt91VariogramModelParametersWenowlookathowparametersofavariogram(covariance)modelaffectthe(OrdinaryKriging,OK)weightsSill,shape,nugget,range,andanisotropy整理ppt92VariogramModelParametersWenTheEffectofSillWithanyrescalingofthevariogram,neithertheKriging

weightsnortheestimatearechangedwhilethevarianceincreasesbythesamefactorusedtoscalethevariogram.整理ppt93TheEffectofSillWithanyresSillSillof10vs.20Sill=10Sill=20整理ppt94SillSillof10vs.20Sill=10TheEffectofShape高斯變異函數模型分配給近樣本點更大的權重。屏幕效應(screeneffect)–一個樣本位于另一個離未知點更近的樣本后面。具有小（或負的）權重，

如樣本5vs.6.高斯模型比指數模型具有更強的屏幕效應。Weightsthatarelessthan0orgreaterthan1canproduceestimateslargerthanthelargestsamplevalueorsmallerthanthesmallest.Weightswithin[0,1]produceestimatesonlywithintheminandmaxofsamplevaluesNegativeweightsmayproducenegativeestimates,althoughinmostscienceapplicationsvaluesshouldbepositive整理ppt95TheEffectofShape高斯變異函數模型分配給ShapeExponentialvs.GaussianmodelExponential

Gaussian整理ppt96ShapeExponentialExponentialGTheNuggetEffect塊金效應（nuggeteffect）使權重變得更相似，并導致更大的kriging方差。純塊金效應模型表明完全缺乏空間相關性，或在小于最小抽樣間距的尺度上才具有空間依賴性。整理ppt97TheNuggetEffect塊金效應（nuggeteNuggetNugget=0vs.=1/2sill整理ppt98NuggetNugget=0vs.=1/2silThe

EffectofRange變程減小，kriging方差增大。如果變程太小，那么所有樣本點離待估點的距離呈現相同遠近。于是估計類似于簡單平均，權重為1/n。Rangeofhvs.0.5hRange=10Range=20整理ppt99TheEffectofRange變程減小，kriginEffectofAnisotropy

MoreweightsaregiventosampleslieinthedirectionofmaximumcontinuityWeightsgiventothesamplesinthemaximumspatialcontinuitywouldincreaseastheanisotropyratiobecomeslarger整理ppt100EffectofAnisotropyMoreweigAnisotropyDirectionalvariogramsandcovariancefunctions整理ppt101AnisotropyDirectionalvariogrKriging在ArcGIS中的實現：探索性數據分析整理ppt102Kriging在ArcGIS中的實現：探索性數據分析整理ppKriging在ArcGIS中的實現：探索性數據分析整理ppt103Kriging在ArcGIS中的實現：探索性數據分析整理ppKriging在ArcGIS中的實現：探索性數據分析整理ppt104Kriging在ArcGIS中的實現：探索性數據分析整理ppKriging在ArcGIS中的實現：探索性數據分析整理ppt105Kriging在ArcGIS中的實現：探索性數據分析整理ppKriging在ArcGIS中的實現整理ppt106Kriging在ArcGIS中的實現整理ppt106Kriging在ArcGIS中的實現整理ppt107Kriging在ArcGIS中的實現整理ppt107Kriging在ArcGIS中的實現整理ppt108Kriging在ArcGIS中的實現整理ppt108Kriging在ArcGIS中的實現整理ppt109Kriging在ArcGIS中的實現整理ppt109Idaho州降雨量等直線圖Kriging在ArcGIS中的實現整理ppt110Idaho州降雨量等直線圖Kriging在ArcGIS中的實誤差圖Kriging在ArcGIS中的實現整理ppt111誤差圖Kriging在ArcGIS中的實現整理ppt111KrigingIDW插值趨勢面：5次多項式整理ppt112KrigingIDW插值趨勢面：5次多項式整理ppt112IDW插值普通克里格插值整理ppt113IDW插值普通克里格插值整理ppt113地形插值比較Voronoi圖TINIDWKriging樣條規則樣條整理ppt114地形插值比較Voronoi圖整理ppt114SimpleKrigingOdinaryKrigingUniversalKriging整理ppt115SimpleKrigingOdinaryKrigingUOdinaryKriging誤差圖UniversalKrigingSimpleKrigingSimpleKrigingOdinaryKrigingUniversalKriging整理ppt116OdinaryKriging誤差圖UniversalKrArcGIS：普通克里格的四張圖整理ppt117ArcGIS：普通克里格的四張圖整理ppt117整理ppt118整理ppt118參考資料O’SullivanD.andUnwinD.GeographicInformationAnalysis[M].JohnWiley&Sons,Inc.2003.p.45-49,p.265-281.劉湘南等.編著.GIS空間分析原理與方法（第二版）[M].科學出版社，2008.pp.195-209.王遠飛等編著.空間數據分析方法[M].科學出版社，2007.pp.152-175.ESRI：UsingArcGIS—GeostatisticalAnalysist.翻譯本：生態學空間分析原理與技術[M].科學出版社，2008.張仁鐸：空間變異理論與應用[M].科學出版社，2005.1.deSmitM.J.,etal.著.杜培軍等譯.地理空間分析——原理、技術與軟件工具[M].電子工業出版社，2008.整理ppt119參考資料O’SullivanD.andUnwinD.思考題IDW插值的基本思想是什么？以普通Kriging插值為例，闡述Kriging插值的基本思想。整理ppt120思考題IDW插值的基本思想是什么？整理ppt120Demo：packagegstat/

整理ppt121/整理ppMeusedataset155samplestakenonasupportof10x10mfromthetop0-20cmofalluvialsoilsina5x2kmpartthefloodplainoftheMaas(Meuse)nearStein(NL).id pointnumberx,y coordinatesEandN,inmcadmium concentrationinthesoil,inmgkg-1copper concentrationinthesoil,inmgkg-1lead concentrationinthesoil,inmgkg-1zinc concentrationinthesoil,inmgkg-1elev elevationabovelocalreferencelevel,inmom organicmatterlossonignition,inpercentffreq floodfrequencyclass, 1:annual,2:2-5years,3:every5yearssoil soilclass,codedlime hasthelandherebeenlimed?0or1=ForTlanduse landuse,codeddist.m distancefrommainRiverMaaschannel,inm整理ppt122Meusedataset155samplestakelibrary(gstat)data(meuse)summary(meuse)data(meuse.grid)g<-gstat(formula=log(zinc)~1,locations=~x+y,data=meuse,model=vgm(1,"Exp",300))x<-predict(g,meuse.grid)image(x,4,main="krigingvarianceanddatapoints")points(meuse$x,meuse$y,pch="+")coordinates(meuse)=~x+yplot(variogram(log(zinc)~1,meuse,cloud=TRUE))plot(variogram(log(zinc)~1,meuse),col=2)x<-variogram(log(zinc)~1,coordinates(meuse),data=meuse,cloud=TRUE)##commandsthatrequireinteraction:plot(plot(x,identify=TRUE),meuse)plot(plot(x,digitize=TRUE),meuse)vgm1<-variogram(log(zinc)~1,meuse)plot(vgm1)model.1<-fit.variogram(vgm1,vgm(1,"Sph",300,1))plot(vgm1,model=model.1)plot(vgm1,plot.numbers=TRUE,pch="+")vgm2<-variogram(log(zinc)~1,meuse,alpha=c(0,45,90,135))plot(vgm2)model.2<-vgm(.59,"Sph",926,.06,anis=c(0,0.3))#thefollowingdemonstratesplottingofdirectionalmodels:plot(vgm2,model=model.2)gridded(meuse.grid)=~x+ym<-vgm(.59,"Sph",874,.04)x<-krige(log(zinc)~1,meuse,meuse.grid,model=m)#ordinarykrigingspplot(x["var1.pred"],main="ordinarykrigingpredictions")spplot(x["var1.var"],main="ordinarykrigingvariance")x<-krige(log(zinc)~1,meuse,meuse.grid,model=m,beta=5.9)#simplekrigingm<-vgm(.4,"Sph",954,.06)#residualvariogramx<-krige(log(zinc)~x+y,meuse,meuse.grid,model=m,block=c(40,40))#universalblockkrigingspplot(x["var1.pred"],main="universalkrigingpredictions")m<-vgm(.59,"Sph",874,.04)x<-krige.cv(log(zinc)~1,meuse,m,nmax=40,nfold=5)#five-foldcrossvalidationbubble(x,"residual",main="log(zinc):5-foldCVresiduals")#demo(cokriging)整理ppt123library(gstat)整理ppt123第八講空間插值與地統計一階方法倒距離權重(IDW)插值趨勢面分析二階方法區域化變量協方差函數,變異函數,變異函數模型交叉驗證地統計：克里金(Kriging)方法概述：理解不同的克里金方法幾種不同的克里金方法

普通克里金(OK)、簡單克里金(SK)、泛克里金(UK)

指示克里金（IK）、協克里金（Ck）整理ppt124第八講空間插值與地統計一階方法整理ppt1一、一階方法倒距離權重（IDW）插值趨勢面分析整理ppt125一、一階方法倒距離權重（IDW）插值整理ppt2IDW插值方法假定每個輸入點都有著局部影響，這種影響隨著距離的增加而減弱。步驟:計算未知點到所有點的距離；計算每個點的權重：

權重是距離倒數的函數。計算結果：

1.倒距離權重(IDW)插值整理pptzhaoy@

zhaoy@IDW插值方法假定每個輸入點都有著局部影響，這種影響隨∑

=14

最大6

最小示例：IDW插值（求圖中0點的值）整理pptzhaoy@

zhaoy@∑=1示例：IDW插值（求圖中0點的值）整理ppt整理pptzhaoy@

zhaoy@整理pptzhaoy@zhIDW插值的一般模型

對所有或選定的i進行計算典型地，α

的取值是1或2除以距離權重的和，保證了權重加起來等于1整理pptzhaoy@

zhaoy@IDW插值的一般模型整理pptzhaoy@henu.e432100160200例：倒距離加權(IDW)插值結果整理pptzhaoy@

zhaoy@432100160200例：倒距離加權(IDW)插值結果整理IDW插值的缺點IDW不能得到大于樣本最大值或小于樣本最小值的估計。對于高程表面，這‘抹平’了峰和谷。（除非它們的高點和低點是樣本的一部分。因為估計值為均值，得到的表面將不通過樣本點整理ppt131IDW插值的缺點IDW不能得到大于樣本最大值或小于樣本最小值IDW插值在ArcGIS中的實現整理ppt132IDW插值在ArcGIS中的實現整理ppt9IDW插值在ArcGIS中的實現整理ppt133IDW插值在ArcGIS中的實現整理ppt10Idaho州降雨量等直線圖IDW插值在ArcGIS中的實現整理ppt134Idaho州降雨量等直線圖IDW插值在ArcGIS中的實現整2.趨勢面分析趨勢面分析(trendsurfaceanalysis)：用數學模型來模擬（或擬合）地理數據的空間分布及其區域性變化趨勢的方法。整理ppt1352.趨勢面分析趨勢面分析(trendsurfacea確定性插值函數擬合最普通的形式：多項式e.gy=ax2+bx+c局部（local）：分段地全局（global）local1storderglobal1storder整理ppt136確定性插值函數擬合local1storderglobal趨勢面分析是一種整體內插法，該方法假設一般趨勢與局部變化無關，并利用曲面方程來模擬待定點附近地形表面的一般趨勢。通常使用的是1次、2次、3次趨勢面，過高次的趨勢面不利于反映空間趨勢，并可能存在趨勢面的“畸變”。其中，2次趨勢面可用待定點附近的6個數據點來計算方程式系數。整理ppt137趨勢面分析是一種整體內插法，該方法假設一般趨勢與局部變化無關DeterministicSolutionsFirstOrderPolynomialInterpolationPredictedModelMeasuredSecondOrder(third,fourth,etc.)

PolynomialInterpolationLocalPolynomialInterpolationRadialBasisFunction(Spline)Interpolation整理ppt138DeterministicSolutionsFirstO趨勢面的性質與特點是一種光滑的數學曲面，能集中地代表地理數據在大范圍內的空間分布變化趨勢。與實際地理曲面不同，它只是實際曲面的一種近似。實際曲面包括趨勢面和剩余曲面兩部分，即：實際曲面＝趨勢面+剩余曲面整理ppt139趨勢面的性質與特點整理ppt16設Zi(xi，yi)表示某一地理特征值在空間上的分布。其中(xi，yi)為平面上點的坐標。任一觀測點Zi可分解為兩個部分，即：趨勢面剩余面整理ppt140設Zi(xi，yi)表示某一地理特征值在空間上的分布。其中(一階趨勢面FirstOrderTrendSurface一階趨勢面殘差ResidualsfromFirstOrderTrendSurface整理ppt141一階趨勢面一階趨勢面殘差整理ppt18趨勢面參數的確定（最小二乘法）使每一個觀測值與趨勢值的殘差平方和最小，即按建立多元線性方程的方法，使Q對系數b0，b1，…，bn求偏導，并令這些偏導數等于零，得趨勢面的正規方程組，解正規方程組，即可求出系數，從而得到趨勢面方程。整理ppt142趨勢面參數的確定（最小二乘法）整理ppt19因為任何函數在一定范圍內總可以用多項式來逼近，并可調整多項式的次數來滿足趨勢面分析的需要，一般來說，多項式的次數越高則趨勢值越接近于觀測值，而剩余值越小。整理ppt143因為任何函數在一定范圍內總可以用多項式來逼近，并可調整多項式多項式趨勢面的數學模型整理ppt144多項式趨勢面的數學模型整理ppt21整理ppt145整理ppt22二階趨勢面SecondOrderTrendSurface(1,x,y,x2,y2,xy)T三階趨勢面ThirdOrderTrendSurface

(1,x,y,x2,y2,xy,x3,y3,x2y,xy2)T整理ppt146二階趨勢面三階趨勢面整理ppt23趨勢面的具體計算方法與步驟：原始數據列表；等間隔選取縱橫坐標網，將原始數據點入坐標；按多元線性回歸分析方法求出趨勢面的正規方程組，解出參數；從趨勢值等值線圖中，獲得地理要素的區域性變化規律；用F分布對趨勢面進行擬合程度檢驗。整理ppt147趨勢面的具體計算方法與步驟：原始數據列表；整理ppt24即用雙三次多項式擬合趨勢面。雙三次多項式（樣條函數）插值整理ppt148即用雙三次多項式擬合趨勢面。雙三次多項式（樣條函數）插值整理該曲面模型有16個待定系數（Cij,i,j=0,1,2,3;）。通常用4個數據點（規則格網的4個頂點）的4個函數值組成的4×4方程組求解（如圖）。這4個函數值是高程Z、x方向斜率R、y方向斜率S，以及扭矩T：56781234916151011121314整理ppt149該曲面模型有16個待定系數（Cij,i,j=0,1,2,3其中Z保證曲面通過格網的4個數據點，R、S、T保證曲面在這4個數據點處光滑連續。雙三次多項式（樣條函數）內插法是規則格網插值的常用方法之一。這種方法通過一系列曲面片段來拼接地形表面，最終得到一個1階、2階連續的表面。該方法屬于局部插值，計算負擔中等；對于平滑表面擬合效果最好，對于起伏的表面擬合效果最差。整理ppt150其中Z保證曲面通過格網的4個數據點，R、S、T保證曲面在這4趨勢面分析在ArcGIS中的實現整理ppt151趨勢面分析在ArcGIS中的實現整理ppt28趨勢面分析在ArcGIS中的實現整理ppt152趨勢面分析在ArcGIS中的實現整理ppt29趨勢面分析在ArcGIS中的實現Idaho州降雨量等直線圖整理ppt153趨勢面分析在ArcGIS中的實現Idaho州降雨量等直線圖整二、二階方法IDW插值和趨勢面方法的缺陷：IDW：距離權重函數的選擇和“鄰居”的定義是其“致命缺陷”（Achilles’heel）。趨勢面分析：Inasense,trendsurfaceanalysisletsthedataspeakforthemselves,whereasIDWinterpolationforcesasetstructureontothem.整理ppt154二、二階方法IDW插值和趨勢面方法的缺陷：整理ppt31二、二階方法區域化變量協方差函數&

半變異函數整理ppt155二、二階方法區域化變量整理ppt32地統計學是以區域化變量理論為基礎，以變異函數為主要工具，研究那些在空間分布上既有隨機性又有結構性，或空間相關和依賴性現象的學科。協方差函數和變異函數是以區域化變量理論為基礎建立起來的地統計學的兩個最基本函數。地統計學的主要方法之一，克里金方法(Kriging)就是建立在變異函數理論和結構分析基礎之上的。整理ppt156地統計學是以區域化變量理論為基礎，以變異函數為主要工具，研究當一個變量的取值與其空間位置有關時，就稱為區域化變量（regionalizedvariable）。區域化變量常常反映某種空間現象的特征，用它來描述的現象稱之為區域化現象。區域化變量，亦稱區域化隨機變量，Matheron(1963)將它定義為以空間點x的三個直角坐標為自變量的隨機場區域化變量具有兩個最顯著，也是最重要的特征：隨機性和結構性。整理ppt157當一個變量的取值與其空間位置有關時，就稱為區域化變量（reg隨機變量隨機函數隨機過程隨機場區域化變量與時間有關的隨機函數帶有多個(2個以上)自變量的隨機函數以空間點的三個直角坐標為自變量整理ppt158隨機變量隨機函數隨機過程隨機場區域化變量與時間有關的隨機函數http://cg.ensmp.fr/Presentation/Matheron/Matheron_en.shtml

ProfessorGeorgesMatheron(1930-2000.8.7)法國數學家和地質學家

整理ppt159http://cg.ensmp.fr/Presentatio區域化變量的功能：

由于區域化變量是一種隨機函數，因而能同時反映空間變量的結構性和隨機性。一方面，當空間點x固定后，Z(x)就是一個隨機變量，這體現了其隨機性。另一方面，在空間兩個不同點x與x+h處的區域化變量值具有某種程度的相關性，這體現了其結構性。整理ppt160區域化變量的功能： 由于區域化變量是一種隨機函數，因而能同時區域化變量的組成部分

數據點結構性可以用均值和常數趨勢表示空間相關數據通常呈現正空間相關性

隨機性測量誤差，其他誤差

整理ppt161區域化變量的組成部分數據點結構性整理ppt38distanceelevation結構性隨機性實際值整理ppt162distanceelevation結構性隨機性實際值整協方差函數與變異函數數學期望、方差和協方差數學期望：一階原點矩方差：二階中心矩協方差：二階混合中心矩整理ppt163協方差函數與變異函數數學期望、方差和協方差整理ppt40協方差函數 類似地，當Z(x)是區域化變量時，對于任意兩點si和sj，空間隨機過程的協方差函數為： 相關系數和方差分別定義為：整理ppt164協方差函數整理ppt41整理ppt165整理ppt42若，則過程是二階平穩的，即均值與方差獨立于空間位置并在研究區域上是常數。于是有：

稱為協方差圖或過程的協方差函數，稱為相關圖或相關函數。顯然，協方差函數僅依賴于向量差h，當h=0時，整理ppt166若若獨立性僅是距離的函數，與方向無關，則空間過程是各向同性(isotropy)的。協方差函數就只依賴于距離向量h：

若在給定距離和方向上，不同位置數值差異的均值和方差為常數，則下式成立：半變異函數，也稱半方差函數。(semi-variogram)整理ppt167若獨立性僅是距離的函數，與方向無關，則空間過程是各向同性(i半變異函數(semi-variogram)—區域化變量的基本研究工具

—半變異函數就是區域化變量增量平方的數學期望之半區域化變量在i、i+h點的值步長為h的樣品對數步長(h)：在一定方向上，距離為h的矢量整理ppt168半變異函數(semi-variogram)區域化變量在i、i方差（Variance）變異函數（Variogram）h對比：方差&變異函數整理ppt169方差（Variance）h對比：方差&變異函數整理ppthaVetordistanceh建立經驗半變異函數(Semivariogram)測度空間變異。對于間隔距離為h的每一對Z(x)和Z(x+h)，測度它們之間差的平方。??h??h??h??h??h??h整理ppt170haVetordistanceh建立經驗半變異函數(Se半變異云圖半變異云圖：數據集中所有點對的“半變異函數—距離圖”。

如果空間相關性存在：

互相靠近的已知點趨向于具有較小的半方差

互相遠離的已知點趨向于具有較大的半方差整理ppt171半變異云圖半變異云圖：數據集中所有點對的“半變異函數—距離圖半變異云圖（semivariogramcloud）整理ppt172半變異云圖（semivariogramcloud）整理pp分組(Binning)如果數據很大，樣點對的數目將迅速增加并且變得難以操作。

因此，可以將樣點對分組，也就是步長(h)分組。如，可將距離在0～1之