地球物理勘探中幾種二維插值方法的誤差分析

王兆國，程順有，劉 財(cái)

1.西北大學(xué)大陸動(dòng)力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，西安 710069

2.西北大學(xué)地質(zhì)學(xué)系，西安 710069

3.吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

在地球物理工作以及地質(zhì)勘探中，為確?？蓪?shí)現(xiàn)性和可操作性，往往對(duì)研究區(qū)進(jìn)行離散觀測(cè)，然后對(duì)離散的數(shù)據(jù)經(jīng)過一系列處理形成時(shí)間或空間上的離散數(shù)據(jù)體，在離散數(shù)據(jù)體的圖像呈現(xiàn)上常采用二維插值方法進(jìn)行處理。二維插值方法已經(jīng)應(yīng)用于許多領(lǐng)域，如：葛志廣等[1]用不同的網(wǎng)格化方法對(duì)高精度磁法數(shù)據(jù)進(jìn)行了分析；張美根等[2]和曾閩山等[3]分別就非規(guī)則地震數(shù)據(jù)和海量散點(diǎn)地震數(shù)據(jù)進(jìn)行了網(wǎng)格化研究；陳曉軍等[4]對(duì)油藏?cái)?shù)據(jù)進(jìn)行了Voronoi網(wǎng)格化研究等。對(duì)于各種二維插值方法，如反距離加權(quán)插值法、克里金插值法、最小曲率法、改進(jìn)的謝別德法、自然鄰點(diǎn)插值法、最近鄰點(diǎn)插值法、多元回歸法、徑向基函數(shù)法、線性插值三角網(wǎng)法、移動(dòng)平均法、局部多項(xiàng)式插值和反插值法，在前人的研究中就其基本原理和效果圖的對(duì)比都有過論述[5-11]。研究圖件要反映整個(gè)研究區(qū)異常特征，但它是離散采樣的結(jié)果，因此就形成了圖件與真實(shí)情況差異問題，即誤差問題。成像誤差的影響因素包含3個(gè)方面的內(nèi)容：一是地質(zhì)體本身大小的影響；二是采樣密度的影響；三是網(wǎng)格化插值方法的影響。此前的研究主要考慮選取不同插值方法的影響，但沒有研究各種插值方法具體的誤差問題，而對(duì)于采樣密度影響幾乎沒有論述。

筆者主要針對(duì)采樣密度和二維插值方法對(duì)成像的誤差大小進(jìn)行研究。為了確定其影響情況，模擬了2個(gè)重力異常場(chǎng)，然后通過不同的采樣間距和網(wǎng)格化方法對(duì)模擬數(shù)據(jù)進(jìn)行不同插值方法的插值，對(duì)比插值數(shù)值與真實(shí)值之間的誤差，從而確定插值方法本身的精度情況。

1 插值方法的基本思想

網(wǎng)格化方法實(shí)質(zhì)上是一種利用已知點(diǎn)值進(jìn)行二維空間插值的方法，各種不同方法的主要區(qū)別在于利用的采樣點(diǎn)范圍的不同和已知點(diǎn)權(quán)重的不同。采樣點(diǎn)范圍的不同表現(xiàn)為全局插值和局部插值，而已知點(diǎn)權(quán)重的不同在于權(quán)重函數(shù)（基函數(shù)）的不同。由于多元回歸法是一種趨勢(shì)面作圖法，局部多項(xiàng)式插值法也可以看作是局部趨勢(shì)面法，即用一個(gè)曲面去擬合已知數(shù)據(jù)，可以采用不同冪次，主要用來區(qū)分區(qū)域場(chǎng)與局域場(chǎng)，因此方法本身就對(duì)精細(xì)的構(gòu)造起到一種去除作用，筆者不再討論。移動(dòng)平均法，通過“窗”進(jìn)行搜索，擬合程度對(duì)“窗”的大小非常依賴，“窗”選擇不合適往往造成圖的畸變，當(dāng)采用不同采樣間距時(shí)，主要是“窗”的選取問題，因此筆者也不再討論。筆者主要討論反距離加權(quán)插值法、克里金插值法、最小曲率法、改進(jìn)的謝別德法、自然鄰點(diǎn)法、最近鄰點(diǎn)法、徑向基函數(shù)法和線性插值三角網(wǎng)法：反距離加權(quán)法容易受到數(shù)據(jù)點(diǎn)集群的影響，計(jì)算結(jié)果中常出現(xiàn)孤立點(diǎn)數(shù)據(jù)明顯高于周圍數(shù)據(jù)點(diǎn)的現(xiàn)象；克里金法利用變差函數(shù)描述函數(shù)值的區(qū)域變化，不僅考慮了已知點(diǎn)對(duì)待插點(diǎn)的影響，而且考慮了已知點(diǎn)之間的相互影響，網(wǎng)格化精度效果的好壞，依賴于理論變差函數(shù)擬合的好壞；最小曲率法主要考慮曲面的光滑性，容易超出最大值和最小值的范疇；改進(jìn)的謝別德法是對(duì)反距離加權(quán)法的改進(jìn)，利用節(jié)點(diǎn)函數(shù)的二次曲面擬合值代替離散點(diǎn)值，提高了內(nèi)插值精度和曲面光滑度，權(quán)函數(shù)只在局部起作用，并使用最小二乘法克服了反距離加權(quán)的“牛眼”缺點(diǎn)；自然鄰點(diǎn)法是基于Voronoi結(jié)構(gòu)的一類插值方法；最近鄰點(diǎn)法假設(shè)任意網(wǎng)格點(diǎn)的屬性值都用距離它最近位置點(diǎn)的屬性值，更適合于變化不大且測(cè)量點(diǎn)密度大的情況；徑向基函數(shù)法是多個(gè)數(shù)據(jù)插值方法組合的一種多形式網(wǎng)格化方法，具有很強(qiáng)的擬合數(shù)據(jù)點(diǎn)、產(chǎn)生光滑曲面的能力；線性插值的三角網(wǎng)法計(jì)算速度快，網(wǎng)格化結(jié)果不光滑，精度和效果較差，計(jì)算可能會(huì)出現(xiàn)不穩(wěn)定。

反距離加權(quán)插值是利用平面上已知的一系列離散點(diǎn)，把離散點(diǎn)與插值點(diǎn)之間距離的倒數(shù)作為權(quán)重，進(jìn)行插值[5-7]，公式如下：

其中：di（x，y）＝，表示離散點(diǎn)（xi，yi）與插值點(diǎn)（xp，yp）間的距離；Zi是各離散點(diǎn)的屬性值；Zp為插值點(diǎn)p的屬性值。

改進(jìn)的謝別德插值法是對(duì)反距離加權(quán)插值的改進(jìn)，主要改進(jìn)有2個(gè)方面：一是把權(quán)重函數(shù)由全局插值改為局部插值；二是引入節(jié)點(diǎn)函數(shù)代替離散點(diǎn)的屬性值[5]，但權(quán)重依然采用距離的倒數(shù)?？死锝鸩逯捣ㄊ且环N對(duì)空間分布數(shù)據(jù)求最優(yōu)、線性、無偏內(nèi)插估計(jì)量的方法。它的優(yōu)點(diǎn)在于不僅考慮了各已知數(shù)據(jù)點(diǎn)的空間相關(guān)性，而且在給出待估計(jì)點(diǎn)數(shù)值的同時(shí)，還能給出表示估計(jì)精度的方差?？死锝鸱匠探M是在無偏性條件和估計(jì)方差最小條件下給出的[5-8]。最小曲率法構(gòu)造具有最小曲率的曲面，使其穿過空間場(chǎng)的每一點(diǎn)，在盡量嚴(yán)格地尊重?cái)?shù)據(jù)的同時(shí)，生成盡可能圓滑的曲面，最小曲率法主要考慮曲面的光滑性，因此插值容易超出最大和最小值的范疇[5]。自然鄰點(diǎn)插值法是基于Voronoi結(jié)構(gòu)的一類插值方法，其權(quán)重函數(shù)形式[5，12]為

其中：li（x）是與節(jié)點(diǎn)關(guān)聯(lián)的Voronoi邊的長度；hi（x）是節(jié)點(diǎn)到插值點(diǎn)的垂直距離。

最近鄰點(diǎn)插值法中隱含的假設(shè)是插值點(diǎn)的屬性值使用距離它最近的節(jié)點(diǎn)屬性值表示，區(qū)域變量平均屬性形式[5]為

其中，Si是由所有網(wǎng)格點(diǎn)連線的垂直平分線切割成的多邊形面積。

徑向基函數(shù)法是多個(gè)數(shù)據(jù)插值方法的組合，其基函數(shù)由單變量函數(shù)構(gòu)成，復(fù)二次基函數(shù)法被認(rèn)為是其中最佳的。線性插值三角網(wǎng)法是基于delaunay三角網(wǎng)剖分的插值方法[13-14]，要保證2個(gè)前提：所有三角形的邊不能相交；保證三角形最小內(nèi)角和為最大。這2個(gè)條件保證盡可能避免生成小內(nèi)角的長薄單元，使三角形能夠接近等角或等邊[15]。

2 模型選取與計(jì)算

為了研究不同采樣間距和插值方法對(duì)成圖誤差的影響，首先確定一個(gè)簡單的模型（模型一）：地質(zhì)體為一長方體，上表面位于10km深度，下底面位于15km深度，水平2個(gè)方向長度相同，分別取0.5，2.0，4.0，6.0，8.0，10.0，15.0，20.0km 的長度，地質(zhì)體在地表的投影中心為坐標(biāo)原點(diǎn)。用如圖1所示的模型，計(jì)算了不同地質(zhì)體尺度（C）分別在0.5，1.0，2.5，4.0，5.0，8.0，10.0，20.0km 的 采樣 間 距和不同插值方法下的絕對(duì)誤差的均方根值（圖2）。

圖1 模型一Fig.1 Model 1

從圖2中可以看出：所有插值方法在相同地質(zhì)體下，隨著采樣間距的增大，誤差均方根值總體趨勢(shì)是增大的，也就是采用小采樣間距時(shí)地質(zhì)體誤差較??；所有方法在相同采樣間距下，隨著目標(biāo)體的增大，誤差均方根值總體是增大的，產(chǎn)生此現(xiàn)象的原因是大目標(biāo)體產(chǎn)生的場(chǎng)的影響范圍大，插值產(chǎn)生的誤差的范圍增大。

圖2 模型一不同插值方法采樣間距與地質(zhì)體尺度關(guān)系的絕對(duì)誤差均方根值曲線圖Fig.2 Curve graph of absolute error mean square value of relationship of sampling spacing and geological body scale of different interpolation methods in model 1

為了進(jìn)一步比較不同網(wǎng)格化方法誤差的差異，選用一個(gè)更一般的重力模型（模型二，圖3）：底部是2層正剩余密度層，被一正剩余密度和一負(fù)剩余密度的薄板狀異常體所隔斷；薄板狀異常體的傾角為45°，兩薄板異常體的寬度均為10km；左側(cè)的2層地層厚度分別為5km和10km，上頂面位于10km深度，下底面位于25km深度；右側(cè)2層地層厚度分別為5km和10km，上頂面位于15km深度，下底面位于30km深度。這些異常是作為整個(gè)區(qū)域的區(qū)域異常體，其在深度-x面上的投影如圖3b所示。2層正剩余密度層上方在薄板狀體兩側(cè)各有一長方體負(fù)剩余密度體，為了更好地研究局部異常的影響，在右側(cè)的長方體負(fù)剩余密度體左側(cè)，建立了3個(gè)正方體正剩余密度體。這些異常作為局部異常體，5個(gè)局部異常體的厚度均為5km，其在x-y面上的投影如圖3c所示。模型二原始數(shù)據(jù)特征如圖4所示，上部為立體圖，下部為投影圖。

圖3 模型二Fig.3 Model 2

圖4 模型二重力異常特征Fig.4 Characteristics of gravity anomaly of model 2

對(duì)于模型二，筆者采用不同的采樣間距計(jì)算了不同網(wǎng)絡(luò)化插值方法的誤差大小，如表1所示。

對(duì)模型二，固定地質(zhì)體尺度，采用不同的采樣間距，計(jì)算不同插值方法的均方根值，來反映不同方法插值的好壞程度。由圖5可以看出：徑向基函數(shù)法在插值成圖時(shí)誤差最小，其次為改進(jìn)的謝別德法和克里金插值法，自然鄰點(diǎn)法和線性插值三角網(wǎng)法在圖上的折線幾乎重合，誤差大小幾乎相同，誤差較大的為反距離加權(quán)插值法、最小曲率法和最近鄰點(diǎn)法；當(dāng)采樣間距小于4km（最小異常地質(zhì)體尺度）時(shí)，誤差較大的后3種方法絕對(duì)誤差均方根值由小到大的順序是：反距離加權(quán)插值法、最近鄰點(diǎn)法、最小曲率法；當(dāng)采樣間距大于4km時(shí)，這3種方法絕對(duì)誤差均方根值由小到大的順序是：最小曲率法、最近鄰點(diǎn)法、反距離加權(quán)插值法。

表1 模型二不同插值方法的絕對(duì)誤差均方根值Table1 Absolute error root mean square values of different interpolation methods of model 2

圖5 模型二不同插插方法的絕對(duì)誤差均方根值曲線圖Fig.5 Curve graph of absolute error mean square value of different interpolation methods of model 2

絕對(duì)誤差均方根值反映了方法總體插值誤差的大小，為了進(jìn)一步觀測(cè)插值時(shí)重力異常任何一個(gè)節(jié)點(diǎn)處的絕對(duì)誤差，對(duì)模型二用采樣間距為0.5km的插值節(jié)點(diǎn)的絕對(duì)誤差來反映內(nèi)部節(jié)點(diǎn)處的誤差（圖6）。從圖6可以看出：徑向基函數(shù)法（圖6g）、改進(jìn)的謝別德法（圖6d）和克里金插值法（圖6b），所有節(jié)點(diǎn)上誤差幾乎相同，繪圖時(shí)對(duì)于地質(zhì)體邊界數(shù)據(jù)的影響最小；而自然鄰點(diǎn)法（圖6e）、線性插值三角網(wǎng)法（圖6h）、反距離加權(quán)插值法（圖6a）、最小曲率法（圖6c），在局部地質(zhì)體邊上范圍內(nèi)絕對(duì)誤差相對(duì)較大，對(duì)于繪圖時(shí)地質(zhì)體邊界的影響大；而最近鄰點(diǎn)法（圖6f），出現(xiàn)局部異常體形態(tài)的明暗相間的灰度圖，最暗（負(fù)）和最明（正）都意味著誤差大，因此最近鄰點(diǎn)法對(duì)繪圖時(shí)誤差影響大。

3 結(jié)論

1）對(duì)于同一插值方法而言，存在小間距絕對(duì)誤差均方根值小于大間距絕對(duì)誤差均方根值的關(guān)系。對(duì)不同的插值方法而言，當(dāng)采樣間距小于4.0km時(shí)，絕對(duì)誤差均方根值由小到大的順序是：徑向基函數(shù)法、改進(jìn)的謝別德法、克里金插值法、自然鄰點(diǎn)法、反距離加權(quán)插值法、最近鄰點(diǎn)法、最小曲率法，并且線性插值三角網(wǎng)法與自然鄰點(diǎn)法具有幾乎相同的數(shù)值；采樣間距大于4.0km時(shí)，絕對(duì)誤差均方根值由小到大的順序是：徑向基函數(shù)法、改進(jìn)的謝別德法、克里金插值法、自然鄰點(diǎn)法、最小曲率法、最近鄰點(diǎn)法、反距離加權(quán)插值法，并且線性插值三角網(wǎng)法和自然鄰點(diǎn)法具有幾乎相同的數(shù)值。從絕對(duì)誤差均方值看，徑向基函數(shù)方法、改進(jìn)的謝別德方法和克里金方法較小，其中徑向基函數(shù)值絕對(duì)誤差均方根值最小。

2）從節(jié)點(diǎn)處絕對(duì)誤差值來看，徑向基函數(shù)方法、克里金方法、改進(jìn)的謝別德方法相對(duì)其他插值方法具有更小的誤差，不存在局部誤差較小或較大的情況，是相對(duì)較好的插值方法，并且徑向基函數(shù)方法是最好的。

圖6 模型二不同插值方法節(jié)點(diǎn)處絕對(duì)誤差灰度圖Fig.6 Grey-scale map of absolute error of nodes of different interpolation methods of model 2

（References）：

[1]葛志廣，宋俊杰.高精度磁法數(shù)據(jù)網(wǎng)格化方法的選取[J].工程地球物理學(xué)報(bào)，2010，7（2）：169-172.Ge Zhiguang，Song Junjie.The Griding Methods of Selection About the Data of High-Precision Magnetic[J].Chinese Journal of Engineering Geophysic，2010，7（2）：169-172.

[2]張美根，烏達(dá)巴拉，王妙月.非規(guī)則網(wǎng)帶斷層地震數(shù)據(jù)的網(wǎng)格化[J].地球物理學(xué)進(jìn)展，1999，14（1）：69-77.Zhang Meigen，Wudabala，Wang Miaoyue.Gridding of Inregular Seismic Data with Faults[J].Progress in Geophysics，1999，14（1）：69-77.

[3]曾閩山，侯巖松.海量地震數(shù)據(jù)網(wǎng)格化算法分析與研究[J].石油天然氣學(xué)報(bào)：江漢石油學(xué)院學(xué)報(bào)，2006，28（2）：72-75.Zeng Minshan，Hou Yansong.The Analysis and Research on Gridding Algorithm for Massive Seismic Data[J].Jounral of Oil and Gas Technology：J JPI，2006，28（2）：72-75.

[4]陳曉軍，陳偉，段永剛，等.油藏Voronoi網(wǎng)格化的研究[J].西南石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，32（1）：121-124.Chen Xiaojun，Chen Wei，Duan Yonggang，et al.The Research on Reservoir Voronoi Grid[J].Journal of Southwest Petroleum University： Science ＆Technology Edition，2010，32（1）：121-124.

[5]陳歡歡，李星，丁文秀.Surfer 8.0等值線繪制中的十二種插值方法[J].工程地球物理學(xué)報(bào)，2007，4（1）：52-57.Chen Huanhuan，Li Xing，Ding Wenxiu.Twelve Kinds of Gridding Methods of Surfer 8.0in Isoline Drawing [J].Chinese Journal of Engineering Geophysics，2007，4（1）：52-57.

[6]蔡玉華.達(dá)里亞工區(qū)速度體建立中幾種網(wǎng)格化方法的對(duì)比分析[J].石油物探，1995，34（4）：100-108.Cai Yuhua.Gridding Methods Used for Constructing Velocity Volume in Daliya：Comparative Analysis[J].Geophysical Prospecting for Petrol，1995，34（4）：100-108.

[7]劉兆平，楊進(jìn)，武煒.地球物理數(shù)據(jù)網(wǎng)格化方法的選取[J].物探與化探，2010，34（1）：93-97.Liu Zhaoping，Yang Jin，Wu Wei.The Choice of Gridding Methods for Geophysical Data [J].Geophysical and Geochemical Exploration，2010，34（1）：93-97.

[8]郭良輝，孟小紅，郭志宏，等.地球物理不規(guī)則分布數(shù)據(jù)的空間網(wǎng)格化法[J].物探與化探，2005，29（5）：438-442.Guo Lianghui，Meng Xiaohong，Guo Zhihong，et al.Gridding Methods of Geophysical Irregular Data in Space Domain [J].Geophysical and Geochemical Exploration，2005，29（5）：438-442.

[9]程紅杰，胡祥云，田米瑪，等.地震數(shù)據(jù)網(wǎng)格化方法研究[J].工程地球物理學(xué)報(bào)，2006，3（1）：28-32.Cheng Hongjie，Hu Xiangyun，Tian Mima，et al.The Study of Seismic Data Gridding Methods[J].Chinese Journal of Engineering Geophysics，2006，3（1）：28-32.

[10]郭良輝，孟小紅，郭志宏，等.反插值法實(shí)現(xiàn)地球物理數(shù)據(jù)快速網(wǎng)格化[J].地球物理學(xué)進(jìn)展，2005，20（3）：671-676.Guo Lianghui，Meng Xiaohong，Guo Zhihong，et al.Fast Gridding of Geophysical Data with Inverse Interpolation[J].Progress in Geophysics，2005，20（3）：671-676.

[11]馬英蓮，彭樹宏，錢靜.基于Surfer軟件的兩種數(shù)據(jù)插值方法研究[J].測(cè)繪通報(bào)，2010（8）：54-57.Ma Yinglian，Peng Shuhong，Qian Jing.A Study of Two Data Interpolation Methods Based on Surfer Software[J].Bulletin of Surveying and Mapping，2010（8）：54-57.

[12]王兆清，馮偉.高度不規(guī)則網(wǎng)格多邊形單元的有理函數(shù)插值格式[J].固體力學(xué)學(xué)報(bào)，2005，26（2）：199-202.Wang Zhaoqing，F(xiàn)eng Wei.Rational Function Interpolation Scheme of Polygonal Elements Based on Highly Irregular Grids[J].Acta Mechanica Solida Sinica，2005，26（2）：199-202.

[13]Green P J，Sibson R.Computing Direchlet Tesselations in the Plane[J].Computer Journal，1978，21（2）：168-173.

[14]Watson D F.Computing the N-Dimensional Delaunay Tessellation with Application to Voronoi Polytopes[J].Computer Journal，1981，24（2）：167-172.

[15]張嶺，郝天珧.基于Delaunay剖分的二維非規(guī)則重力建模及重力計(jì)算[J].地球物理學(xué)報(bào)，2006，49（3）：877-884.Zhang Ling，Hao Tianyao.2-D Irregular Gravity Modeling and Computation of Gravity Based on Delaunay Triangulation [J].Chinese Journal of Geophysics，2006，49（3）：877-884.