地學統(tǒng)計-第五章課件

第五章空間插值與克里格法華中農(nóng)業(yè)大學楊勇8/5/20231第五章空間插值與克里格法華中農(nóng)業(yè)大學楊勇7/31/20對自然各種屬性的測量只能得到有限樣本點的值，不可能對每個點都進行采樣，然而我們總是想知道未測點的值，因此，就需要根據(jù)實測得到的離散數(shù)據(jù)，對未知點進行預測。8/5/20232華中農(nóng)業(yè)大學資源與環(huán)境學院對自然各種屬性的測量只能得到有限樣本點的值，不可能對每個點都算法分析8/5/20233華中農(nóng)業(yè)大學資源與環(huán)境學院算法分析7/31/20233華中農(nóng)業(yè)大學資源與環(huán)境學院要想獲得上圖，就必須有各種空間分析手段支持，特別是各種空間插值方法。大多數(shù)空間插值（預測）方法，都可被看做是數(shù)據(jù)的加權平均數(shù)，因此，有下面的通用預測公式：如何分配權重是關鍵問題，多數(shù)空間插值方法只考慮到系統(tǒng)的或確定的變異，而沒有考慮到其誤差，地統(tǒng)計的克里格方法，不但能描述空間變異的分布特征，而且能夠表達預測誤差。8/5/20234華中農(nóng)業(yè)大學資源與環(huán)境學院要想獲得上圖，就必須有各種空間分析手段支持，特別是各種空間插5.1主要空間插值法簡介分類：確定性方法：基于實測數(shù)據(jù)的相似性程度或平滑程度，利用數(shù)學函數(shù)進行插值（如逆距離加權法）地統(tǒng)計方法：利用實測數(shù)據(jù)的統(tǒng)計特性來量化其空間自相關程度，生產(chǎn)插值面并評價預測的不確定性8/5/20235華中農(nóng)業(yè)大學資源與環(huán)境學院5.1主要空間插值法簡介分類：7/31/20235華中農(nóng)業(yè)5.1主要空間插值法簡介分類：整體插值法：利用整個實測數(shù)據(jù)集來預測局部插值法：在大面積的研究區(qū)域上選取較小的空間單元，利用預測點周圍的臨近樣點來進行預測8/5/20236華中農(nóng)業(yè)大學資源與環(huán)境學院5.1主要空間插值法簡介分類：7/31/20236華中農(nóng)業(yè)5.1主要空間插值法簡介分類：精確性插值：實測點的預測值等于采得的實測值不精確插值：實測點的預測值不等于采得的值8/5/20237華中農(nóng)業(yè)大學資源與環(huán)境學院5.1主要空間插值法簡介分類：7/31/20237華中農(nóng)業(yè)5.1.1空間整體插值法1、全局多項式插值法（趨勢面分析法）：即用數(shù)學公式表達感興趣區(qū)域上的一種漸變的趨勢。平面：曲面：多項式中的參數(shù)系數(shù)往往用最小二乘法求解。但該方法是不精確的插值方法，很少有實測點剛好在生產(chǎn)的插值面上，而是或高或低于插值面，高低數(shù)值相加，之和近似為0。8/5/20238華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.1空間整體插值法1、全局多項式插值法（趨勢面分析法5.1.1空間整體插值法全局多項式插值法的插值結果往往呈條帶狀（左圖），適合于描述那些呈明顯趨勢分布的屬性，不適合描述那些空間分布波動較大（較破碎，右圖）的自然屬性8/5/20239華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.1空間整體插值法全局多項式插值法的插值結果往往呈條5.1.1空間整體插值法2、變換函數(shù)插值法：根據(jù)一個或多個空間參量的經(jīng)驗方程進行整體空間插值，這種經(jīng)驗方程稱為變換函數(shù)。即用與被預測屬性相關的其他屬性建立回歸方程，進行空間預測：b0,b1,b2為回歸系數(shù)，p1，p2為獨立空間變量，z(x)為被預測屬性8/5/202310華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.1空間整體插值法2、變換函數(shù)插值法：根據(jù)一個或多個5.1.2空間局部插值法1、泰森多邊形插值：由一組連續(xù)多邊形組成，多邊形的邊界是由相鄰兩點直線的垂直平分線組成。特性：（1）每個多邊形內(nèi)僅包含一個離散數(shù)據(jù)點。（2）在多邊形內(nèi)的任一點k(x，y)同Pi(xi，yi)之間距離總小于它同其它離散點Pj(xj，yj)之間距離。（3）泰森多邊形的任意一個頂點必有三條邊與它連接，這些邊是相鄰三個泰森多邊形兩兩拼接的公共邊。（4）泰森多邊形的任意一個頂點周圍存在三個離散點，將其連成三角形后其外接圓的圓心即為該頂點，該三角形稱泰森三角形8/5/202311華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.2空間局部插值法1、泰森多邊形插值：由一組連續(xù)多邊5.1.2空間局部插值法各泰森多邊形內(nèi)的每一點屬性均由各多邊形內(nèi)的已知點確定，若求數(shù)據(jù)域內(nèi)任意一點數(shù)據(jù)屬性Z（xi，yi），則需首先判斷待求點所落入的多邊形，然后再由控制該多邊形的已知點Z（x，y）推算得到。8/5/202312華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.2空間局部插值法各泰森多邊形內(nèi)的每一點屬性均由各多5.1.2空間局部插值法2、三角測量插值法：將采樣點用直線與其相鄰點連接成三角形，三角形內(nèi)部包括任何樣點，形成一個包括多個傾斜三角板的多面體（TIN）未測點只可能在三角形內(nèi)或三角形邊線上，利用線性插值即可求得缺點是每個預測值只是根據(jù)三個實測值得到，且有時會產(chǎn)生突變現(xiàn)象8/5/202313華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.2空間局部插值法2、三角測量插值法：將采樣點用直線5.1.2空間局部插值法3、逆距離加權法（IDW）：利用被預測區(qū)域點周圍的實測值來預測未采樣點的值，實測點離預測點越近，則對插值的結果影響越大。其中p為實測值對預測值的影響級，若p=0，則每一個權重是一樣的，預測值是所有實測值的平均值，當p增加時，相距較遠的點的權重迅速減小，2最為常用。由于IDW方法只考慮距離進行權重分配，所以臨近實測點的貢獻往往很大，而造成空間分布的多點中心現(xiàn)象。8/5/202314華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.2空間局部插值法3、逆距離加權法（IDW）：利用被5.1.2空間局部插值法4、局部多項式插值法（移動內(nèi)插法）：多項式插值法將整個區(qū)域考慮成一個平面或曲面，而局部多項式插值法是在劃定的領域內(nèi)（窗口內(nèi)）用其中的實測數(shù)據(jù)來擬合不同次數(shù)的多項式。8/5/202315華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.2空間局部插值法4、局部多項式插值法（移動內(nèi)插法）5.1.2空間局部插值法5、簡單移動平均法：6、樣條插值法：8/5/202316華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.2空間局部插值法5、簡單移動平均法：6、樣條插值法5.1.2空間局部插值法7、克里格方法：和IDW一樣，也是一種局部估計的加權平均，但是它對各實測點權重的確定是通過半方差分析獲取的，可分為線性克里格法和非線性克里格法。（1）普通克里格（6）概率克里格（2）簡單克里格（7）貝葉斯克里格（3）泛克里格（8）普通協(xié)同克里格（4）指示克里格（5）析取克里格8/5/202317華中農(nóng)業(yè)大學資源與環(huán)境學院5.1.2空間局部插值法7、克里格方法：和IDW一樣，也是克里格法實質(zhì)上是利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結構特點，對未采樣點的區(qū)域化變量的取值進行線性無偏最優(yōu)估計的一種方法，從數(shù)學角度講就是一種對空間分布的數(shù)據(jù)求線性最優(yōu)無偏內(nèi)插估計量的一種方法。是根據(jù)待估樣點有限領域內(nèi)若干已測定的樣點數(shù)據(jù)，在考慮樣點形狀、大小和空間相互位置關系，它們與待估樣點相互空間位置關系，以及變異函數(shù)提供的結構信息之后，對該待估樣點進行的一種線性無偏最優(yōu)估計8/5/202318華中農(nóng)業(yè)大學資源與環(huán)境學院克里格法實質(zhì)上是利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結構特點5.2線性預測克里格法5.2.1普通克里格法：假定Z(x)是滿足本證假設的一個隨機過程，該隨機過程有n個觀測值z(xi)，要預測未采樣點x0處的值，則線性預測值Z*(x0)可以表示如下：（5.2.1）Kriging是在使預測無偏并有最小方差的基礎上，去確定最優(yōu)的權重值，滿足以下兩個條件：（1）無偏性條件（5.2.2）（2）最優(yōu)條件：（5.2.3）8/5/202319華中農(nóng)業(yè)大學資源與環(huán)境學院5.2線性預測克里格法5.2.1普通克里格法：假定Z(x普通克里格法利用式5.2.1取代式5.2.2的左邊部分，有：（5.2.4）根據(jù)本征假設，上式可以進一步表示為：要使上式成立，則必須滿足以下條件：8/5/202320華中農(nóng)業(yè)大學資源與環(huán)境學院普通克里格法利用式5.2.1取代式5.2.2的左邊部分，有：本征假設條件1：條件2：r(h)稱為半方差函數(shù)，也叫變異函數(shù)本征假設是地統(tǒng)計學中對隨機函數(shù)的基本假設在研究區(qū)域內(nèi)，區(qū)域化變量Z(x)的增量的數(shù)學期望對任意x和h存在且等于0在研究區(qū)域內(nèi)，區(qū)域化變量的增量[Z(x)-Z(x+h)]的方差對任意x和h存在且平穩(wěn)8/5/202321華中農(nóng)業(yè)大學資源與環(huán)境學院本征假設條件1：r(h)稱為半方差函數(shù)，也叫變異函數(shù)在研究區(qū)普通克里格法在本證假設條件下，公式（5.2.3）左邊的式子可以表示為：根據(jù)方差最小原則，借助拉格朗日乘子，普通克里格的預測方程組為：預測方差為：8/5/202322華中農(nóng)業(yè)大學資源與環(huán)境學院普通克里格法在本證假設條件下，公式（5.2.3）左邊的式子可普通克里格法克里格公式也可以用矩陣的形式表示，對點狀克里格，有：8/5/202323華中農(nóng)業(yè)大學資源與環(huán)境學院普通克里格法克里格公式也可以用矩陣的形式表示，對點狀克里格，普通克里格法實例8/5/202324華中農(nóng)業(yè)大學資源與環(huán)境學院普通克里格法實例7/31/202324華中農(nóng)業(yè)大學資源與普通克里格法實例8/5/202325華中農(nóng)業(yè)大學資源與環(huán)境學院普通克里格法實例7/31/202325華中農(nóng)業(yè)大學資源與8/5/202326華中農(nóng)業(yè)大學資源與環(huán)境學院7/31/202326華中農(nóng)業(yè)大學資源與環(huán)境學院8/5/202327華中農(nóng)業(yè)大學資源與環(huán)境學院7/31/202327華中農(nóng)業(yè)大學資源與環(huán)境學院8/5/202328華中農(nóng)業(yè)大學資源與環(huán)境學院7/31/202328華中農(nóng)業(yè)大學資源與環(huán)境學院

？8/5/202329華中農(nóng)業(yè)大學資源與環(huán)境學院7/31/202329華中農(nóng)業(yè)大學資源與環(huán)境學院對于上圖，由于數(shù)據(jù)構型沒有變，因此，新的未知點對克里格矩陣沒有影響，對矩陣B有影響，因而值需要重新計算矩陣B即可，這樣可以節(jié)省計算過程，同時也是為什么特別強調(diào)野外實際取樣設計應盡可能采用規(guī)則格網(wǎng)數(shù)據(jù)結構的重要原因。8/5/202330華中農(nóng)業(yè)大學資源與環(huán)境學院對于上圖，由于數(shù)據(jù)構型沒有變，因此，新的未知點對克里格矩陣沒新的估計點x0的克里格方程組為：新的估計點的克里格估計值計算公式為：估計方差為：8/5/202331華中農(nóng)業(yè)大學資源與環(huán)境學院新的估計點x0的克里格方程組為：7/31/202331華中農(nóng)不規(guī)則網(wǎng)絡取樣數(shù)據(jù)實例在研究區(qū)域內(nèi)，區(qū)域化變量Z(x)滿足二階平穩(wěn)和本征假設，其協(xié)方差函數(shù)和變異函數(shù)存在，且均為各向同性的指數(shù)模型：8/5/202332華中農(nóng)業(yè)大學資源與環(huán)境學院不規(guī)則網(wǎng)絡取樣數(shù)據(jù)實例在研究區(qū)域內(nèi)，區(qū)域化變量Z(x)滿足二區(qū)域中位置點x0，在其領域內(nèi)有7個已知不規(guī)則分布的樣點，每一個樣點的坐標，樣點的觀測值與位置點x0的距離如表所示：8/5/202333華中農(nóng)業(yè)大學資源與環(huán)境學院區(qū)域中位置點x0，在其領域內(nèi)有7個已知不規(guī)則分布的樣點，每一7個已知樣點間的距離8/5/202334華中農(nóng)業(yè)大學資源與環(huán)境學院7個已知樣點間的距離7/31/202334華中農(nóng)業(yè)大學資源8/5/202335華中農(nóng)業(yè)大學資源與環(huán)境學院7/31/202335華中農(nóng)業(yè)大學資源與環(huán)境學院未知點x0的估計值為：其克里格估計方差為：8/5/202336華中農(nóng)業(yè)大學資源與環(huán)境學院未知點x0的估計值為：7/31/202336華中農(nóng)業(yè)大學資塊段普通克里格法（BlockKriging）如圖，如果估計的不是x0點，而是以x0為中心快段的平均值，這時要用塊段普通克里格法，由于與上例數(shù)據(jù)結構不變，因此，克里格矩陣不變，只是矩陣B變?yōu)?/5/202337華中農(nóng)業(yè)大學資源與環(huán)境學院塊段普通克里格法（BlockKriging）如圖，如果估計其中V為塊段面積，但積分經(jīng)常很難解，因此，常常（1）將待估塊段V離散化成若干點，計算已知樣點與塊段V內(nèi)被離散化的若干個點之間的協(xié)方差函數(shù)C(h)或變異函數(shù)r(h)，最后求出平均數(shù)；（2）除了將待估塊段V離散化外，如果已知樣點也是一個小塊段的支撐構成，將這塊段也離散化成一只的若干點，計算點與點之間的協(xié)方差函數(shù)C(h)或變異函數(shù)r(h)x0x01x02x03x048/5/202338華中農(nóng)業(yè)大學資源與環(huán)境學院其中x0x01x02x03x047/31/202338華中農(nóng)離散化后，待估塊段V中4個離散化點孩子見的平均協(xié)方差函數(shù)為：離散化是一個近似的方法，顯然，離散化越多，得到的平均協(xié)方差函數(shù)值或變異函數(shù)值精度越高，但計算的工作量很大，因此，在離散化時，要考慮精度和工作量的平衡。8/5/202339華中農(nóng)業(yè)大學資源與環(huán)境學院離散化后，待估塊段V中4個離散化點孩子見的平均協(xié)方差函數(shù)為：有關克里格法的說明在克里格方程組中，克里格矩陣只取決于已知樣點的相對集合特征，而與待估樣點或塊段無關，因此，只要兩個數(shù)據(jù)構型相同時，其克里格矩陣也相同。計算過程中，只需求一次逆矩陣，就可以得到另一個接的列矩陣，如果兩個待估塊段或樣點的幾個特征也相同，那么所得克里格權重系數(shù)的解得列矩陣也一定相同，鑒于這方面的原因，空間取樣設計時應盡量保證數(shù)據(jù)構形的系統(tǒng)性和規(guī)則性。這樣，就可將單獨一個克里格方案在整個估計中重復使用。8/5/202340華中農(nóng)業(yè)大學資源與環(huán)境學院有關克里格法的說明在克里格方程組中，克里格矩陣只取決于已知樣克里格估計的可靠性克里格估計量是一種線性、無偏最優(yōu)估計量，而且可以給出估計的方差，它在已知樣點上的估計等于樣點值本身。因而它是真正的空間局部內(nèi)插法，具有較高的可靠性。許多研究表明，克里格法的精度明顯高于多邊形法、三角形法、局部平均法和距離倒數(shù)平法加權法。8/5/202341華中農(nóng)業(yè)大學資源與環(huán)境學院克里格估計的可靠性7/31/202341華中農(nóng)業(yè)大學資源與克里格權重系數(shù)的對稱性如圖，每個樣點的權重系數(shù)分別為k1,k2…k7，當區(qū)域化變量Z(X)具有各向同性結構時，在圈中系數(shù)之間有k2=k3,k4=k6.這主要是因為樣點x2與x3，x4與x6對待估點x0幾何位置是對稱的，因而它們之間的克里格權重系數(shù)也具有對稱性。8/5/202342華中農(nóng)業(yè)大學資源與環(huán)境學院克里格權重系數(shù)的對稱性如圖，每個樣點的權重系數(shù)分別為k1,k叢聚效應降低克里格權重系數(shù)在上圖中，即使區(qū)域化變量Z(x)是各向同性結構的，x1與x5對待估樣點x0幾何對稱，也不能存在k1=k5，因為在樣點x5附近還存在樣點x4,x6和x7，這三個樣點與x5叢聚在一起，而x1是單獨的一個樣點，計算克里格權重系數(shù)的結果使k1>k5。x4,x5和x7的叢聚作用降低了x5對待估樣點x0的影響。而x1不存在叢聚效應。因此，在克里格估計中，不會由于一些樣點叢聚在一起而增大其權重系數(shù)，這也正是克里格法估計的優(yōu)點。8/5/202343華中農(nóng)業(yè)大學資源與環(huán)境學院叢聚效應降低克里格權重系數(shù)在上圖中，即使區(qū)域化變量Z(x)是叢聚效應降低克里格權重系數(shù)8/5/202344華中農(nóng)業(yè)大學資源與環(huán)境學院叢聚效應降低克里格權重系數(shù)7/31/202344華中農(nóng)業(yè)大學屏蔽效應當塊金值很小或不存在時，已知樣點的克里格權重系數(shù)的大小受屏蔽效應影響，如圖，已知樣點x5雖然與樣點x1到待估樣點x0的距離相等，但是x1的克里格權重系數(shù)k1確大于x5的權重系數(shù)k5，這主要是因為樣點x5受x4的屏蔽效應影響。8/5/202345華中農(nóng)業(yè)大學資源與環(huán)境學院屏蔽效應當塊金值很小或不存在時，已知樣點的克里格權重系數(shù)的大屏蔽效應如圖，待估點x0附近有x1,x2,……x12，共12個已知樣點，由數(shù)據(jù)構形可知，x1,x2,x3,x4與x0的幾何位置對稱且相等，x5,x6…x12與x0的幾何位置也對稱且相等。內(nèi)圈的每個已知樣點的克里格權重系數(shù)是(1-k)/4，外圈的克里格權重系數(shù)是k/8，總和為1，因此，內(nèi)圈的權重系數(shù)明顯大于外圈樣點的權重系數(shù)，這是由于內(nèi)圈樣點屏蔽了外圈樣點的緣故。8/5/202346華中農(nóng)業(yè)大學資源與環(huán)境學院屏蔽效應如圖，待估點x0附近有x1,x2,……x12，共12屏蔽效應8/5/202347華中農(nóng)業(yè)大學資源與環(huán)境學院屏蔽效應7/31/202347華中農(nóng)業(yè)大學資源與環(huán)境學院屏蔽效應屏蔽效應還與塊金常數(shù)有很大關系，當塊金常數(shù)增大時，屏蔽效應減弱，當為純塊金效應時，所有樣點之間相互獨立，協(xié)方差函數(shù)為0，變異函數(shù)等于外延方差，即基臺值，則待估樣點x0與周圍任何已知樣點的克里格權重系數(shù)均相同，此時屏蔽效應消失，任何一點上的克里格線性無偏最優(yōu)估計量都是所有樣點的算術平均值。8/5/202348華中農(nóng)業(yè)大學資源與環(huán)境學院屏蔽效應屏蔽效應還與塊金常數(shù)有很大關系，當塊金常數(shù)增大時，屏理論模型對克里格估值的影響數(shù)據(jù)由8個樣點和1個待估點組成，結構分析后，采用球狀模型作為理論模型樣點位置、屬性值及與待估點的距離8/5/202349華中農(nóng)業(yè)大學資源與環(huán)境學院理論模型對克里格估值的影響數(shù)據(jù)由8個樣點和1個待估點組成，結理論模型對克里格估值的影響球狀模型：8/5/202350華中農(nóng)業(yè)大學資源與環(huán)境學院理論模型對克里格估值的影響球狀模型：7/31/202350華1）尺度對克里格估值的影響若令r2(h)=0.5r1(h)，球狀模型尺度變化后，屏蔽效應減弱，方差較小，精度提高8/5/202351華中農(nóng)業(yè)大學資源與環(huán)境學院1）尺度對克里格估值的影響若令r2(h)=0.5r1(h)，（1）尺度對克里格估值的影響r1(h)r2(h)8/5/202352華中農(nóng)業(yè)大學資源與環(huán)境學院（1）尺度對克里格估值的影響r1(h)r2(h)7/31/2（2）變程對克里格估計的影響原球狀模型變程為40，如果變?yōu)?0，則新的球狀模型為，則內(nèi)圈樣點的克里格權重系數(shù)增大，而外圈樣點的克里格權重系數(shù)減小，增加了屏蔽效應，同時也增大了克里格估計方差，使估計的精度降低。8/5/202353華中農(nóng)業(yè)大學資源與環(huán)境學院（2）變程對克里格估計的影響原球狀模型變程為40，如果變?yōu)?（2）變程對克里格估計的影響8/5/202354華中農(nóng)業(yè)大學資源與環(huán)境學院（2）變程對克里格估計的影響7/31/202354華中農(nóng)業(yè)大（3）塊金效應對克里格估值的影響由變異函數(shù)性質(zhì)可知，，基臺值相當樣點之間自相關消失，相互獨立時樣本方差。對于一個已抽取的空間樣本，其方差已確定，因此基臺值=常數(shù)，當C0增加時，供高要降低，當C0增加到與樣點獨立時的方差時，C=0，此時基臺值就是塊金常數(shù)C0，因此，塊金常數(shù)C0的增大，導致樣點之間的相關性降低而獨立性加大。8/5/202355華中農(nóng)業(yè)大學資源與環(huán)境學院（3）塊金效應對克里格估值的影響由變異函數(shù)性質(zhì)可知，（3）塊金效應對克里格估值的影響若把原來的球狀模型塊金值擴大一倍，即C0=0.34，基臺值不變，則模型的曲線圖：8/5/202356華中農(nóng)業(yè)大學資源與環(huán)境學院（3）塊金效應對克里格估值的影響若把原來的球狀模型塊金值擴大（3）塊金效應對克里格估值的影響塊金值變大后，降低了內(nèi)圈樣點的權重系數(shù)，增大了外圈樣點的權重系數(shù)，而克里格方差也增大了近一倍，這說明塊金效應增大可使屏蔽效應降低，顯然，當C0=0.52時，8個樣點的克里格權重系數(shù)均為0.125，此時的塊金效應相當于純塊金效應。塊金值不變時塊金值增大1倍8/5/202357華中農(nóng)業(yè)大學資源與環(huán)境學院（3）塊金效應對克里格估值的影響塊金值變大后，降低了內(nèi)圈樣點（4）理論模型的種類對克里格估計的影響如果換成線性有基臺值模型：兩種模型變異函數(shù)曲線（左）和協(xié)方差函數(shù)曲線（右）對比8/5/202