克里金法時一種用於空間插值的地學統計方程式設計客棧法。
克里金法用半變異測定空間要素，要素即自相關要素。
半變異公式為：
其中(h)是已知點xi和xj的半變異，***h***表示這兩個點之間的距離，z是屬性值。
假設不存在漂移，普通克里金程式設計客棧法重點考慮空間相關因素，並用擬合的半變異直接進行插值。
估算某測量點z值的通用方程為：
式中，z0是待估計值，zx是已知點x的值，wx是每個已知點關聯的權重，s是用於估計的已知點數目。
權重可以由一組矩陣方程得到。
此程式對半變異進行擬合時採用的時最簡單的正比例函式擬合
資料為csv格式
儲存格式如下：
第一行為第乙個點以此類推
最後一行是待求點座標，其中z為未知值，暫且假設為0
**如下：
import numpy as np
from math import*
from numpy.linalg import *
h_data=np.loadtxt(open('高程點資料.csv'),delimiter=",",skiprows=0)
print('原始資料如下(x,y,z):\n未知點高程初值設為0\n',h_data)
def dis(p1,p2):
a=pow((pow((p1[0]-p2[0]),2)+pow((p1[1]-p2[1]),2)),0.5)
return a
def rh(z1,z2):
r=1/2*pow((z1[2]-z2[2]),2)
return r
def proportional(x,y):
xx,xy=0,0
for i in range(len(x)):
xx+=pow(x[i],2)
xy+=x[i]*y[i]
k=xy/xx
return k
r=;pp=;p=;
for i in range(len(h_data)):
pp.append(h_data[i])
for i in range(len(pp)):
for j in range(len(pp)):
p.append(dis(pp[i],pp[j]))
r.append(rh(pp[i],pp[j]))
r=np.array(r).reshape(len(h_data),len(h_data))
r=np.delete(r,len(h_data)-1,axis =0)
r=np.delete(r,len(h_data)-1,axis =1)
h=np.array(p).reshape(len(h_data),len(h_data))
h=np.delete(h,len(h_data)-1,axis =0)
oh=h[:,len(h_data)-1]
h=np.delete(h,len(h_data)-1,axis =1)
hh=np.triu(h,0)
rr=np.triu(r,0)
r0=;h0=;
for i in程式設計客棧 range(len(h_data)-1):
for j in range(len(h_data)-1):
if hh[i][j] !=0:
a=h[i][j]
h0.append(a)
if rr[i][j] !=0:
a=rr[i][j]
r0.append(a)
k=proportional(h0,r0)
hnew=h*k
a2=np.ones((1,len(h_data)-1))
a1=np.ones((len(h_data)-1,1))
a1=np.r_[a1,[[0]]]
hnew=np.r_[hnew,a2]
hnew=np.c_[hnew,a1]
print('半方差聯立矩陣：\n',hnew)
oh=np.array(k*oh)
oh=np.r_[oh,[1]]
w=np.dot(inv(hnew),oh)
print('權陣運算結果：\n',w)
z0,s2=0,0
for i in range(len(h_data)-1):
z0=w[i]*h_data[i][2]+z0
s2=w[i]*oh[i]+s2
s2=s2+w[len(h_data)-1]
print('未知點高程值為：\n',z0)
print('半變異值為：\n',pow(s2,0.5))
input()
運算結果
python初學，為了完成作業寫了個小程式來幫助計算，因為初學知識有限，有很多地方寫的很複雜，可以優化的地方很多。 還望讀者諒解，歡迎斧正謝謝！
(美)張康聰 著;陳健飛等譯. 地理資訊系統導論（第三版）. 北京：清華大學出版社, 2009.04.
本文標題: python 普通克里金（kriging）法的實現
本文位址: /jiaoben/python/293863.html
普通克里金插值
最近因為專案需要，研究了下克里金插值演算法。在地質學中，克里金插值演算法是一種使用的空間屬性估計技術，克里金插值說到底是個回歸問題，且依據的因素只有兩個位置之間的距離。克里金插值演算法又分為很多中，比如普通克里金插值，簡單克里金插值等，不同的克里金插值演算法只是假設條件不同。下面以普通克里金為例來說...
克里金插值
由於用supermap objects 沒有解決插值範圍的問題 見本版帖子 求助！哪位大俠在用supermap objects，請教乙個插值區域的問題 改用arcgis engine來做，現在遇到同樣的問題。使用iinterpolationop的krige方法已經實現了插值，但範圍侷限於氣象站點的外...
克里金插值c程式 克里金插值方法的原理
克里金插值方法原理 步驟1 4用來說明半方差函式的構建，步驟5說明了 模型，即如何求取未知點的數值 半方差函式訓練樣本的獲取公式 2.構建散點圖 x軸 距離，y軸 半方差值 3.根據已有的函式擬合經驗半方差圖 arcgis中提供了五種函式 circular spherical exponential...