王美琪 李 建

(西南石油大學計算機科學學院 四川 成都 610500)

0 引 言

克里金插值法是一種在地統(tǒng)計學中常用的空間插值方法[1]。該插值方法最優(yōu)無偏，效果優(yōu)于距離加權(quán)法和趨勢面法，目前已應用于許多領(lǐng)域，幾乎是20世紀中后期以來應用最為廣泛的地理空間插值方法[2]，如地質(zhì)學、土壤學、生態(tài)學、氣象學、遙感甚至其他研究“時空變量”的領(lǐng)域[3-5]?？死锝鸩逯捣ㄖ饕譃閮纱蟛糠郑?1) 量化已知數(shù)據(jù)的空間相關(guān)性即變差函數(shù)，一般來說，距離較近的點屬性值越接近，距離較遠的點屬性值相差較多；(2) 確定鄰域范圍并搜索鄰域點，確定克里金方程組，求解克里金方程組得到權(quán)重系數(shù)，利用權(quán)重系數(shù)加權(quán)求和得到待插值點的屬性值。

在石油行業(yè)中，通過提取三維地震目的層的屬性數(shù)據(jù)[6]，計算得到屬性數(shù)據(jù)的分布平面圖，如振幅、頻率、屬性聚類等含油氣預測平面圖，為找油找氣提供目標方向。三維地震屬性平面圖的形成需要插值，由于三維地震屬性數(shù)據(jù)巨大，在實際應用中難以全區(qū)域計算，一般采用局部鄰域為每一個待插值點搜索出合適的鄰域點作為輸入。目前鄰域點搜索多采用指定半徑的圓或橢圓，將在圓或橢圓內(nèi)的點作為鄰域點輸入，該方法須遍歷，并計算所有已知點到待插值點的距離。文獻[7]漂移模式和搜索范圍參數(shù)的影響較小,可忽略不計；文獻[8]改進了鄰域克里金插值法；文獻[9]采用Delaunay三角劃分搜索和固定距離的方式搜索點；文獻[10]引入KD-Tree的二維索引結(jié)構(gòu)，提出改進J鄰近點搜索方法，效率有所提高，但該類搜索仍需計算已知點到待插值點的距離，這無形中增加了計算量。本文提出了一種無須距離計算、無須遍歷已知點而直接利用待插值點位置的增減實現(xiàn)鄰域點的選取算法，快速實現(xiàn)克里金插值。經(jīng)過實驗證明，在同精度下，該算法比指定距離半徑的圓或橢圓搜索算法快，特別對于海量數(shù)據(jù)，該算法更顯優(yōu)勢。

1 VAOS基本原理

數(shù)據(jù)插值是通過已知鄰近數(shù)據(jù)點獲得未知數(shù)據(jù)點，使數(shù)據(jù)規(guī)則化、網(wǎng)格化的過程。如何快速確定待內(nèi)插網(wǎng)格點周圍的已知鄰近數(shù)據(jù)點是提高克里金插值效率的關(guān)鍵。

1.1 坐標到網(wǎng)格的轉(zhuǎn)換

假設空間數(shù)據(jù)點的坐標為P(x,y,z)，可轉(zhuǎn)換為網(wǎng)格位置P′(ndx,ndy,ndz)來表達，n為網(wǎng)格數(shù)，dx、dy、dz為網(wǎng)格間距。將已知數(shù)據(jù)坐標統(tǒng)一到待插值網(wǎng)格，通過待插值網(wǎng)格位置的增減毫不費時找到已知點P(x,y,z)，如圖1所示。而已知點原始位置未變，不影響計算，因坐標轉(zhuǎn)換取整有可能幾點成為同一點，說明這幾點相近。

圖1 坐標網(wǎng)格統(tǒng)一

假設：(1) 已知點平面位置表示為P(x,y)，x、y最小值分別為xmin、ymin；(2) 待插值點網(wǎng)格方向x、y起始位置表示為nx0、ny0；(3) 待插值點方向網(wǎng)格步長表示為dx、dy。

則已知點P(x,y)在待插值點的網(wǎng)格位置可表示為P′(nx,ny)，其中：

已知點統(tǒng)一到在待插值網(wǎng)格位置后，只需通過插值點網(wǎng)格位置的增減查找已知點，而已知點的原始位置并未變化，參加克里金插值時不受影響。

坐標到網(wǎng)絡轉(zhuǎn)換的偽代碼如下：

輸入：已知點P(x,y)個數(shù)為N；待插值點網(wǎng)格方向起始位置表示為nx0、ny0,待插值點方向網(wǎng)格步長表示為dx、dy。

(1) 已知點N中找到最小的x、y值xmin、ymin。

for (i=0;i

{

x>xmin;xmin=x;

//最小的x

y>ymin;ymin=y;

//最小的y

}

(2) 已知點轉(zhuǎn)換為待插值得網(wǎng)格位置nx、ny。

for (i=0;i

{

nx=(x-xmin)/dx+nx0

ny=(x-xmin)/dx+nx0

flag(nx,ny)=1

//標示該位置有已知點值

}

1.2 待插值點位置的增減搜索鄰域點算法

假設搜索的待插值點為w(i,j)，如圖2所示。

圖2 鄰近點搜索算法

(1) 待插值點w(i,j)，x、y方向網(wǎng)格位置為i、j。

(2) 待插值點w(i,j)，x、y方向增減網(wǎng)格總數(shù)nbx、nby。

(3) 待插值點w(i,j)，x、y方向網(wǎng)格間距ndx=1,ndy=1。

(4) 待插值點網(wǎng)格x、y方向起始位置表示為：nx0=0,ny0=0。

(5) 待插值點x、y方向終止位置表示為：nxn,nyn。

則待插值點w(i,j)已知鄰近點的搜索范圍實際變?yōu)橛蒳-nbx、i+nbx、j-nby、j+nby網(wǎng)格構(gòu)成的邊長為2nbx、2nby的四邊形。判斷四邊形內(nèi)的P′(nx,ny)是否有數(shù)據(jù)。

這樣通過待插值點下標位置變量(i,j)的增減循環(huán)快速實現(xiàn)對所有待插值點鄰近已知點的搜索。該過程偽代碼如下：

//為所有待插值點搜索鄰近已知點

for (i=0;i

//x方向位置變化

{

for (j=0;j

//y方向位置變化

{

//為待插值點w(i,j)搜索鄰近已知點

//設鄰近已知點：N=0

for (m=i-nbx;m

//i方向

//判斷從i-nbx到i+nbx范圍

{

m<0

//進行下次循環(huán)

m>nxn

//進行下次循環(huán)

for (k=j-nby;k

//j方向

//判斷從j-nby,j+nby范圍

{

k<0

//進行下次循環(huán)

k>nxn

//進行下次循環(huán)

如P′(m,n)有數(shù)據(jù)n=n+1;

}

}

//克里金插值

}

}

1.3 VAOS算法流程

VAOS算法流程如圖3所示。

圖3 VAOS算法流程

2 驗證與應用

為驗證本計算方法，用實際三維地震屬性數(shù)據(jù)進行插值，并與傳統(tǒng)指定距離半徑搜索鄰近域克里金插值算法、Surfer軟件插值算法[11-12]的結(jié)果進行對比。

該工區(qū)主測線290-420(131條)，聯(lián)絡線320-400(81點)，樣本點數(shù)總點數(shù)為N=131×81=10 611，待插值網(wǎng)格點數(shù)為M=200×200=40 000。

2.1 對比分析驗證

(1) 效果圖對比。三種算法插值的效果圖對比如圖4所示。

(a) Surfer插值結(jié)果 (b) 距離半徑 (c) VAOS方法圖4 算法插值效果

(2) 時間、復雜度對比。距離半徑搜索算法：進行一次最鄰近點w(i,j)搜索的時間復雜度為O(N)，并計算待插值點到N點距離，并判斷是否在距離半徑內(nèi)，本次半徑距離為3，獲得的鄰近已知最小5個點，最大12個點，插值時間為135.6 s。VAOS搜索算法：進行一次最鄰近點搜索的時間復雜度為O(nbx×nby)，不需要計算待插值點到N點距離，只需要判斷是否有值。本次給定nbx=3，nby=3，獲得的鄰近已知最小5個點，最大12個點，插值時間為2.79 s。

表1 算法插值時間對比表

2.2 三維地震油氣預測應用

預測目的層如圖5(a)所示；提取振幅、頻率等多種屬性進行聚類，本文算法克里金插值繪制油氣預測結(jié)果圖如圖5(b)所示。

(a) 預測目的層

3 結(jié) 語

通過對實際三維地震屬性數(shù)據(jù)插值效果圖、時間和復雜度進行對比，可以看出：插值的效果圖基本一致，距離半徑搜索算法復雜度遠高于VAOS搜索算法和Surfer，指定距離半徑搜索算法需遍歷所有已知，并涉及距離計算，十分費時。VAOS算法無須距離計算、無須遍歷已知點，直接利用待插值點位置的增減實現(xiàn)鄰域點的選取，速度極快。

VAOS算法將遍歷數(shù)據(jù)搜索算法變?yōu)閳D像位置的增減可謂是一個創(chuàng)新，該算法搜索速度快，實用性強。