Delaunay三角剖分插值算法在MT成圖中的應(yīng)用①

楊利容，簡興祥

（1.成都理工大學(xué)，四川 成都 610059；2.四川省核工業(yè)地質(zhì)局，四川 成都 610021）

0 引言

在大地電磁（MT）反演迭代過程中，需要實時生成多種圖件提供給反演解釋人員，以動態(tài)了解反演擬合的情況，如反演結(jié)果二維等值線圖、反射系數(shù)二維等值線圖、實測數(shù)據(jù)二維等值線圖以及當(dāng)前模型正演的擬合結(jié)果二維等值線圖等。但由于數(shù)據(jù)采集技術(shù)、工區(qū)客觀條件上的種種限制以及計算機的計算能力等原因，地球物理學(xué)中的大多數(shù)數(shù)據(jù)都是不規(guī)則的離散數(shù)據(jù)。繪制以上各種數(shù)據(jù)的等值線圖時，首先需要對離散數(shù)據(jù)進行網(wǎng)格化處理，將稀疏的、不規(guī)則分布的數(shù)據(jù)插值加密為規(guī)則分布的數(shù)據(jù)，以適應(yīng)成像繪圖的需要。在MT反演迭代過程中對數(shù)據(jù)進行網(wǎng)格化處理必然對插值算法的可靠性和實時性要求比較高。網(wǎng)格化方法的好壞不僅直接影響到網(wǎng)格化數(shù)據(jù)的質(zhì)量、精度和可信程度，而且還將進一步影響到數(shù)據(jù)解釋處理圖件的質(zhì)量、效果和可靠性。為此有必要研究一種可靠性高、精度能滿足要求、針對大數(shù)據(jù)量快速實時的二維網(wǎng)格化方法。

目前在地球物理領(lǐng)域普遍使用的網(wǎng)格化方法主要有最小曲率插值法（minimum curvature splines）、Kriging插值法、加權(quán)反距離插值法等［1］。各種方法各有優(yōu)缺點。如最小曲率插值法速度快，適合處理數(shù)據(jù)量大的情況，由于其主要考慮曲面的光滑性，不能達到精確的插值結(jié)果，容易超出最大值和最小值的范圍；Kriging插值法根據(jù)原數(shù)據(jù)所隱含的趨勢特征，以區(qū)域化變量理論為基礎(chǔ)，以變差函數(shù)為主要工具，在保證研究對象的估計值滿足無偏性條件和最小方差條件的前提下求得估計值，是一種精度比較高的網(wǎng)格化插值算法，但該方法隨著數(shù)據(jù)量的增大計算速度較慢［2］。所以上述常用的插值方法不能滿足實時成像的要求。基于Delaunay三角剖分的快速插值算法由于采用所有的數(shù)據(jù)點去構(gòu)造最佳Delaunay三角網(wǎng)，原始數(shù)據(jù)點在插值后保持不變，其他待插點只跟它所處的三角形的三個頂點相關(guān)，因而是一種局部插值方法；另一方面，該方法并不進行外插，因此適合于網(wǎng)格化處理帶地形特征的數(shù)據(jù)。由于該方法能夠很好的擬合地形，因此被廣泛的應(yīng)用到數(shù)字高程模擬（DEM）中，并出現(xiàn)了很多算法［11－13］。本文擬實現(xiàn)一種基于 Delaunay三角剖分的線性插值的算法，并將其應(yīng)用到自主開發(fā)MT二維反演軟件快速成像二維網(wǎng)格化處理中。

1 Delaunay三角剖分

如何把一個散點集合剖分成不規(guī)則的三角形網(wǎng)格（Triangular Irregular Network，TIN），這就是散點集的三角剖分問題。在所有生成TIN的方法中，Delaunay三角網(wǎng)最優(yōu)，它盡可能避免了病態(tài)三角形的出現(xiàn)，常常被用來生成TIN。對于給定的初始點集P，有多種三角網(wǎng)剖分方式，而Delaunay三角網(wǎng)有以下特性：其Delaunay三角網(wǎng)是唯一的；三角網(wǎng)的外邊界構(gòu)成了點集P的凸多邊形“外殼”；如果將三角網(wǎng)中的每個三角形的最小角進行升序排列，則Delaunay三角網(wǎng)的排列得到的數(shù)值最大，從這個意義上講，Delaunay三角網(wǎng)是“最接近于規(guī)則化”的三角網(wǎng)。

要滿足Delaunay三角剖分的定義，必須符合兩個重要的準(zhǔn)則：（1）空圓特性，任何一個Delaunay三角形的外接圓的內(nèi)部不能包含其它任何點；（2）最大化最小角特性，每兩個相鄰的三角形構(gòu)成的凸四邊形的對角線，在相互交換后，六個內(nèi)角的最小角不再增大。

理論上為了構(gòu)造Delaunay三角網(wǎng)，Lawson提出局部優(yōu)化過程LOP（Local Optimization Procedure），一般三角網(wǎng)經(jīng)過LOP處理，即可確保成為Delaunay三角網(wǎng)。如圖1所示，先求出包含新插入點p的外接圓的三角形，這種三角形稱為影響三角形，刪除影響三角形的公共邊，將p與全部影響三角形的頂點連接，完成p點在原Delaunay三角形中的插入。

圖1 局部優(yōu)化過程LOP方法示意Fig.1 Sketch of local optimization process method.

Delaunay剖分是一種三角剖分的標(biāo)準(zhǔn)，實現(xiàn)它有多種算法，考慮到編碼簡易性以及算法執(zhí)行的效率，本文對平面點集的Delaunay三角剖分采用以下的算法（圖2）。（1）根據(jù)點集的坐標(biāo)范圍，求出點集的凸多邊形外殼；（2）選擇一個點，構(gòu)造一個初始Delaunay三角網(wǎng)；（3）從集合中選擇一個點，采用局部優(yōu)化過程修改生成新的Delaunay三角網(wǎng)；（4）重復(fù)步驟3，直到所有點都計算完。

2 基于Delaunay三角剖分線性插值算法

基于Delaunay三角剖分的線性插值算法首先將輸入的散點集合構(gòu)建一個Delaunay三角網(wǎng)，然后循環(huán)將待插點加入已存在的Delaunay三角網(wǎng)從而實現(xiàn)插值運算。對于給定的一個Delaunay三角形的三個頂點的值fi（xi，yi），其中i＝1，2，3，待插點（x，y）處的值f（x，y）可以通過下式求得［6］：

其中φi（x，y）是一個2－D的基函數(shù)，它是一個在頂點（xi，yi）處（此時值為1）線性變化到另一個頂點（xj，yj）（j≠i）（此時值為0）的函數(shù)。在實際計算中，可以將式（1）變換為

為了求取系數(shù)c＝ ［c1，c2，c3］T，可以根據(jù)已知的三個頂點fi（xi，yi）｛i＝1，2，3｝來聯(lián)立線性方程組來求得：

圖2 平面上Delaunay三角剖分算法示意圖Fig.2 Schematic diagram of Delaunay triangulations algorithm on the plane.

求解式（3）的線性方程組得出系數(shù)c＝ ［c1，c2，c3］T的值，然后利用式（2）便可求出待插點（x，y）處的值。

基于Delaunay三角剖分線性插值算法實現(xiàn)步驟如下：（1）輸入平面上散點數(shù)據(jù)，構(gòu)建Delaunay三角網(wǎng)；（2）輸入待插點pi（x，y），查找pi所處的Delaunay三角形，計算c1，c2，c3，計算f（x，y）；（3）重復(fù)步驟（2）直到計算完所有待插點；（4）輸出計算結(jié)果用于成圖并結(jié)束程序。

3 應(yīng)用實例

由于本文所實現(xiàn)的方法是基于Delaunay三角剖分的線性插值算法，因此它的主要優(yōu)勢在于插值效率，另一方面由于該算法并不外插，因此它另外一個優(yōu)勢在于便于網(wǎng)格化帶地形特征的數(shù)據(jù)。Surfer中的自然鄰點插值算法（NNI）也是基于Delaunay三角剖分的一種局部插值算法，它們具有較好的可比性，因此本文應(yīng)用一個帶地形的野外實測MT數(shù)據(jù)的反演結(jié)果作為例子，從插值精度上作定性的比較，插值效率上作定量的比較。

3.1 應(yīng)用實例1

圖3中的例子中原始數(shù)據(jù)共有7 960個散點數(shù)據(jù)，分別用本文算法（圖3（a））和Surfer NNI算法（圖3（b））用200×200來進行網(wǎng)格化處理。從圖3中可以看出，本文的算法在插值精度和光滑度上要比Surfer NNI插值算法要略差一些，但是總體形態(tài)即插值結(jié)果基本一致，因此在既要兼顧到插值精度又要滿足實時性的情況下，本文的算法是基本符合要求的。

圖3 實例1本文算法和Surfer NNI算法插值結(jié)果比較Fig.3 Comparing of interpolation results between the algorithm of this paper and Surfer NNI in case 1.

為了定量的說明本文的插值算法的插值速度，對圖3所示的例子分別用不同的網(wǎng)格密度對該數(shù)據(jù)進行網(wǎng)格化處理并和商業(yè)軟件Surfer中速度較快的NNI插值算法進行對比，如表1所示。

表1 本文算法與Surfer NNI插值算法對比

從表1中可以看出，本文算法在速度上要優(yōu)于Surfer NNI算法。同時可以看出，隨著插值網(wǎng)格的密度增大，算法的耗時接近線性增加，因此本文的算法的速度相對NNI算法來講更適合于實時性要求較高的場合。

3.2 應(yīng)用實例2

圖4中的例子中原始數(shù)據(jù)共有720個散點數(shù)據(jù)，分別用本文算法（圖4（a））和Surfer NNI算法（圖4（b））用400×400來進行網(wǎng)格化處理。從圖4中可以看出，本文的算法插值結(jié)果在光滑度上要比Surfer NNI插值算法要略差一些，但是總體形態(tài)即插值結(jié)果基本一致。另外，在對帶凸地形的插值效果上，兩種插值算法也基本上一致，都能很好的模擬地形的起伏，都達到了滿意的效果。但在算法速度上，本文的算法僅耗時0.02s，而Surfer NNI卻花費4.95s。因此在對算法實時性要求較高，但精度基本能滿足要求的情況下本文的算法具有較大的優(yōu)勢。

圖4 實例2本文算法和Surfer NNI插值結(jié)果比較Fig.4 Comparing of interpolation results between the algorithm of this paper and Surfer NNI in case 2.

4 結(jié)論

本文實現(xiàn)的基于Denaulay三角剖分的插值算法穩(wěn)定性好，速度快，并具有保真度高、便于網(wǎng)格化帶地形特征的數(shù)據(jù)等優(yōu)勢，適用于MT等不規(guī)則分布數(shù)據(jù)的實時、快速網(wǎng)格化處理，是地球物理數(shù)據(jù)實時二維網(wǎng)格化最佳的插值方法之一。

［1］Watson D F.Contouring：aguide to the analysis and display of spatial data［M］.［S.l.］：Pergamon Press，1992.

［2］劉兆平，楊進，武煒.地球物理數(shù)據(jù)網(wǎng)格化方法的選?。跩］.物探與化探，2010，34（1）：93－97.

［3］郭良輝，孟小紅，等.地球物理不規(guī)則分布數(shù)據(jù)的空間網(wǎng)格化法［J］.物探與化探，2005，29（5）：438－442.

［4］趙文芳.離散點集Delaunay三角網(wǎng)生成算法改進與軟件開發(fā)［J］.測繪工程，2003，12（04）：22－25.

［5］凌海濱，吳兵.改進的自連接Delaunay三角網(wǎng)生成算法［J］.計算機應(yīng)用，1999，19（12）：10－12.

［6］文偉，楊耀權(quán)，于希寧.用Visual C語言實現(xiàn)的Delaunay三角剖分算法［J］.華北電力大學(xué)學(xué)報，2000，27（4）：54－58.

［7］徐青，常歌，楊力.基于自適應(yīng)分塊的TIN三角網(wǎng)建立算法［J］.中國圖象圖形學(xué)報，2000，5（6）：461－465.

［8］潘榮江，屠長河，孟祥旭，等.基于均勻網(wǎng)格的Delaunay三角網(wǎng)算法在隨機聚合網(wǎng)屏中的應(yīng)用［J］.中國圖象圖形學(xué)報 ，2002，7（5）：495－500.

［9］Q Fan，A Efrat，V Koltun，et al.Hardware－Assisted Natural Neighbor Interpolation［A］∥Proceedings of ALENEX［C］.2005：111－120.

［10］Sambridge M，Braun J，McQueen H.Geophysical parameterization and interpolation of irregular data using natural neighbors［J］.Geophysical Journal International 1995，12（2）：837－857.

［11］Lee D T，Schacheer B J.Two algorithms for constructing a Delaunay triangulation［J］.International Journal of Computer and Information Science，1980，9（3）：219－242.

［12］Lawson C L.Software for C surface interpolation［M］.New York：Academic Press，1997.

［13］Watson D F.Computing the n－dimension Delaunay tessellation with application to Voronoi polygons［J］.Computer Journal，1981，24（2）：167－172.

［14］徐凱軍，石雙虎，周家惠.三維大地電磁激電效應(yīng)特征研究［J］.西北地震學(xué)報，2009，31（1）：31－34.

［15］李希亮，劉希強，董曉娜，等.高階統(tǒng)計量方法在地球物理學(xué)中的應(yīng)用與展望［J］.西北地震學(xué)報，2010，32（2）：201－205.