1、第八講 空間插值與地統(tǒng)計,一階方法 倒距離權重(IDW)插值 趨勢面分析 二階方法 區(qū)域化變量 協(xié)方差函數(shù), 變異函數(shù), 變異函數(shù)模型 交叉驗證 地統(tǒng)計：克里金(Kriging)方法 概述：理解不同的克里金方法 幾種不同的克里金方法 普通克里金(OK)、簡單克里金(SK)、泛克里金(UK) 指示克里金（IK）、協(xié)克里金（Ck）,一、一階方法,倒距離權重（IDW）插值 趨勢面分析,IDW 插值方法假定每個輸入點都有著局部影響，這種影響隨著距離的增加而減弱。步驟: 計算未知點到所有點的距離； 計算每個點的權重： 權重是距離倒數(shù)的函數(shù)。 計算結果：,1. 倒距離權重(IDW)插值, =1 4 最大

2、6 最小,示例：IDW插值（求圖中0點的值）,IDW插值的一般模型 對所有或選定的 i 進行計算 典型地， 的取值是1或2 除以距離權重的和，保證了權重加起來等于1,2020/9/22,河南大學環(huán)境與規(guī)劃學院 ,例：倒距離加權(IDW)插值結果,IDW插值的缺點,IDW不能得到大于樣本最大值或小于樣本最小值的估計。對于高程表面，這抹平了峰和谷。（除非它們的高點和低點是樣本的一部分。因為估計值為均值，得到的表面將不通過樣本點,IDW插值在ArcGIS中的實現(xiàn),IDW插值在ArcGIS中的實現(xiàn),Idaho州降雨量等直線圖,IDW插值在ArcGIS中的實現(xiàn),2. 趨勢面分析,趨勢面分析(trend

3、surface analysis)：用數(shù)學模型來模擬（或擬合）地理數(shù)據的空間分布及其區(qū)域性變化趨勢的方法。,確定性插值,函數(shù)擬合 最普通的形式：多項式 e.g y = ax2+bx+c 局部（local）：分段地 全局（global）,趨勢面分析是一種整體內插法，該方法假設一般趨勢與局部變化無關，并利用曲面方程來模擬待定點附近地形表面的一般趨勢。 通常使用的是1次、2次、3次趨勢面，過高次的趨勢面不利于反映空間趨勢，并可能存在趨勢面的“畸變”。其中，2次趨勢面可用待定點附近的6個數(shù)據點來計算方程式系數(shù)。,Deterministic Solutions,First Order Polynomia

4、l Interpolation,Predicted Model Measured,Second Order (third, fourth, etc.) Polynomial Interpolation,Local Polynomial Interpolation,Radial Basis Function (Spline) Interpolation,趨勢面的性質與特點 是一種光滑的數(shù)學曲面，能集中地代表地理數(shù)據在大范圍內的空間分布變化趨勢。 與實際地理曲面不同，它只是實際曲面的一種近似。 實際曲面包括趨勢面和剩余曲面兩部分，即： 實際曲面 趨勢面 + 剩余曲面,設Zi(xi，yi)表示某一地

5、理特征值在空間上的分布。其中(xi，yi)為平面上點的坐標。任一觀測點Zi可分解為兩個部分，即：,趨勢面,剩余面,一階趨勢面 First Order Trend Surface,一階趨勢面殘差 Residuals from First Order Trend Surface,趨勢面參數(shù)的確定（最小二乘法） 使每一個觀測值與趨勢值的殘差平方和最小，即 按建立多元線性方程的方法，使Q對系數(shù)b0，b1，bn求偏導，并令這些偏導數(shù)等于零，得趨勢面的正規(guī)方程組，解正規(guī)方程組，即可求出系數(shù)，從而得到趨勢面方程。,因為任何函數(shù)在一定范圍內總可以用多項式來逼近，并可調整多項式的次數(shù)來滿足趨勢面分析的需要，一般

6、來說，多項式的次數(shù)越高則趨勢值越接近于觀測值，而剩余值越小。,多項式趨勢面的數(shù)學模型,二階趨勢面 Second Order Trend Surface (1, x, y, x2, y2 , xy)T,三階趨勢面 Third Order Trend Surface (1, x, y, x2, y2 , xy, x3, y3 , x2y, xy2)T,趨勢面的具體計算方法與步驟：,原始數(shù)據列表； 等間隔選取縱橫坐標網，將原始數(shù)據點入坐標； 按多元線性回歸分析方法求出趨勢面的正規(guī)方程組，解出參數(shù)； 從趨勢值等值線圖中，獲得地理要素的區(qū)域性變化規(guī)律； 用F分布對趨勢面進行擬合程度檢驗。,即用雙三次多項

7、式擬合趨勢面。,雙三次多項式（樣條函數(shù)）插值,該曲面模型有16個待定系數(shù)（Cij, i,j=0,1,2,3;）。通常用4個數(shù)據點（規(guī)則格網的4個頂點）的4個函數(shù)值組成的44方程組求解（如圖）。這4個函數(shù)值是高程Z、 x方向斜率R、 y方向斜率S，以及扭矩T：,其中Z保證曲面通過格網的4個數(shù)據點，R、S、T保證曲面在這4個數(shù)據點處光滑連續(xù)。 雙三次多項式（樣條函數(shù)）內插法是規(guī)則格網插值的常用方法之一。這種方法通過一系列曲面片段來拼接地形表面，最終得到一個1階、2階連續(xù)的表面。該方法屬于局部插值，計算負擔中等；對于平滑表面擬合效果最好，對于起伏的表面擬合效果最差。,趨勢面分析在ArcGIS中的實現(xiàn)

8、,趨勢面分析在ArcGIS中的實現(xiàn),趨勢面分析在ArcGIS中的實現(xiàn),Idaho州降雨量等直線圖,二、二階方法,IDW插值和趨勢面方法的缺陷： IDW： 距離權重函數(shù)的選擇和“鄰居”的定義是其“致命缺陷”（Achilles heel）。 趨勢面分析：In a sense, trend surface analysis lets the data speak for themselves, whereas IDW interpolation forces a set structure onto them.,二、二階方法,區(qū)域化變量 協(xié)方差函數(shù) & 半變異函數(shù),地統(tǒng)計學是以區(qū)域化變量理論為基礎，

9、以變異函數(shù)為主要工具，研究那些在空間分布上既有隨機性又有結構性，或空間相關和依賴性現(xiàn)象的學科。 協(xié)方差函數(shù)和變異函數(shù)是以區(qū)域化變量理論為基礎建立起來的地統(tǒng)計學的兩個最基本函數(shù)。地統(tǒng)計學的主要方法之一，克里金方法(Kriging)就是建立在變異函數(shù)理論和結構分析基礎之上的。,當一個變量的取值與其空間位置有關時，就稱為區(qū)域化變量（regionalized variable）。區(qū)域化變量常常反映某種空間現(xiàn)象的特征，用它來描述的現(xiàn)象稱之為區(qū)域化現(xiàn)象。 區(qū)域化變量，亦稱區(qū)域化隨機變量，Matheron(1963)將它定義為以空間點x的三個直角坐標為自變量的隨機場 區(qū)域化變量具有兩個最顯著，也是最重要的特

10、征：隨機性和結構性。,隨機變量,隨機函數(shù),隨機過程,隨機場,區(qū)域化變量,與時間有關的隨機函數(shù),帶有多個(2個以上)自變量的隨機函數(shù),以空間點的三個直角坐標為自變量,http:/cg.ensmp.fr/Presentation/Matheron/Matheron_en.shtml,Professor Georges Matheron (1930-2000.8.7) 法國數(shù)學家和地質學家,區(qū)域化變量的功能：,由于區(qū)域化變量是一種隨機函數(shù)，因而能同時反映空間變量的結構性和隨機性。 一方面，當空間點 x 固定后，Z(x)就是一個隨機變量，這體現(xiàn)了其隨機性。 另一方面，在空間兩個不同點 x 與 x+h

11、處的區(qū)域化變量值具有某種程度的相關性，這體現(xiàn)了其結構性。,區(qū)域化變量的組成部分,數(shù)據點,結構性 可以用均值和常數(shù)趨勢表示 空間相關 數(shù)據通常呈現(xiàn)正空間相關性 隨機性 測量誤差，其他誤差,distance ,elevation,結構性,隨機性,實際值,協(xié)方差函數(shù)與變異函數(shù),數(shù)學期望、方差和協(xié)方差 數(shù)學期望：一階原點矩 方差：二階中心矩 協(xié)方差：二階混合中心矩,協(xié)方差函數(shù) 類似地，當Z(x)是區(qū)域化變量時，對于任意兩點si和sj ，空間隨機過程的協(xié)方差函數(shù)為： 相關系數(shù)和方差分別定義為：,若 ，則過程是二階平穩(wěn)的，即均值與方差獨立于空間位置并在研究區(qū)域上是常數(shù)。于是有： 稱為協(xié)方差圖或過程的協(xié)方差

12、函數(shù)， 稱為相關圖或相關函數(shù)。顯然，協(xié)方差函數(shù)僅依賴于向量差h，當h=0時，,若獨立性僅是距離的函數(shù)，與方向無關，則空間過程是各向同性(isotropy)的。協(xié)方差函數(shù)就只依賴于距離向量h：,若在給定距離和方向上，不同位置數(shù)值差異的均值和方差為常數(shù)，則下式成立：,半變異函數(shù)，也稱半方差函數(shù)。 (semi-variogram),半變異函數(shù)(semi-variogram) 區(qū)域化變量的基本研究工具 半變異函數(shù)就是區(qū)域化變量增量平方的 數(shù)學期望之半,區(qū)域化變量在i、i+h點的值,步長為h的樣品對數(shù),步長(h)：在一定方向上，距離為h的矢量,方差（Variance） 變異函數(shù)（Variogram）,對

13、比：方差 & 變異函數(shù),h,a,Vetor distance h,建立經驗半變異函數(shù)(Semivariogram),測度空間變異。 對于間隔距離為h的每一對Z(x) 和Z(x+h)，測度它們之間差的平方。,半變異云圖,半變異云圖：數(shù)據集中所有點對的“半變異函數(shù)距離圖”。 如果空間相關性存在： 互相靠近的已知點趨向于具有較小的半方差 互相遠離的已知點趨向于具有較大的半方差,半變異云圖（semivariogram cloud）,分組(Binning),如果數(shù)據很大，樣點對的數(shù)目將迅速增加并且變得難以操作。 因此，可以將樣點對分組，也就是步長(h)分組。如，可將距離在01之間的作為第一組；h在12之

14、間的作為第二組；依此類推。 Binning is a process that averages semivariance data by distance and direction (hey, its weighted!),Binned semivariogram,模型擬合,半變異函數(shù)必須擬合一個數(shù)學函數(shù)或模型，以用于估計任何給定距離的半方差。,理論變異函數(shù)模型,實踐中，常用的是變異函數(shù)圖：,not related anymore,變程范圍內才有結構性變化（有規(guī)律的變化）,反映隨機性大?。?主要來源于區(qū)域化變量Z(x)在小于抽樣尺度h時所具有的內部變異；另外還有抽樣分析誤差。,變異函數(shù)是一

15、個單調不減函數(shù)。當h超過某一個范圍，例如變程，變異函數(shù)不再增大，而是趨于一個極限值，即為基臺值。實際上等于區(qū)域化變量的先驗方差。即，,即基臺值與塊金值之差，表示數(shù)據中存在空間相關性引起的方差變化范圍。,The sill is the value at which the semivariogram levels off (its asymptotic value),The range is the distance at which the semivariogram levels off (the spatial extent of structure in the data),The nu

16、gget is the semivariance at a distance 0.0, (the y intercept),A semivariogram is a plot of the structure function that, like autocorrelation, describes the relationship between measurements taken some distance apart. Semivariograms define the range or distance over which spatial dependence exists.,A

17、rcGIS：Geostatistical Analyst provides the following functions to choose from to model the empirical semivariogram: Circular Spherical Tetraspherical Pentaspherical Exponential Gaussian Rational Quadratic Hole Effect K-Bessel J-Bessel Stable,理論半方差模型的類型,例： 假設某地區(qū)降水量Z(x)（單位：mm）是二維區(qū)域化隨機變量，滿足二階平穩(wěn)假設，其觀測值的空

18、間正方形網格數(shù)據如圖1所示（點與點之間的距離為h=1 km）。試計算其南北方向及西北和東南方向的變異函數(shù)。,圖1 空間正方形網格數(shù)據（點間距h=100m）,從圖1可以看出，空間上有些點，由于某種原因沒有采集到。如果沒有缺失值，可直接對正方形網格數(shù)據結構計算變異函數(shù)；在有缺失值的情況下，也可以計算變異函數(shù)。只要“跳過”缺失點位置即可（圖2）。,首先計算南北方向上的變異函數(shù)值，由變異函數(shù)的計算公式可得 =385/72=5.35,圖2 缺失值情況下樣本數(shù)對的組成和計算過程為缺失值,同樣計算出 最后，得到南北方向和西北東南方向上的變異函數(shù)（下表）。 同樣可以計算東西方向上的變異函數(shù)。,semivari

19、ogram 南北向,semivariogram 西北東南向,3. 交叉驗證（Cross Validation）,對每種插值方法重復下面的步驟，實現(xiàn)對不同插值方法的比較： 從數(shù)據集中除去一個已知點的測量值； 用剩余的點估計除去點的值； 比較原始值和估計值，計算出估計值的預測誤差。 針對每個已知點，進行上述步驟，然后評價不同插值方法的精確度。常用的評價指標是均方根（RMS）：,交叉驗證,計算每個點的RMS,計算某種插值方法的平均RMS,?,?,?,?,?,1,2,3,選擇某種插值方法,4,5,Kriging內插和外推背后的理論由法國數(shù)學家 Georges Matheron基于Daniel G. K

20、rige的碩士論文而發(fā)展， Krige是距離加權平均的開創(chuàng)者，用于當時南非的金礦勘探。 The English verb is to krige and the most common noun is kriging.,三、克里金插值方法,克里金（kriging）插值：根據隨機場中待測值鄰居位置的觀測值對待測值進行插值的一組地統(tǒng)計技術。 (e.g.,把 高程 z 作為地理位置的一個函數(shù)。),克里金（Kriging）插值法，又稱空間局部估計或空間局部插值法，是地統(tǒng)計學的主要內容之一?？死锝鸱ń⒃谧儺惡瘮?shù)理論及結構分析基礎之上，它是在有限區(qū)域內對區(qū)域化變量進行無偏最優(yōu)估計的一種方法。 其實質是利

21、用區(qū)域化變量的原始數(shù)據和變異函數(shù)的結構特點，對未采樣點的區(qū)域化變量值進行線性無偏最優(yōu)估計。 克里金法的適用條件：如果變異函數(shù)和相關分析的結果表明區(qū)域化變量存在空間相關性。,假設： 所有的隨機誤差都具有二階平穩(wěn)性，即，隨機誤差的均值為零，并且任意兩個隨機誤差之間的協(xié)方差只與兩者之間的距離和方向有關，而與它們的具體位置無關。,克里金(Kriging)方法,克里金方法的基本形式,確定趨勢項,空間相關隨機誤差項,理解不同的克里金模型,對誤差項的假設：期望值為0，并且 和 之間的自相關不取決于s點的位置，而取決于位移量h。為確保自相關方差有解，必須允許某兩點間的自相關可以相等。如，下面有箭頭相連的兩對位

22、置點假設具有相同的自相關性。,趨勢值 可以被簡單地賦予一個常量，即，在任何位置處 。,理解不同的克里金模型,如果在任何時候趨勢 已知，無論趨勢是否是常量，都形成簡單克里金模型。,趨勢也可以表示為： 若趨勢中的系數(shù)未知，就是泛克里金模型。,如果 未知，就是普通克里金模型。,理解不同的克里金模型,可以對Z(s)進行一下變換。例如，可以把它換為一個指示變量，若Z(s)低于一定的閾值（如空氣中的溴氧濃度值0.12ppm），則變?yōu)?；若Z(s)高于一定的閾值，則變?yōu)?。然后對高于閾值的情況進行預測。便形成了指示克里金模型。,現(xiàn)在看一下方程 的左邊。,普通克里金（Ordinary Kriging，OK）,

23、簡單克里金（Simple Kriging，SK）,泛克里金（Universal Kriging，UK）,指示克里金（Indicator Kriging，IK）,協(xié)克里金（Cokriging，CK）,假設 和 是未知常數(shù)。,感興趣的變量是 Z1 , 同時利用 Z1的自相關以及與Z2之間的交叉相關。,Types of Kriging,Simple kriging is optimal estimation of a random field, e.g. T(x), with a known mean, m(x), and a known covariance PTT(x,x). Ordinary

24、kriging is optimal estimation of a random field, e.g. T(x), with an unknown constant or linearly trending mean, but a known semivariogram gTT(x,x). Universal kriging is optimal estimation of a random field, e.g. T(x), with an unknown polynomial trending mean, but a known semivariogram gTT(x,x).,運用克里

25、金(Kriging)方法進行插值的過程可以分為兩步： 第一步：變異函數(shù)分析。即，樣點的空間結構量化分析，是指對樣點數(shù)據擬合一個空間獨立模型。 第二步：對未知點值進行預測。利用第一步擬合的變異函數(shù)、樣點數(shù)據的空間分布以及樣點數(shù)據值對未知點進行預測。,普通克里格(OK)模型為：,其中，s=(x, y)為空間位置。如，觀測點(1, 5)，Z(1, 5)=100。,未知均值,首先假設區(qū)域化變量 滿足二階平穩(wěn)假設，其數(shù)學期望為m，協(xié)方差函數(shù) 及變異函數(shù) 存在。即： 假設在待估計點(x)的鄰域內共有n個實測點，即x1，x2，xn，其樣本值為Z(xi) 。那么，普通克里格(OK)的插值公式為：,第i個觀測點

26、的權重,預測點,觀測點,這個簡單的模型類似于倒距離權重法(IDW)，不同的是IDW的權重只與距預測點的距離有關，而普通克里格(OK)方法的權重則取決于半變異圖、距預測點的距離和預測點周圍觀測值的空間關系。,最優(yōu)(Best)：要求 盡可能小，從而預測值將盡可能的接近未知值 線性(Linear) : weighted linear combinations of the data 無偏性(Unbiased)：一些點的預測值比真值高，一些比真值低，平均起來預測值和真值的差為0。要求權重系數(shù)值和為1。 Estimation,普通克里金（OK）的目標：BLUE,在無偏估計的限制條件下，通過最小化可得到克

27、里金方程：,包含所有樣點對的半變異值,基于點i和j之間距離的半變異值,包含每一個觀測點與預測點之間的半變異值,基于第i個觀測點與預測點之間距離的半變異值,待定,普通克里金（OK）的計算過程,計算樣點對之間的距離&方差，理論半變異值&平均理論半變異值 擬合模型 計算權系數(shù)與預測 克里金方差,計算樣點對之間的距離和方差,半變異值 = 0.5方差,示例：普通克里金方法,如果數(shù)據很大，樣點對的數(shù)目將迅速增加并且變得難以操作。因此，可以將樣點對分組，也就是將步長(h)分組。在該例中，首先將距離在12之間的作為第一組；h在23之間的作為第二組；依此類推（下表）。,擬合模型,為了預測，需要用理論半方差模型擬

28、合經驗半方差，以描述其變程、基臺值和塊金(range, sill, & nugget )。 四個常用模型：,Gaussian:,Linear:,Spherical:,Exponential:,Where: c0 = nugget b = regression slope a = range c0+ c = sill,Assumes no sill or range,擬合模型 用平均半變異值為縱坐標，步長（抽樣間距）為橫坐標，做理論半變異圖。,半變異值 = 斜率距離 = 13.5h,生成伽瑪矩陣：如，樣點(1,5)與樣點(3,4)的半變異值為：13.52.236=30.19,伽瑪矩陣：,最后一行

29、的1和0根據無偏估計的限制條件求得。,伽瑪矩陣的逆矩陣：,計算權系數(shù)與預測 普通克里金的權系數(shù)矩陣為：,向量g由未知點生成。如用點(1,4)來計算，計算該點與每個觀測點的距離，如(1,4)與(1,5)、 (3,4)、 (1,3)、 (4,5)、 (5,1)的距離。預測點(1,4)的g向量計算結果如下表。,半變異值 = 13.5h,計算權系數(shù)：,正如預期，權系數(shù)隨距離的增加而減小，但由于將各點的空間分布也考慮進去，其結果比直接的距離權重要好。,權系數(shù)為負值的處理,Deutsch Clayton V. Correcting for negative weights in ordinary krig

30、ingJ. Computers & Geosciences, 1996, 22(7): 765-773.,克里格方差,為了估計預測結果的不確定性，將權系數(shù)向量的每一行和g向量的每一行相乘，然后將結果相加，得到預測克里格方差，其平方根就是克里格標準差。,如果假設誤差是正態(tài)分布的，則95%的置信區(qū)間為： 克里格預測值1.96克里格標準差 = 95.49, 109.75,這僅是一個小例子，但從中可以看出幾個重要特點： Kriging是計算密集型的方法。 需要一個合適的軟件。雖然一些GIS軟件提供半方差估計、建模和Kriging，該鄰域大多嚴肅的工作者使用特殊軟件如GSLIB、Variowin，或GS

31、+。 所有的結果依賴于為從樣點數(shù)據估計的半方差所擬合的模型，以及相關假設。其中包含一些主觀決定（多少個距離組合？擬合的模型是什么？基臺值和塊金值取多少？）。 除了所討論的普通克里金，還有其他不同類型的克里金方法，比如，當存在均值漂移（drift）時采用Universal Kriging；CoKriging擴展到同時考慮兩個或更多變量。 ,Variogram Model Parameters,We now look at how parameters of a variogram (covariance) model affect

32、 the (Ordinary Kriging, OK) weights Sill, shape, nugget, range, and anisotropy,The Effect of Sill,With any rescaling of the variogram, neither the Kriging weights nor the estimate are changed while the variance increases by the same factor used to scale the variogram.,Sill,Sill of 10 vs. 20,Sill = 1

33、0,Sill = 20,The Effect of Shape,高斯變異函數(shù)模型分配給近樣本點更大的權重。 屏幕效應(screen effect) 一個樣本位于另一個離未知點更近的樣本后面。具有小（或負的）權重， 如樣本 5 vs. 6. 高斯模型比指數(shù)模型具有更強的屏幕效應。 Weights that are less than 0 or greater than 1 can produce estimates larger than the largest sample value or smaller than the smallest. Weights within 0,1 produ

34、ce estimates only within the min and max of sample values Negative weights may produce negative estimates, although in most science applications values should be positive,Shape,Exponential vs. Gaussian model,Exponential,Gaussian,The Nugget Effect,塊金效應（nugget effect）使權重變得更相似，并導致更大的kriging方差。 純塊金效應模型表

35、明完全缺乏空間相關性，或在小于最小抽樣間距的尺度上才具有空間依賴性。,Nugget,Nugget = 0 vs. =1/2 sill,The Effect of Range,變程減小，kriging方差增大。 如果變程太小，那么所有樣本點離待估點的距離呈現(xiàn)相同遠近。于是估計類似于簡單平均，權重為1/n。,Range of h vs. 0.5h,Range = 10,Range = 20,Effect of Anisotropy,More weights are given to samples lie in the direction of maximum continuity Weights

36、 given to the samples in the maximum spatial continuity would increase as the anisotropy ratio becomes larger,Anisotropy,Directional variograms and covariance functions,Kriging在ArcGIS中的實現(xiàn)：探索性數(shù)據分析,Kriging在ArcGIS中的實現(xiàn)：探索性數(shù)據分析,Kriging在ArcGIS中的實現(xiàn)：探索性數(shù)據分析,Kriging在ArcGIS中的實現(xiàn)：探索性數(shù)據分析,Kriging在ArcGIS中的實現(xiàn),Krig

37、ing在ArcGIS中的實現(xiàn),Kriging在ArcGIS中的實現(xiàn),Kriging在ArcGIS中的實現(xiàn),Idaho州降雨量等直線圖,Kriging在ArcGIS中的實現(xiàn),誤差圖,Kriging在ArcGIS中的實現(xiàn),Kriging,IDW插值,趨勢面：5次多項式,IDW插值,普通克里格插值,地形插值比較,Voronoi圖 TIN IDW Kriging 樣條 規(guī)則樣條,Simple Kriging,Odinary Kriging,Universal Kriging,Odinary Kriging,誤差圖,Universal Kriging,Simple Kriging,Simple Krig

38、ing,Odinary Kriging,Universal Kriging,ArcGIS：普通克里格的四張圖,參考資料,OSullivan D. and Unwin D. Geographic Information Analysis M. John Wiley & Sons, Inc. 2003. p.45-49, p.265-281. 劉湘南等. 編著. GIS空間分析原理與方法（第二版）M. 科學出版社，2008. pp.195-209. 王遠飛等編著. 空間數(shù)據分析方法M. 科學出版社，2007. pp.152-175. ESRI：Using ArcGISGeostatistical

39、Analysist. 翻譯本：生態(tài)學空間分析原理與技術M. 科學出版社，2008. 張仁鐸：空間變異理論與應用M. 科學出版社，2005.1. de Smit M. J., et al. 著. 杜培軍等譯. 地理空間分析原理、技術與軟件工具M. 電子工業(yè)出版社，2008.,思考題,IDW插值的基本思想是什么？ 以普通Kriging插值為例，闡述Kriging插值的基本思想。,Demo： package gstat,/,Meuse data set,155 samples taken on a support of 10 x10 m from the

40、 top 0-20 cm of alluvial soils in a 5x2 km part the floodplain of the Maas (Meuse) near Stein (NL). id point number x, y coordinates E and N, in m cadmium concentration in the soil, in mg kg-1 copper concentration in the soil, in mg kg-1 lead concentration in the soil, in mg kg-1 zinc concentration

41、in the soil, in mg kg-1 elev elevation above local reference level, in m om organic matter loss on ignition, in percent ffreq flood frequency class, 1: annual, 2: 2-5 years, 3: every 5 years soil soil class, coded lime has the land here been limed? 0 or 1 = F or T landuse land use, coded dist.m distance from main River Maas channel, in m,library(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=kriging variance and data points