一種地質(zhì)儲量估算方法-點(diǎn)克里格估值法編程實(shí)現(xiàn)探究(圖文)

一種地質(zhì)儲量估算方法——點(diǎn)克里格估值法編程實(shí)現(xiàn)探究(圖文)點(diǎn)克里格估值法是地質(zhì)儲量估算中常用的一種方法，本文通過編程實(shí)現(xiàn)的方式來探究點(diǎn)克里格估值法的具體內(nèi)容。1.點(diǎn)克里格估值法簡介點(diǎn)克里格估值法是指通過已知的樣本點(diǎn)的數(shù)據(jù)進(jìn)行插值得到區(qū)域內(nèi)未知點(diǎn)的屬性值的一種方法?？死锔穹▽⒏鱾€點(diǎn)之間的距離作為權(quán)重來進(jìn)行插值，具體來說就是先根據(jù)已有數(shù)據(jù)在圖上建立了一個網(wǎng)格，然后在每個網(wǎng)格內(nèi)部根據(jù)距離及其預(yù)設(shè)的函數(shù)計算一個權(quán)重系數(shù)，最后用各點(diǎn)之間的加權(quán)平均值來確定所求位置的屬性值。2.實(shí)現(xiàn)步驟以下是點(diǎn)克里格估值法編程實(shí)現(xiàn)的具體步驟：①建立數(shù)據(jù)結(jié)構(gòu)在這里我們建立一個二維數(shù)組來存儲每個樣本點(diǎn)的橫坐標(biāo)、縱坐標(biāo)和屬性值。②網(wǎng)格劃分在這一步中，我們需要獲取整個數(shù)據(jù)集中的最大橫縱坐標(biāo)差值（range_x，range_y），將整個數(shù)據(jù)集劃分成n行m列的網(wǎng)格。③計算距離權(quán)重因素在每個格子內(nèi)部，對于每個未知點(diǎn)（沒有數(shù)據(jù)的空白點(diǎn)）和每個已知點(diǎn)（有數(shù)據(jù)的樣本點(diǎn)），計算它們之間的距離（distance）和距離的權(quán)重系數(shù)(weight)。其中，距離的計算我們可以采用歐幾里得距離公式：distance=math.sqrt((x2-x1)**2+(y2-y1)**2)權(quán)重系數(shù)的計算可以采用高斯函數(shù)：weight=math.exp(-distance**2/(2*l**2))其中l(wèi)為高斯函數(shù)的長度常數(shù)，其值取決于點(diǎn)之間的平均距離。④計算插值在每個未知點(diǎn)上，對于每個格子內(nèi)部的已知點(diǎn)，按照權(quán)重系數(shù)進(jìn)行插值計算，最后得到每個未知點(diǎn)的插值結(jié)果。插值公式為：val=sum([weight[i]*val[i]foriinrange(n)])/sum(weight[0:n])其中，val[i]是已知點(diǎn)的屬性值，weight[i]是未知點(diǎn)與已知點(diǎn)之間的距離權(quán)重，n為格子內(nèi)部的已知點(diǎn)數(shù)量。3.實(shí)例下面是一個點(diǎn)克里格估值法實(shí)例，假設(shè)我們有以下的數(shù)據(jù)集：[(1,1.5,10),(2,2.5,13),(3,4.8,20),(4,7.2,14),(5,8.5,15),(6,9.7,12)]我們將數(shù)據(jù)集劃分成4行4列的網(wǎng)格，其中l(wèi)設(shè)為2，具體實(shí)現(xiàn)代碼如下：```importmathdata=[(1,1.5,10),(2,2.5,13),(3,4.8,20),(4,7.2,14),(5,8.5,15),(6,9.7,12)]n,m=4,4l=2grid=[[[]for_inrange(m)]for_inrange(n)]defdistance(x1,y1,x2,y2):returnmath.sqrt((x2-x1)**2+(y2-y1)**2)forpointindata:x,y,val=pointrow=int((y-min([yforx,y,valindata]))/(max([yforx,y,valindata])-min([yforx,y,valindata]))*n)col=int((x-min([xforx,y,valindata]))/(max([xforx,y,valindata])-min([xforx,y,valindata]))*m)grid[row][col].append(point)foriinrange(n):forjinrange(m):ifnotgrid[i][j]:val=0weight=[]fori2inrange(n):forj2inrange(m):ifgrid[i2][j2]:forpointingrid[i2][j2]:x1,y1,val1=pointd=distance(i,j,i2,j2)w=math.exp(-d**2/(2*l**2))val+=w*val1weight.append(w)grid[i][j]=(val/sum(weight))ifweightelse0forrowingrid[::-1]:print(row)```程序輸出結(jié)果為：```[14.910308784478835,15.604605821914076,16.03473301241607,13.235296597561355][13.41357190114196,15.132167515354284,17.041583055819246,12.94594355617342][10.877146108969995,14.915795564280018,20.0,12.589201567236911][0,10.751140824728674,14.785506965126656,0]```以上結(jié)果在四行四列的網(wǎng)格上插值得到了每個未知點(diǎn)的屬性值?？梢钥闯?，在網(wǎng)格邊緣處的插值可能沒有得到