插值算法在物探數(shù)據(jù)處理中的應(yīng)用

冉根領(lǐng)，王 磊

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

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

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

到目前為止，國內(nèi)外插值研究的相關(guān)著作多達(dá)4000多項(xiàng)，許多文獻(xiàn)對離散數(shù)據(jù)的插值做過深入研究，并且提出許多行之有效的算法。作者沿著前人的腳步，對物探數(shù)據(jù)處理中的插值和物探解釋中三維成像進(jìn)行了一些研究，編寫程序?qū)崿F(xiàn)插值算法并應(yīng)用于高密度數(shù)據(jù)處理和物探數(shù)據(jù)的三維地質(zhì)解釋中。

1 方法原理

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

Voronoi圖的簡單類比理解就是：將平面看做一水平面、數(shù)據(jù)點(diǎn)的集合看作n個水面上的振動源點(diǎn)，源點(diǎn)振動水波開始以相同的速度向四面?zhèn)鞑ィ煌c(diǎn)振動所引起的水波相遇時著振動停止，最后水波紋相遇的地方所形成的線構(gòu)成的網(wǎng)格就是Voronoi圖。如圖1。

圖1 Voronoi圖

圖2 Voronoi圖與Delaunay三角剖分

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

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

Delaunay三角剖分所具備的優(yōu)異特性[4]：

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

（2）唯一性：不論從區(qū)域何處開始構(gòu)建，最終都將得到一致的結(jié)果。

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

（4）最規(guī)則：如果將三角網(wǎng)中的每個三角形的最小角進(jìn)行升序排列，則Delaunay三角網(wǎng)的排列得到的數(shù)值最大。

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

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

同理Delaunay四面體剖分也具有相應(yīng)的性質(zhì)。利用Delaunay三角剖分、四面體剖分來進(jìn)行離散數(shù)據(jù)的網(wǎng)格化可獲得穩(wěn)定、唯一的網(wǎng)格。

圖3 Delaunay四面體剖分

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

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

（2）隨機(jī)的對點(diǎn)集P中的數(shù)據(jù)點(diǎn)進(jìn)行排序；

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

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

Delaunay四面體剖分與三角形剖分類似，其規(guī)范網(wǎng)格的原則為空外接球準(zhǔn)則，即不多于4點(diǎn)共球。

1.2 數(shù)據(jù)插值

1.2.1 判斷待插值點(diǎn)的位置

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

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

1.2.2 插值算法

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

最近鄰插值：將距離待插值點(diǎn)最近且已知的數(shù)據(jù)點(diǎn)值作為插值結(jié)果。

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

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

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

針對四面體插值，作者只實(shí)現(xiàn)了線性插值，插值網(wǎng)格剖分如圖4所示。

圖4 網(wǎng)格化及插值點(diǎn)圖示

2 實(shí)際應(yīng)用

2.1 平面數(shù)據(jù)插值

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

圖5 原始數(shù)據(jù)視電阻率圖

圖6 虛假點(diǎn)位置

圖7 數(shù)據(jù)插值后視電阻率圖

2.2三維數(shù)據(jù)插值

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

圖8 三維數(shù)據(jù)切片