VAOS克里金插值在三維地震屬性平面圖的應用

王美琪 李 建

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

0 引 言

克里金插值法是一種在地統計學中常用的空間插值方法[1]。該插值方法最優無偏，效果優于距離加權法和趨勢面法，目前已應用于許多領域，幾乎是20世紀中后期以來應用最為廣泛的地理空間插值方法[2]，如地質學、土壤學、生態學、氣象學、遙感甚至其他研究“時空變量”的領域[3-5]。克里金插值法主要分為兩大部分：(1) 量化已知數據的空間相關性即變差函數，一般來說，距離較近的點屬性值越接近，距離較遠的點屬性值相差較多；(2) 確定鄰域范圍并搜索鄰域點，確定克里金方程組，求解克里金方程組得到權重系數，利用權重系數加權求和得到待插值點的屬性值。

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

1 VAOS基本原理

數據插值是通過已知鄰近數據點獲得未知數據點，使數據規則化、網格化的過程。如何快速確定待內插網格點周圍的已知鄰近數據點是提高克里金插值效率的關鍵。

1.1 坐標到網格的轉換

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

圖1 坐標網格統一

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

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

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

坐標到網絡轉換的偽代碼如下：

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

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

for (i=0;i

{

x>xmin;xmin=x;

//最小的x

y>ymin;ymin=y;

//最小的y

}

(2) 已知點轉換為待插值得網格位置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方向網格位置為i、j。

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

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

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

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

則待插值點w(i,j)已知鄰近點的搜索范圍實際變為由i-nbx、i+nbx、j-nby、j+nby網格構成的邊長為2nbx、2nby的四邊形。判斷四邊形內的P′(nx,ny)是否有數據。

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

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

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

//進行下次循環

m>nxn

//進行下次循環

for (k=j-nby;k

//j方向

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

{

k<0

//進行下次循環

k>nxn

//進行下次循環

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

}

}

//克里金插值

}

}

1.3 VAOS算法流程

VAOS算法流程如圖3所示。

圖3 VAOS算法流程

2 驗證與應用

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

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

2.1 對比分析驗證

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

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

(2) 時間、復雜度對比。距離半徑搜索算法：進行一次最鄰近點w(i,j)搜索的時間復雜度為O(N)，并計算待插值點到N點距離，并判斷是否在距離半徑內，本次半徑距離為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)所示；提取振幅、頻率等多種屬性進行聚類，本文算法克里金插值繪制油氣預測結果圖如圖5(b)所示。

(a) 預測目的層

3 結 語

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

VAOS算法將遍歷數據搜索算法變為圖像位置的增減可謂是一個創新，該算法搜索速度快，實用性強。