顧及空間各向異性的IDW插值算法

孫夢(mèng)楠，劉少華，劉京城

(長(zhǎng)江大學(xué) 地球科學(xué)學(xué)院，湖北 武漢 430100)

0 引 言

近年來(lái)，隨著科學(xué)技術(shù)的快速發(fā)展，對(duì)空間數(shù)據(jù)的處理有了更高精度的要求，空間插值方法的優(yōu)化和應(yīng)用越來(lái)越受到人們的重視。由于空間數(shù)據(jù)具有一定的空間依賴(lài)性與相關(guān)性，為了對(duì)待插值數(shù)據(jù)進(jìn)行最優(yōu)內(nèi)插計(jì)算，我們需要考慮數(shù)據(jù)空間分布的相關(guān)性和變異性，使得插值結(jié)果更趨于實(shí)際的要求。距離平方反比插值方法是眾多插值方法之一，是以已知點(diǎn)與待插值點(diǎn)之間的距離為基礎(chǔ)，在搜索鄰域內(nèi)對(duì)待插點(diǎn)進(jìn)行無(wú)偏最優(yōu)插值的一種方法[1]。

現(xiàn)實(shí)生活中，事物具有一定的相關(guān)性，同時(shí)也存在著各向異性。各向異性主要是指物體在不同方向上屬性值變化的速度具有一定的差異性，甚至屬性值有可能突然變化的情況[2,3]。各向異性是地質(zhì)地層研究過(guò)程中常見(jiàn)的性質(zhì)。例如，長(zhǎng)江黃河中河道方向與垂直方向水體污染物的濃度分布規(guī)律，大氣環(huán)境污染物分布受到風(fēng)向和重力的影響，礦產(chǎn)儲(chǔ)量分布受到地質(zhì)構(gòu)造環(huán)境的影響等等。在實(shí)際應(yīng)用中，空間插值必須考慮物體在三維空間中的空間結(jié)構(gòu)性和各向異性特征，才能構(gòu)建與實(shí)體相統(tǒng)一的模型。

本文基于傳統(tǒng)的距離平方反比插值法，參考克里格方法[4]，結(jié)合物體的各向異性特征，探討并建立了基于方向的各向異性插值優(yōu)化算法，通過(guò)對(duì)模型進(jìn)行不同的參數(shù)設(shè)置，對(duì)相關(guān)結(jié)果進(jìn)行比較，驗(yàn)證了算法在空間插值方面的有效性。

1 傳統(tǒng)的距離平方反比插值法

距離平方反比插值算法的主要原理是，根據(jù)空間各已知點(diǎn)的坐標(biāo) (Xi，Yi)(i=1，2，…，n)， 和各已知點(diǎn)屬性值Zi(i=1，2，…，n)， 從而假設(shè)待插值點(diǎn)的坐標(biāo)為 (x，y)， 插值得到的結(jié)果為Z(x，y)。Z(x，y) 可理解為各已知點(diǎn)屬性值與待求結(jié)果的距離平方加權(quán)平均，如式(1)所示

(1)

該方法的算法簡(jiǎn)單，計(jì)算量小，計(jì)算時(shí)消耗計(jì)算機(jī)的空間、花費(fèi)的時(shí)間相對(duì)較小，時(shí)間與空間復(fù)雜度低，但存在如下問(wèn)題[5]。

(1)只是一種純幾何的加權(quán)。若鄰域內(nèi)多個(gè)點(diǎn)到待插點(diǎn)的距離相同，則它們的權(quán)重系數(shù)相同，這樣就忽視了不同方位的已知點(diǎn)由于各向異性對(duì)預(yù)測(cè)點(diǎn)造成的影響，即沒(méi)有考慮空間數(shù)據(jù)的各向異性特征。

(2)搜索半徑r值以及樣本點(diǎn)的選取無(wú)標(biāo)準(zhǔn)。插值結(jié)果較大程度上決定于鄰近點(diǎn)的選取和搜索范圍的確定。

2 基于方向的各向異性插值算法

2.1 插值法中參數(shù)的確定

各向異性是一種常見(jiàn)的現(xiàn)象，是物質(zhì)的全部或部分化學(xué)、物理等性質(zhì)隨方向的改變而變化的特性[6]。考慮到待插值數(shù)據(jù)具有各向異性特征，優(yōu)化方法將搜索區(qū)域由原來(lái)的圓形改為橢圓，如圖1所示，使得在搜索區(qū)域內(nèi)的點(diǎn)更具有代表性，能夠?qū)Υ逯迭c(diǎn)有所影響。對(duì)比圖1(a)、圖1(b)兩個(gè)圖形，發(fā)現(xiàn)區(qū)域內(nèi)取得的建模點(diǎn)有所變化。

圖1 兩種方法的搜索區(qū)域?qū)Ρ?/p>

通過(guò)設(shè)置模型變量來(lái)求取橢圓公式，如圖2所示，其中a為橢圓長(zhǎng)軸的值，b為橢圓短軸的值，θ為橢圓長(zhǎng)軸延長(zhǎng)線與x軸的夾角。模型參數(shù)可通過(guò)最小二乘法進(jìn)行函數(shù)匹配來(lái)確定。由于最小二乘法橢圓擬合算法可能包含誤差較大的樣本點(diǎn)，該樣本點(diǎn)參與運(yùn)算會(huì)使橢圓擬合的結(jié)果產(chǎn)生偏差[7,8]。針對(duì)這種情況，采用基于最小二乘法的橢圓擬合改進(jìn)算法確定橢圓公式。隨機(jī)選取6個(gè)樣本點(diǎn)計(jì)算橢圓公式，統(tǒng)計(jì)與此公式匹配的所有樣本點(diǎn)個(gè)數(shù)，重復(fù)上述步驟一定次數(shù)。其中，匹配樣本點(diǎn)多的橢圓即為最優(yōu)橢圓，從而確定橢圓公式中參數(shù)的值[9,10]。

圖2 優(yōu)化模型變量三

一般方程式橢圓方程，如式(2)所示

Ax2+Bxy+Cy2+Dx+Ey+F=0

(2)

由于解的任意整數(shù)倍都是同一橢圓方程，且為了避免零值，設(shè)置約束條件如式(3)所示

A+C=1

(3)

通過(guò)求取目標(biāo)函數(shù)(4)的最小值來(lái)確定各系數(shù)

(4)

根據(jù)極值原理，欲使目標(biāo)函數(shù)值最小，必有式(5)

(5)

由此可得一個(gè)線性方程組，利用全主元高斯消去法并結(jié)合約束條件，求得方程系數(shù)A、B、C、D、E、F的值[11]。對(duì)橢圓方程進(jìn)行化解，得到橢圓的標(biāo)準(zhǔn)方程，同時(shí)確定參數(shù)θ、a、b的值。

根據(jù)橢圓一般方程式和參數(shù)θ、a、b的值，可計(jì)算得到旋轉(zhuǎn)角度為θ的標(biāo)準(zhǔn)斜橢圓方程，在計(jì)算對(duì)待插值有影響的點(diǎn)時(shí)，得到橢圓式(6)

(6)

如果樣本點(diǎn)代入式(6)，其值小于1，則為感興趣點(diǎn)，代入計(jì)算。

2.2 距離平方歸一化處理

如圖3所示，通過(guò)點(diǎn)O(X0，Y0)、P(X1，Y1)、Q(X，Y)， 求取直線OP公式，如式(7)所示

圖3 歸一化點(diǎn)位

(7)

為實(shí)現(xiàn)歸一化，需要求OP與OQ長(zhǎng)度的比值，即k=(x1-x0)/(x-x0) 的值[12,13]。結(jié)合橢圓公式，可得到以下方程組(8)

(8)

(9)

(10)

(11)

b2[(x1-x0)(x-x0)cosθ-(y1-y0)(x-x0)sinθ]2+
a2[(x1-x0)(x-x0)sinθ+(y1-y0)(x-x0)cosθ]2=
a2b2(x1-x0)2

(12)

b2{(x-x0)[(x1-x0)cosθ-(y1-y0)sinθ]}2+
a2{(x-x0)2[(x1-x0)sinθ+(y1-y0)cosθ]2}=
a2b2(x1-x0)2

(13)

(14)

(15)

推導(dǎo)公式如式(9)至式(15)所示，從而得到式(16)，取得k值如下

(16)

根據(jù)k值，將原方法中距離d值歸一化，得d′=k*b(b為橢圓短軸)。考慮到物體的各向異性特征，通過(guò)異向比壓縮距離軸，使之成為各向同性的模型，即構(gòu)成一個(gè)新的球狀模型，如圖4中虛線圓。將原方法各點(diǎn)距離平方權(quán)重改變，最后d′代入原方法d計(jì)算即可。

圖4 歸一化球狀模型

3 算法驗(yàn)證與分析

為驗(yàn)證本文算法的可行性與精確度，本文選取某油田油井分層數(shù)據(jù)C1層的30個(gè)點(diǎn)層位數(shù)據(jù)作為測(cè)試數(shù)據(jù)。

3.1 交叉驗(yàn)證

將測(cè)試數(shù)據(jù)劃分為插值點(diǎn)和檢查點(diǎn)，選取10個(gè)點(diǎn)作為檢查點(diǎn)，對(duì)距離平方反比插值法和改進(jìn)算法進(jìn)行交叉驗(yàn)證。實(shí)際值與兩種方法預(yù)測(cè)值對(duì)比如圖5、圖6所示，發(fā)現(xiàn)原方法中實(shí)際值與預(yù)測(cè)值繪制的點(diǎn)偏離直線較遠(yuǎn)，插值結(jié)果與實(shí)際值的離散程度以及插值誤差的離散程度均大于本文方法。以誤差標(biāo)準(zhǔn)差、平均誤差、中誤差和方差作為插值結(jié)果的評(píng)價(jià)標(biāo)準(zhǔn)，評(píng)價(jià)標(biāo)準(zhǔn)值越小，表明插值結(jié)果越準(zhǔn)確。數(shù)據(jù)結(jié)果見(jiàn)表1，后者的評(píng)價(jià)標(biāo)準(zhǔn)值均小于前者，其中距離平方反比插值法的實(shí)際值與預(yù)測(cè)值的方差較本文方法更大[14,15]。從而分析得出，本算法和傳統(tǒng)算法相比，精度相對(duì)較高，可行性較好，插值結(jié)果較準(zhǔn)確。

圖5 實(shí)際值與原方法預(yù)測(cè)值對(duì)比

圖6 實(shí)際值與本方法預(yù)測(cè)值對(duì)比

統(tǒng)計(jì)參數(shù)距離平方反比插值法改進(jìn)方法誤差標(biāo)準(zhǔn)差16.859 081 239.848 568 424平均誤差1597.4631594.18中誤差44.486 960 8542.719 062 8方差1979.089 6851824.918 326

3.2 對(duì)比驗(yàn)證

根據(jù)最小二乘法的橢圓擬合改進(jìn)算法確定橢圓參數(shù)，將橢圓旋轉(zhuǎn)不同的角度做空間插值計(jì)算，并與距離平方反比插值法進(jìn)行比較。圖7為結(jié)果對(duì)比，其中圖7(a)為距離平方反比插值法的結(jié)果，圖7(d)為本算法的結(jié)果，圖7(b)、圖7(c)分別是將斜橢圓公式的θ值設(shè)為0°、60°的結(jié)果。

圖7 不同方向插值結(jié)果對(duì)比

鑒于傳統(tǒng)方法存在的問(wèn)題，在進(jìn)一步研究本算法時(shí)，首先考慮如何最大化利用數(shù)據(jù)方位的各向異性來(lái)提高建模精度。驗(yàn)證方案提取的4種方式結(jié)果，從砂體地層面表面重構(gòu)效果來(lái)看，傳統(tǒng)方法和優(yōu)化的方案在模型表面特征方面，存在一定的相似性，但在細(xì)節(jié)上差異較大。

從宏觀上看，4種方案預(yù)測(cè)的砂體地層面表面模型特征基本一致，砂體地層面的表面模型與已知層位數(shù)據(jù)關(guān)系明顯，局部區(qū)域的最值基本都是已知樣本點(diǎn)所在的位置，極值分布規(guī)律清晰，與地質(zhì)、測(cè)井?dāng)?shù)據(jù)基本一致。在細(xì)節(jié)方面，4種方案的差異較大。圖7(a)為傳統(tǒng)方法，相對(duì)于圖7(b)、圖7(c)、圖7(d)來(lái)說(shuō)，建模模型為圓形區(qū)域，屏蔽了不同方位的鄰近點(diǎn)由于各向異性特征對(duì)待估點(diǎn)造成的影響，與實(shí)際砂體地層面不相符，精度較低。圖7(b)、圖7(c)、圖7(d)中，砂體地層面的表面模型具有明顯方向性。圖7(b)在水平方向變異明顯，插值結(jié)果反映東西方位地層走勢(shì)。圖7(c)在與水平方向呈60°夾角的方向上變異明顯，插值結(jié)果為東北-西南向地層走勢(shì)。圖7(d)在與水平方向呈θ度夾角的方向上變異明顯，插值結(jié)果說(shuō)明其為東南-西北向的地層走勢(shì)。實(shí)驗(yàn)數(shù)據(jù)實(shí)際地層走向與圖7(d)非常吻合，圖7(d)的精度高，符合預(yù)期。

通過(guò)對(duì)優(yōu)化方案最終表面的重構(gòu)結(jié)果進(jìn)行測(cè)試對(duì)比可發(fā)現(xiàn)，使用基于方向的距離平方反比插值的改進(jìn)方法，其結(jié)果有一定方向性變化，比較客觀地反映了砂體地層面的方向性特征。考慮到實(shí)際的應(yīng)用中，事物普遍存在著各向異性，如在砂體地層表面模型重建時(shí)，其水流方向性特別明顯，普通算法生成的砂體地層面就與實(shí)際不相符，精度較低。而本算法考慮各向異性，可以很好解決該問(wèn)題，改進(jìn)算法的插值效果明顯。

4 結(jié)束語(yǔ)

本文針對(duì)傳統(tǒng)距離平方反比插值法存在的問(wèn)題，提出了基于方向的各向異性插值算法，有效地克服了傳統(tǒng)方法的各向同性，解決了空間插值普遍存在的各向異性等問(wèn)題。該算法建立模型的重點(diǎn)在于確定合適的搜索區(qū)域參數(shù)，本文對(duì)模型中相關(guān)參數(shù)的求解過(guò)程進(jìn)行了詳細(xì)說(shuō)明，并對(duì)所提算法的插值效果與距離平方反比插值法進(jìn)行了測(cè)試。本方法結(jié)果與實(shí)際情況非常吻合，具有良好的精度與準(zhǔn)確度，更能真實(shí)反映空間數(shù)據(jù)的實(shí)際情況。