插值算法在物探數據處理中的應用

冉根領，王 磊

（1貴州地礦基礎工程有限公司，貴州 貴陽 550000；2重慶市二零八工程檢測有限公司，重慶 400700）

在工程勘察或地質調查中，由于地形地物、記錄及經濟成本等主客觀原因，在工區內只能獲得有限的鉆孔和物探資料，而這些資料往往是一系列離散的、空間上分布不均勻的點數據，對許多地質現象的解釋以及構建地下三維數據模型往往都是基于這些點數據做出的[1]。對已有的離散數據點選用合適的插值算法加密有效數據來最大限度地利用這些數據所蘊涵的信息，這是地質、物探實現二、三維數據可視化的關鍵步驟也是許多地球物理學家和計算機圖形學家一直研究的課題方向[2]。

就物探數據而言，其由于野外客觀因素的影響采集獲得的數據往往不完整。數據的缺失不僅僅意味著地下信息的丟失，而且還會在數據反演處理過程中帶入不必要的噪聲，從而影響物探的分辨率和信噪比[3]。另外，由于現在的許多反演算法都是以均勻網格為基礎開發的，非均勻數據可能造成算法的不能運行或運行結果不可預測。因此，需要對不規則的離散數據進行插值來實現網格重建。

到目前為止，國內外插值研究的相關著作多達4000多項，許多文獻對離散數據的插值做過深入研究，并且提出許多行之有效的算法。作者沿著前人的腳步，對物探數據處理中的插值和物探解釋中三維成像進行了一些研究，編寫程序實現插值算法并應用于高密度數據處理和物探數據的三維地質解釋中。

1 方法原理

1.1 Delaunay三角剖分、四面體剖分

Voronoi圖的簡單類比理解就是：將平面看做一水平面、數據點的集合看作n個水面上的振動源點，源點振動水波開始以相同的速度向四面傳播，不同點振動所引起的水波相遇時著振動停止，最后水波紋相遇的地方所形成的線構成的網格就是Voronoi圖。如圖1。

圖1 Voronoi圖

圖2 Voronoi圖與Delaunay三角剖分

二維平面中，離散點集中無四點共圓，對應的voronoi網格如圖2右所示，每個獨立的單元稱為泰森多邊形，若兩個單元之間共用一條邊，則將兩個單元的中心點連線形成Delaunay三角剖分的網格。

三維空間中，將離散點集看作空間中一些向外發射電磁波的信號源，當不同發射源的電磁信號相遇時就形成了Voronoi圖；離散點集中無五點共球，對應的voronoi網格，若兩個空間單元之間共用一個面，則將兩個空間單元的中心點連線形成Delaunay四面體剖分的網格，如圖3。

Delaunay三角剖分所具備的優異特性[4]：

（1）最接近：以最近的三點形成三角形，且各線段(三角形的邊)皆不相交。

（2）唯一性：不論從區域何處開始構建，最終都將得到一致的結果。

（3）最優性：任意兩個相鄰三角形形成的凸四邊形的對角線如果可以互換的話，那么兩個三角形六個內角中最小的角度不會變大。

（4）最規則：如果將三角網中的每個三角形的最小角進行升序排列，則Delaunay三角網的排列得到的數值最大。

（5）區域性：新增、刪除、移動某一個頂點時只會影響臨近的三角形。

（6）具有凸多邊形的外殼：三角網最外層的邊界形成一個凸多邊形的外殼。

同理Delaunay四面體剖分也具有相應的性質。利用Delaunay三角剖分、四面體剖分來進行離散數據的網格化可獲得穩定、唯一的網格。

圖3 Delaunay四面體剖分

本文Delaunay剖分的算法采用逐點插入法，它由Lawson于1977年提出，隨后由Bowyer，Wattson等人對其進行了改進[5]：

（1）構建初始大三角形(0,3＊max)、(-3＊max,3＊max)、(3＊max,0)包圍整個數據點集P，其中max為數據點集中x、y絕對值的最大值；

（2）隨機的對點集P中的數據點進行排序；

（3）按照隨機排好的順序將點集P中的點依次插入到初始大三角形中，在整個過程中都要隨時檢查、規范并更新三角剖分網格；a、定位點Pr在三角剖分網格中的位置；b、如果點Pr在三角形△PiPjPk內部，則連接點Pr和三角形△PiPjPk的三個頂點，形成三個子三角形；如果Pr在三角形△PiPjPk的一條邊上，則連接Pr和△PiPjPk的第三個頂點，形成兩個新的三角形；c、通過空外接圓準則檢查并規范整個三角形網格。如果第四個頂點在三角形的外接圓之內，那么修正對角線（對角線對調）完成局部優化過程的處理；如果第四個頂點不在三角形的外接圓之內，那么已經符合空外接圓準則不做任何處理。d、做一次對角線修正，更新一次三角形網格。

（4）待點集P中所有的點都插入完全之后移出初始大三角形。

Delaunay四面體剖分與三角形剖分類似，其規范網格的原則為空外接球準則，即不多于4點共球。

1.2 數據插值

1.2.1 判斷待插值點的位置

已知三角形（四面體）的三（四）個頂點坐標，則三角形（四面體）的面積（體積）可利用公式求出。判斷待插值點所在的位置。

可利用待插值點與某個三角形（四面體）形成的三個子三角形（四面體）的面積（體積）之和是否等于此三角形的面積（體積）來判斷點在網格中的那個三角形（四面體）中。查找過程中可使用二分法來減少計算量。

1.2.2 插值算法

三角形、四面體插值算法通常有3種：最近鄰插值(nearest neighbor interpolation)、自然鄰點插值(Nearest neighbor interpolation)、線性插值(Linear interpolation)。

最近鄰插值：將距離待插值點最近且已知的數據點值作為插值結果。

自然鄰點插值：與待插值點泰森多邊形相交的泰森多邊形中的樣本點被用來參與插值，它們對待插值點的影響權重和它們所處泰森多邊形與待插值點新生成的泰森多邊形相交的面積成正比。

線性插值：待插值點與三角形（四面體）各個頂點的連線將剖分好的三角形（四面體）分為2或3個新的小三角形（四面體），以新的小三角形（四面體）的面積（體積）為權重來進行插值。

作者編寫了平面三角形網格插值的三種方法并進行了對比，數據來源于野外實際的瞬變電磁，如圖2.4所示，三種插值方法中：線性插值的效果最好，最近鄰插值次之，自然鄰點插值最差。

針對四面體插值，作者只實現了線性插值，插值網格剖分如圖4所示。

圖4 網格化及插值點圖示

2 實際應用

2.1 平面數據插值

在野外數據采集中電極安放好后，如果某一個電極的接地條件不好或線頭接觸不良，對于供電回路，直接影響著供電電流的大小，從而影響著電位差的測量精度；對于測量回路，則會產生數據出現數據不穩定甚至虛假點。單個電極的接地條件不好或線頭接觸不良，在整個電阻率擬斷面上反應出“八”、“”或“/”字型的假異常，如圖5所示。采用半二階差分、曲率公式導出的濾波公式對“壞點”進行識別[6]，如圖6所示，之后對“壞點”進行插值來處理野外生產過程中所造成的數據缺失，如圖7所示。

圖5 原始數據視電阻率圖

圖6 虛假點位置

圖7 數據插值后視電阻率圖

2.2三維數據插值

實際工作中的瞬變電磁方法反演得到的數據為一系列離散的三維坐標和電阻率，其剖面的視電阻率圖無法直觀的反應地層、斷裂走向，經過三維插值做“切片”之后，如圖8，可以清晰的推斷出地下地質體的走向傾向。

圖8 三維數據切片