克里金法時(shí)一種用于空間插值的地學(xué)統(tǒng)計(jì)方法。

成都創(chuàng)新互聯(lián)公司長(zhǎng)期為千余家客戶(hù)提供的網(wǎng)站建設(shè)服務(wù)，團(tuán)隊(duì)從業(yè)經(jīng)驗(yàn)10年，關(guān)注不同地域、不同群體，并針對(duì)不同對(duì)象提供差異化的產(chǎn)品和服務(wù)；打造開(kāi)放共贏平臺(tái)，與合作伙伴共同營(yíng)造健康的互聯(lián)網(wǎng)生態(tài)環(huán)境。為犍為企業(yè)提供專(zhuān)業(yè)的網(wǎng)站制作、做網(wǎng)站，犍為網(wǎng)站改版等技術(shù)服務(wù)。擁有10多年豐富建站經(jīng)驗(yàn)和眾多成功案例,為您定制開(kāi)發(fā)。

克里金法用半變異測(cè)定空間要素，要素即自相關(guān)要素。

半變異公式為：

其中γ(h) 是已知點(diǎn) xi 和 xj 的半變異，***h***表示這兩個(gè)點(diǎn)之間的距離，z是屬性值。

假設(shè)不存在漂移，普通克里金法重點(diǎn)考慮空間相關(guān)因素，并用擬合的半變異直接進(jìn)行插值。

估算某測(cè)量點(diǎn)z值的通用方程為：

式中，z0是待估計(jì)值，zx是已知點(diǎn)x的值，Wx是每個(gè)已知點(diǎn)關(guān)聯(lián)的權(quán)重，s是用于估計(jì)的已知點(diǎn)數(shù)目。
權(quán)重可以由一組矩陣方程得到。

此程序?qū)Π胱儺愡M(jìn)行擬合時(shí)采用的時(shí)最簡(jiǎn)單的正比例函數(shù)擬合

數(shù)據(jù)為csv格式

保存格式如下：

第一行為第一個(gè)點(diǎn)以此類(lèi)推

最后一行是待求點(diǎn)坐標(biāo)，其中z為未知值，暫且假設(shè)為0

代碼如下：

import numpy as np
from math import*
from numpy.linalg import *
h_data=np.loadtxt(open('高程點(diǎn)數(shù)據(jù).csv'),delimiter=",",skiprows=0)
print('原始數(shù)據(jù)如下(x,y,z):\n未知點(diǎn)高程初值設(shè)為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('半方差聯(lián)立矩陣：\n',hnew)
oh=np.array(k*oh)
oh=np.r_[oh,[1]]
w=np.dot(inv(hnew),oh)
print('權(quán)陣運(yùn)算結(jié)果：\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('未知點(diǎn)高程值為：\n',z0)
print('半變異值為：\n',pow(s2,0.5))
input()

另外有需要云服務(wù)器可以了解下創(chuàng)新互聯(lián)scvps.cn，海內(nèi)外云服務(wù)器15元起步，三天無(wú)理由+7*72小時(shí)售后在線(xiàn)，公司持有idc許可證，提供“云服務(wù)器、裸金屬服務(wù)器、高防服務(wù)器、香港服務(wù)器、美國(guó)服務(wù)器、虛擬主機(jī)、免備案服務(wù)器”等云主機(jī)租用服務(wù)以及企業(yè)上云的綜合解決方案，具有“安全穩(wěn)定、簡(jiǎn)單易用、服務(wù)可用性高、性?xún)r(jià)比高”等特點(diǎn)與優(yōu)勢(shì)，專(zhuān)為企業(yè)上云打造定制，能夠滿(mǎn)足用戶(hù)豐富、多元化的應(yīng)用場(chǎng)景需求。

標(biāo)題名稱(chēng)：python普通克里金（Kriging）法的實(shí)現(xiàn)-創(chuàng)新互聯(lián)
文章出自：http://weahome.cn/article/ceoijj.html