孫夢楠,劉少華,劉京城
(長江大學 地球科學學院,湖北 武漢 430100)
近年來,隨著科學技術的快速發(fā)展,對空間數據的處理有了更高精度的要求,空間插值方法的優(yōu)化和應用越來越受到人們的重視。由于空間數據具有一定的空間依賴性與相關性,為了對待插值數據進行最優(yōu)內插計算,我們需要考慮數據空間分布的相關性和變異性,使得插值結果更趨于實際的要求。距離平方反比插值方法是眾多插值方法之一,是以已知點與待插值點之間的距離為基礎,在搜索鄰域內對待插點進行無偏最優(yōu)插值的一種方法[1]。
現(xiàn)實生活中,事物具有一定的相關性,同時也存在著各向異性。各向異性主要是指物體在不同方向上屬性值變化的速度具有一定的差異性,甚至屬性值有可能突然變化的情況[2,3]。各向異性是地質地層研究過程中常見的性質。例如,長江黃河中河道方向與垂直方向水體污染物的濃度分布規(guī)律,大氣環(huán)境污染物分布受到風向和重力的影響,礦產儲量分布受到地質構造環(huán)境的影響等等。在實際應用中,空間插值必須考慮物體在三維空間中的空間結構性和各向異性特征,才能構建與實體相統(tǒng)一的模型。
本文基于傳統(tǒng)的距離平方反比插值法,參考克里格方法[4],結合物體的各向異性特征,探討并建立了基于方向的各向異性插值優(yōu)化算法,通過對模型進行不同的參數設置,對相關結果進行比較,驗證了算法在空間插值方面的有效性。
距離平方反比插值算法的主要原理是,根據空間各已知點的坐標 (Xi,Yi)(i=1,2,…,n), 和各已知點屬性值Zi(i=1,2,…,n), 從而假設待插值點的坐標為 (x,y), 插值得到的結果為Z(x,y)。Z(x,y) 可理解為各已知點屬性值與待求結果的距離平方加權平均,如式(1)所示
(1)
該方法的算法簡單,計算量小,計算時消耗計算機的空間、花費的時間相對較小,時間與空間復雜度低,但存在如下問題[5]。
(1)只是一種純幾何的加權。若鄰域內多個點到待插點的距離相同,則它們的權重系數相同,這樣就忽視了不同方位的已知點由于各向異性對預測點造成的影響,即沒有考慮空間數據的各向異性特征。
(2)搜索半徑r值以及樣本點的選取無標準。插值結果較大程度上決定于鄰近點的選取和搜索范圍的確定。
各向異性是一種常見的現(xiàn)象,是物質的全部或部分化學、物理等性質隨方向的改變而變化的特性[6]??紤]到待插值數據具有各向異性特征,優(yōu)化方法將搜索區(qū)域由原來的圓形改為橢圓,如圖1所示,使得在搜索區(qū)域內的點更具有代表性,能夠對待插值點有所影響。對比圖1(a)、圖1(b)兩個圖形,發(fā)現(xiàn)區(qū)域內取得的建模點有所變化。
圖1 兩種方法的搜索區(qū)域對比
通過設置模型變量來求取橢圓公式,如圖2所示,其中a為橢圓長軸的值,b為橢圓短軸的值,θ為橢圓長軸延長線與x軸的夾角。模型參數可通過最小二乘法進行函數匹配來確定。由于最小二乘法橢圓擬合算法可能包含誤差較大的樣本點,該樣本點參與運算會使橢圓擬合的結果產生偏差[7,8]。針對這種情況,采用基于最小二乘法的橢圓擬合改進算法確定橢圓公式。隨機選取6個樣本點計算橢圓公式,統(tǒng)計與此公式匹配的所有樣本點個數,重復上述步驟一定次數。其中,匹配樣本點多的橢圓即為最優(yōu)橢圓,從而確定橢圓公式中參數的值[9,10]。
圖2 優(yōu)化模型變量三
一般方程式橢圓方程,如式(2)所示
Ax2+Bxy+Cy2+Dx+Ey+F=0
(2)
由于解的任意整數倍都是同一橢圓方程,且為了避免零值,設置約束條件如式(3)所示
A+C=1
(3)
通過求取目標函數(4)的最小值來確定各系數
(4)
根據極值原理,欲使目標函數值最小,必有式(5)
(5)
由此可得一個線性方程組,利用全主元高斯消去法并結合約束條件,求得方程系數A、B、C、D、E、F的值[11]。對橢圓方程進行化解,得到橢圓的標準方程,同時確定參數θ、a、b的值。
根據橢圓一般方程式和參數θ、a、b的值,可計算得到旋轉角度為θ的標準斜橢圓方程,在計算對待插值有影響的點時,得到橢圓式(6)
(6)
如果樣本點代入式(6),其值小于1,則為感興趣點,代入計算。
如圖3所示,通過點O(X0,Y0)、P(X1,Y1)、Q(X,Y), 求取直線OP公式,如式(7)所示
圖3 歸一化點位
(7)
為實現(xiàn)歸一化,需要求OP與OQ長度的比值,即k=(x1-x0)/(x-x0) 的值[12,13]。結合橢圓公式,可得到以下方程組(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)
推導公式如式(9)至式(15)所示,從而得到式(16),取得k值如下
(16)
根據k值,將原方法中距離d值歸一化,得d′=k*b(b為橢圓短軸)??紤]到物體的各向異性特征,通過異向比壓縮距離軸,使之成為各向同性的模型,即構成一個新的球狀模型,如圖4中虛線圓。將原方法各點距離平方權重改變,最后d′代入原方法d計算即可。
圖4 歸一化球狀模型
為驗證本文算法的可行性與精確度,本文選取某油田油井分層數據C1層的30個點層位數據作為測試數據。
將測試數據劃分為插值點和檢查點,選取10個點作為檢查點,對距離平方反比插值法和改進算法進行交叉驗證。實際值與兩種方法預測值對比如圖5、圖6所示,發(fā)現(xiàn)原方法中實際值與預測值繪制的點偏離直線較遠,插值結果與實際值的離散程度以及插值誤差的離散程度均大于本文方法。以誤差標準差、平均誤差、中誤差和方差作為插值結果的評價標準,評價標準值越小,表明插值結果越準確。數據結果見表1,后者的評價標準值均小于前者,其中距離平方反比插值法的實際值與預測值的方差較本文方法更大[14,15]。從而分析得出,本算法和傳統(tǒng)算法相比,精度相對較高,可行性較好,插值結果較準確。
圖5 實際值與原方法預測值對比
圖6 實際值與本方法預測值對比
統(tǒng)計參數距離平方反比插值法改進方法誤差標準差16.859 081 239.848 568 424平均誤差1597.4631594.18中誤差44.486 960 8542.719 062 8方差1979.089 6851824.918 326
根據最小二乘法的橢圓擬合改進算法確定橢圓參數,將橢圓旋轉不同的角度做空間插值計算,并與距離平方反比插值法進行比較。圖7為結果對比,其中圖7(a)為距離平方反比插值法的結果,圖7(d)為本算法的結果,圖7(b)、圖7(c)分別是將斜橢圓公式的θ值設為0°、60°的結果。
圖7 不同方向插值結果對比
鑒于傳統(tǒng)方法存在的問題,在進一步研究本算法時,首先考慮如何最大化利用數據方位的各向異性來提高建模精度。驗證方案提取的4種方式結果,從砂體地層面表面重構效果來看,傳統(tǒng)方法和優(yōu)化的方案在模型表面特征方面,存在一定的相似性,但在細節(jié)上差異較大。
從宏觀上看,4種方案預測的砂體地層面表面模型特征基本一致,砂體地層面的表面模型與已知層位數據關系明顯,局部區(qū)域的最值基本都是已知樣本點所在的位置,極值分布規(guī)律清晰,與地質、測井數據基本一致。在細節(jié)方面,4種方案的差異較大。圖7(a)為傳統(tǒng)方法,相對于圖7(b)、圖7(c)、圖7(d)來說,建模模型為圓形區(qū)域,屏蔽了不同方位的鄰近點由于各向異性特征對待估點造成的影響,與實際砂體地層面不相符,精度較低。圖7(b)、圖7(c)、圖7(d)中,砂體地層面的表面模型具有明顯方向性。圖7(b)在水平方向變異明顯,插值結果反映東西方位地層走勢。圖7(c)在與水平方向呈60°夾角的方向上變異明顯,插值結果為東北-西南向地層走勢。圖7(d)在與水平方向呈θ度夾角的方向上變異明顯,插值結果說明其為東南-西北向的地層走勢。實驗數據實際地層走向與圖7(d)非常吻合,圖7(d)的精度高,符合預期。
通過對優(yōu)化方案最終表面的重構結果進行測試對比可發(fā)現(xiàn),使用基于方向的距離平方反比插值的改進方法,其結果有一定方向性變化,比較客觀地反映了砂體地層面的方向性特征??紤]到實際的應用中,事物普遍存在著各向異性,如在砂體地層表面模型重建時,其水流方向性特別明顯,普通算法生成的砂體地層面就與實際不相符,精度較低。而本算法考慮各向異性,可以很好解決該問題,改進算法的插值效果明顯。
本文針對傳統(tǒng)距離平方反比插值法存在的問題,提出了基于方向的各向異性插值算法,有效地克服了傳統(tǒng)方法的各向同性,解決了空間插值普遍存在的各向異性等問題。該算法建立模型的重點在于確定合適的搜索區(qū)域參數,本文對模型中相關參數的求解過程進行了詳細說明,并對所提算法的插值效果與距離平方反比插值法進行了測試。本方法結果與實際情況非常吻合,具有良好的精度與準確度,更能真實反映空間數據的實際情況。