克里金插值法的詳細(xì)介紹kriging

1、kriging插值作為地統(tǒng)計(jì)學(xué)中的一種插值方法由南非采礦工程師D.G.Krige于1951年首次提出，是 一種求最優(yōu)、線形、無(wú)偏的空間內(nèi)插方法。在充分考慮觀測(cè)資料之間的相互關(guān)系后，對(duì)每一個(gè)觀測(cè)資料賦 予一定的權(quán)重系數(shù)，加權(quán)平均得到估計(jì)值。這里介紹普通Kriging插值方法的基本步驟：1.該方法中衡量各點(diǎn)之間空間相關(guān)程度的測(cè)度是半方 差，其計(jì)算公式為：迥=土伊 3）-z（w）2h為各點(diǎn)之間距離，n是由h分開(kāi)的成對(duì)樣本點(diǎn)的數(shù)量，z是點(diǎn)的屬性值。2.在不同距離的半方差值都計(jì)算出來(lái)后，繪制半方差圖，橫軸代表距離，縱軸代表半方差。半方差圖 中有三個(gè)參數(shù)nugget （表示距離為零時(shí)的半方差）,sill

2、（表示基本達(dá)到恒定的半方差值），range（表示一個(gè) 值域范圍，在該范圍內(nèi)半方差隨距離增加，超過(guò)該范圍，半方差值趨于恒定）。利用做出的半方差圖找出與 之?dāng)M合的最好的理論變異函數(shù)模型（這是關(guān)鍵所在），可用于擬合的模型包括高斯模型、線性模型、球狀 模型、指數(shù)模型、圓形模型。00 二0c+c7ha -球狀模型，球面模型空間相關(guān)隨距離的增長(zhǎng)逐漸衰減，當(dāng)距離大于球面半徑后，空間相關(guān)消失。3.用擬合的模型計(jì)算出三個(gè)參數(shù)。例如球狀模型中nugget為c0，range為a，sill為c。4.利用擬合的模型估算未知點(diǎn)的屬性值，方程為：，Z0為估計(jì)值，zx是已知點(diǎn)的值，wx為權(quán)重，s是用來(lái)估算未知點(diǎn)的已知點(diǎn)的數(shù)目

3、。假如用三個(gè)點(diǎn)來(lái)估算，則有泌o）壯0）%壯0）a _1 _屁）/（偵1（松），（奴）1/（柘 1）/（奴）/（為）1 1110這樣權(quán)重就可以求出，然后估算未知點(diǎn)。(上述內(nèi)容根據(jù)地理信息系統(tǒng)導(dǎo)論(Kang-tsung Chang著；陳健飛等譯，科學(xué)出版社，2003) 第十三章內(nèi)容進(jìn)行總結(jié)，除球狀模型公式外其余公式皆來(lái)自此書)下面是本人自己編寫的利用海洋中斷面上觀測(cè)站點(diǎn)的實(shí)測(cè)溫度值來(lái)估算未觀測(cè)處的溫度的Fortran程 序，利用距離未知點(diǎn)最近的五個(gè)觀測(cè)點(diǎn)來(lái)估算未知點(diǎn)的溫度，選用模型為球狀模型。do ii=1,nx if(tgrid(ii,1)=0.)thendo i=1,dsite(ii)!首先尋

4、找距離最近的五個(gè)已知點(diǎn)位置do j=1,nhif(d(mm(ii),j).ne.0.or.j=1)thenhmie(j)=d(mm(ii),j)-dgrid(i)elsehmie(j)=9999end ifhmid(j)=abs(hmie(j)end dodo j=1,nh do k=j,nhif(hmid(j)hmid(k)thenelsem1=hmid(j)hmid(j)=hmid(k)hmid(k)=m1end ifend doend dodo j=1,5do k=1,nhif(abs(hmie(k)=hmid(j)thenlocat(j)=kend ifend doend dodo j

5、=1,4do k=j+1,5if(locat(j)=locat(k)thendo i3=1,nhif(abs(hmie(i3)=abs(hmie(locat(j).and.i3.ne.locat(j)thenlocat(j)=i3exitend ifenddoendifenddoenddo!然后求各點(diǎn)間距離，并求半方差do j=1,5do k=1,5hij(j,k)=abs(d(mm(ii),locat(j)-d(mm(ii),locat(k)/1000.end doend dodo j=1,5hio(j)=sqrt(hmid(j)*2+(abs(latgrid(ii)-lonlat(mm(i

6、i),2)*llat)*2$ +(abs(longrid(ii)-lonlat(mm(ii),1)*(1.112e5*$cos(0.017*(latgrid(ii)+lonlat(mm(ii),2)/2)*2)/1000.end dodo j=1,5do k=1,5if(hij(j,k).eq.0.)thenrleft(j,k)=0.elserleft(j,k)=sill*(1.5*hij(j,k)/range-0.5*hij(j,k)*3/range*3)end ifif(hio(j).eq.0.)thenrrig(1,j)=0.elserrig(1,j)=sill*(1.5*hio(j)/range-0.5*hio(j)*3/range*3)end ifend doend do rrig(1,6)=1. rleft(6,6)=0.do j=1,5 rleft(6,j)=1. rleft(j,6)=1.end dotry=rleftcall brinv(rleft,nnn,lll,is,js) ty1=matmul(try,rleft)!求權(quán)重wq=matmul(rrig,rleft)!插值所有格點(diǎn)上t,sdo j=1,5tgrid(ii,i)=t