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

王兆國，程順有，劉 財

1.西北大學大陸動力學教育部重點實驗室，西安 710069

2.西北大學地質學系，西安 710069

3.吉林大學地球探測科學與技術學院，長春 130026

0 引言

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

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

1 插值方法的基本思想

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

反距離加權插值是利用平面上已知的一系列離散點，把離散點與插值點之間距離的倒數作為權重，進行插值[5-7]，公式如下：

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

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

其中：li（x）是與節點關聯的Voronoi邊的長度；hi（x）是節點到插值點的垂直距離。

最近鄰點插值法中隱含的假設是插值點的屬性值使用距離它最近的節點屬性值表示，區域變量平均屬性形式[5]為

其中，Si是由所有網格點連線的垂直平分線切割成的多邊形面積。

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

2 模型選取與計算

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

圖1 模型一Fig.1 Model 1

從圖2中可以看出：所有插值方法在相同地質體下，隨著采樣間距的增大，誤差均方根值總體趨勢是增大的，也就是采用小采樣間距時地質體誤差較?。凰蟹椒ㄔ谙嗤蓸娱g距下，隨著目標體的增大，誤差均方根值總體是增大的，產生此現象的原因是大目標體產生的場的影響范圍大，插值產生的誤差的范圍增大。

圖2 模型一不同插值方法采樣間距與地質體尺度關系的絕對誤差均方根值曲線圖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

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

圖3 模型二Fig.3 Model 2

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

對于模型二，筆者采用不同的采樣間距計算了不同網絡化插值方法的誤差大小，如表1所示。

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

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

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

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

3 結論

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

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

圖6 模型二不同插值方法節點處絕對誤差灰度圖Fig.6 Grey-scale map of absolute error of nodes of different interpolation methods of model 2

（References）：

[1]葛志廣，宋俊杰.高精度磁法數據網格化方法的選取[J].工程地球物理學報，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]張美根，烏達巴拉，王妙月.非規則網帶斷層地震數據的網格化[J].地球物理學進展，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]曾閩山，侯巖松.海量地震數據網格化算法分析與研究[J].石油天然氣學報：江漢石油學院學報，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網格化的研究[J].西南石油大學學報：自然科學版，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].工程地球物理學報，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]蔡玉華.達里亞工區速度體建立中幾種網格化方法的對比分析[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].物探與化探，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]郭良輝，孟小紅，郭志宏，等.地球物理不規則分布數據的空間網格化法[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]程紅杰，胡祥云，田米瑪，等.地震數據網格化方法研究[J].工程地球物理學報，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]郭良輝，孟小紅，郭志宏，等.反插值法實現地球物理數據快速網格化[J].地球物理學進展，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軟件的兩種數據插值方法研究[J].測繪通報，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]王兆清，馮偉.高度不規則網格多邊形單元的有理函數插值格式[J].固體力學學報，2005，26（2）：199-202.Wang Zhaoqing，Feng 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剖分的二維非規則重力建模及重力計算[J].地球物理學報，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.