胥凱雄, 劉海飛, 柳建新

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083;2.有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，長沙 410083)

0 引言

受環(huán)境、人員等因素影響，采集到的實(shí)測數(shù)據(jù)經(jīng)常是不均勻分布的，因此必須對數(shù)據(jù)進(jìn)行相應(yīng)地網(wǎng)格化，以便于后期解釋。常用的網(wǎng)格化方法包括最小曲率法、基于三角網(wǎng)的線性插值法、克里格法、多元二次函數(shù)插值法等。最小曲率法[1]是物探位場數(shù)據(jù)的一種常用網(wǎng)格化方法，其速度較快，精度高，但在稀疏控制點(diǎn)的網(wǎng)格化容易出現(xiàn)吉布斯現(xiàn)象。基于三角網(wǎng)的線性插值法簡單易行，計算速度快，但其網(wǎng)格化結(jié)果不光滑，精度不夠，而且計算可能出現(xiàn)不穩(wěn)定?？死锔穹ɡ米儾詈瘮?shù)動態(tài)地決定變量的數(shù)值，是一種在最小估計方差意義上的最佳無偏內(nèi)插法。其網(wǎng)格化精度的高低，依賴于優(yōu)化準(zhǔn)則函數(shù)擬合的好壞，若擬合不好，則其結(jié)果的誤差較大。多元二次函數(shù)法是Hardy[2]提出的一種徑向基函數(shù)插值法，易實(shí)現(xiàn)，網(wǎng)格化精度高。芮小平等[3]將這種方法應(yīng)用到地質(zhì)數(shù)據(jù)解釋中，效果良好。但由于利用全局的已知點(diǎn)來估計每個待插點(diǎn)，而且計算中往往要求解大型線性方程組，故計算時間較慢。劉海飛等[4]提出了方位加權(quán)多元二次函數(shù)曲面插值法，通過八方位搜索取點(diǎn)來提高待插點(diǎn)的特征值，使得補(bǔ)插網(wǎng)格點(diǎn)的真實(shí)值可信，提高了計算效率。筆者基于前人[5-7]對surfer的二次開發(fā)研究，利用MFC和Surfer的結(jié)合，實(shí)現(xiàn)了交互設(shè)計和方位加權(quán)插值算法的等值線圖繪制。通過將可視化平臺與Surfer的簡單結(jié)合，直接進(jìn)行批處理數(shù)據(jù)，為科學(xué)分析和實(shí)際應(yīng)用提供一定的便利。

1 方法原理

1.1 按方位取點(diǎn)

我們采用方位加權(quán)二次插值，即以網(wǎng)格點(diǎn)為中心，把網(wǎng)格區(qū)域分為4個象限，再將每個象限分成若干區(qū)域，從每個區(qū)域中取離這個點(diǎn)最近一點(diǎn)來作加權(quán)計算，求出屬性值(圖1)。用方位加權(quán)時還需選擇象限等分?jǐn)?shù)，即每一個象限分成幾個區(qū)域。而由于實(shí)際測量的原始值是離散的點(diǎn)，呈非均勻分布，無法給出合適的搜索半徑。對此，可將插值區(qū)域的總面積等份分成10份，取每份面積作為初始搜索半徑R，如果搜索的點(diǎn)達(dá)不到設(shè)定的個數(shù)，則以1.5R為半徑繼續(xù)搜索；若已經(jīng)滿足或到達(dá)搜索邊界時，搜索中止；如果每個等分角內(nèi)搜索的已知點(diǎn)數(shù)多于設(shè)定數(shù)，則選取最近點(diǎn)進(jìn)行計算。這樣以來，在增強(qiáng)網(wǎng)格節(jié)點(diǎn)的局部自適應(yīng)性的同時，保持局部補(bǔ)插點(diǎn)的真實(shí)性。

1.2 實(shí)現(xiàn)過程

多重二次曲面插值法是基于二次曲面方程，構(gòu)建區(qū)域內(nèi)一系列不規(guī)則的連續(xù)曲面。這里采用自適性較好的圓錐方程(x2+y2)tan2α=z作為曲面函數(shù)。方程中z為圓錐頂點(diǎn)的標(biāo)高；tanα為圓錐曲面的斜率，相當(dāng)于權(quán)系數(shù)。根據(jù)已知的數(shù)據(jù)點(diǎn)，由此得到多重二次曲面函數(shù)的方程組式(1)。

(1)

式中：i= 1，2，…n；cj代表第j個屬性點(diǎn)相對第i屬性點(diǎn)所張圓錐面的斜率；(xi，yi)和(xj，yj)分別為已知點(diǎn)i和J的坐標(biāo)；zi為已知點(diǎn)的屬性值；n代表局部已知點(diǎn)的個數(shù)。若令

CT=|c1c2…cn|

ZT=|z1z2…zn|

則可將式(1)改寫為矩陣形式為式(2)。

AC=Z

(2)

式(1)中:斜率系數(shù)c為未知數(shù)，可以通過列選主元法求解得出。然后將解向量cj及待插點(diǎn)(xp，yp)代入式(3)。

圖1 方位取點(diǎn)示意圖Fig.1 The sketch of searching data points on orientation

(3)

通過求和可以計算出待插點(diǎn)的屬性值zp。這樣以來，平面區(qū)域內(nèi)所有網(wǎng)格節(jié)點(diǎn)的屬性值即構(gòu)建的圓錐頂點(diǎn)便可以得到，由此所有實(shí)測數(shù)據(jù)的網(wǎng)格化處理完成。

2 可視化軟件設(shè)計

Visual C++平臺是基于Windows環(huán)境下應(yīng)用最廣泛的可視化集成開發(fā)工具，而MFC就是一個龐大的類庫，可以實(shí)現(xiàn)用戶界面的設(shè)計、文件的操作及數(shù)據(jù)庫訪問等強(qiáng)大的功能。本軟件在VS2010環(huán)境下，依托MFC對話框，用c++語言進(jìn)行編程，設(shè)計圖形用戶界面，可視化實(shí)現(xiàn)方位加權(quán)多重二次插值算法。

圖2 軟件思路流程圖Fig.2 The flowchart of software ideas

圖3 初始化的可交互界面圖Fig.3 Initial interactive interface diagram

如圖2所示，是該軟件的設(shè)計思路流程圖。首先，運(yùn)行得到初始化的應(yīng)用程序，點(diǎn)擊讀取數(shù)據(jù)文件按鈕后，可以快速地顯示出文本中XYZ三列數(shù)據(jù)的信息，接著對其進(jìn)行二維網(wǎng)格，最終生成并保存為Surfer的可讀文件*grd格式。由于Surfer軟件自帶的mapping選項(xiàng)包括了數(shù)據(jù)的標(biāo)注，圖形的拾取等強(qiáng)大功能，再調(diào)用其進(jìn)行等值線的繪制和設(shè)置。

圖4 區(qū)域網(wǎng)格化可視圖Fig.4 Regional gridded view

MFC下初始化的應(yīng)用程序圖(圖3)，一共可分為三個模塊：第一模塊是功能區(qū)，具體內(nèi)容包括讀取實(shí)測數(shù)據(jù)，二維網(wǎng)格化以及最終生成可讀文件。通過這一結(jié)合，直接調(diào)用功能按鈕，實(shí)現(xiàn)人機(jī)交互的響應(yīng)，不必在后臺撰寫大量代碼，便于用戶進(jìn)行批處理數(shù)據(jù)文件，解釋地質(zhì)地理信息；第二模塊為顯示區(qū)，用戶讀取數(shù)據(jù)文件后，在界面端可以顯示最值、數(shù)據(jù)個數(shù)等源數(shù)據(jù)信息。顯示區(qū)的編輯框均含成員變量，方便用戶自定義網(wǎng)格化參數(shù)，包括方向上的節(jié)點(diǎn)數(shù)、間距、最值等，實(shí)現(xiàn)所需區(qū)域的網(wǎng)格插值，然后利用UpdateData(true)函數(shù)將外部值返回到成員變量內(nèi)部；第三模塊為繪圖區(qū)，設(shè)計測網(wǎng)時，也方便網(wǎng)格化時，了解區(qū)域測點(diǎn)的分布稀疏程度，設(shè)置合適的網(wǎng)格化參數(shù)，通過界面坐標(biāo)轉(zhuǎn)換，在繪圖區(qū)顯示測點(diǎn)的分布及網(wǎng)格區(qū)域。

2.1 數(shù)據(jù)導(dǎo)入

打開對話框，設(shè)置好相應(yīng)的參數(shù)，當(dāng)用戶點(diǎn)擊控件按鈕時，DoModal函數(shù)進(jìn)行消息映射，負(fù)責(zé)顯示這個模態(tài)的對話框，利用malloc分別動態(tài)分配x、y、z坐標(biāo)數(shù)據(jù)的存儲空間，按行讀入數(shù)據(jù)，以空格為界，讀出行中的各個字符串，分別保存到對應(yīng)的數(shù)組當(dāng)中。讀取過程中，如果數(shù)據(jù)較多時，添加進(jìn)度條，便于動態(tài)的響應(yīng)程序。while(！infile.eof())函數(shù)判斷是否到達(dá)文件尾部，若所有數(shù)據(jù)都讀取完畢，則跳出while循環(huán)，進(jìn)度條到達(dá)末尾，關(guān)閉文件。最終統(tǒng)計出數(shù)據(jù)個數(shù)(即行數(shù))和最值等信息。

2.2 離散點(diǎn)顯示

讀取數(shù)據(jù)信息結(jié)束后，為了便于用戶了解區(qū)域網(wǎng)格測點(diǎn)的分布情況，可以將二維坐標(biāo)點(diǎn)通過界面可視化，方便選擇合適的參數(shù)統(tǒng)計和排布區(qū)域資料。數(shù)據(jù)讀取完成后，利用像素點(diǎn)在界面上進(jìn)行可視化操作。換用FillSolidRect函數(shù)，在指定矩形進(jìn)行填充，便于顯示(圖4)。部分代碼如下：

for(int i=0;i

{point.x = mDrawRect.left;

point.y = mDrawRect.top;

if(nxrange !=0) point.x +=(m_pt[i].x - minvalueX)/nxrange * mwidth;

if(nyrange !=0) point.y +=(mheight - (m_pt[i].y - minvalueY)/nyrange * mheight);

dc.FillSolidRect(point.x-1,point.y-1,2,2,GRB(255,0,0)); }

2.3 二維網(wǎng)格化

利用move to、line to函數(shù)，畫出矩形區(qū)的橫縱軸線，即可視化的區(qū)域網(wǎng)格分布。通過Invalidate實(shí)現(xiàn)Onpaint對窗口的不斷重繪，轉(zhuǎn)換實(shí)際點(diǎn)坐標(biāo)到畫布的像素點(diǎn)，由此畫出二維坐標(biāo)點(diǎn)以及網(wǎng)格區(qū)域，從而完成區(qū)域網(wǎng)格化的界面顯示。為了方便用戶識讀，利用DrawText函數(shù)的指定矩形輸出功能顯示出對應(yīng)的刻度值。然后生成xn個x方向坐標(biāo)，yn個y方向坐標(biāo)，根據(jù)兩個方向的間距，求出區(qū)域內(nèi)各個節(jié)點(diǎn)的坐標(biāo)即待插點(diǎn)的二維坐標(biāo)。節(jié)點(diǎn)坐標(biāo)的最值確定出已知屬性值的最大矩形區(qū)域，定義數(shù)組存放八個方位的屬性值以及屬性值的個數(shù)，按八個方位搜索已知屬性點(diǎn)，計算出這些待插點(diǎn)的z值。

圖5 等值線圖Fig.5 Contour maps(a) 克里格插值法(100*78)；(b) 克里格插值法(200*100)；(c) 多重二次插值法(100*78)；(d)多重二次插值法(200*100)；(e)方位加權(quán)多重二次法(節(jié)點(diǎn)數(shù)100*78)；(f)方位加權(quán)多重二次法(節(jié)點(diǎn)數(shù)200*100)

2.4 網(wǎng)格文件輸出

等值線的繪制需要將數(shù)據(jù)文件格式轉(zhuǎn)化為Surfer可讀的*grd格式文件，方便處理和成圖。ascii碼的直接grd文件格式包括5個文件頭信息，包括行列數(shù)、屬性值等。下面給出它的一般格式：

DSAA

xnyn(x、y方向網(wǎng)格點(diǎn)的個數(shù)xn、yn)

xminxmax(x方向網(wǎng)格點(diǎn)的最小值xmin和最大值xmax)

yminymax(y方向網(wǎng)格點(diǎn)的最小值ymin和最大值ymax)

zminzmax(z方向網(wǎng)格點(diǎn)的最小值zmin和最大值zmax)

grid row 1 (z的值在點(diǎn)(x1,y1),…(xn,y1))

grid row 2 (z的值在點(diǎn)(x1,y2), ,…(xn,y2))

grid row……grid row n (z的值在點(diǎn)(x1,yn), (x2,yn),…(xn,yn))

圖6 克里格網(wǎng)格化后的等值線圖Fig.6 Contour map obtained by krigin interpolation

圖7 多重二次所成等值線圖Fig.7 Contour map obtained by multiple quadric surface algorithm

圖8 方位加權(quán)多重二次曲面法所成等值線圖Fig.8 Contour map obtained by azimuth weighted multiple quadric surface algorithm

通過插值算法求出網(wǎng)格點(diǎn)的屬性值z后，利用VC編程，按行寫入上述信息，即完成了*grd格式的生成。

3 網(wǎng)格化結(jié)果

3.1 算例一

為了驗(yàn)證算法的準(zhǔn)確性，圖5是Surfer軟件中Demogrid數(shù)據(jù)文件，分別經(jīng)過克里格法、多重二次插值法及方位加權(quán)多重二次法插值得到的電阻率等值線圖，內(nèi)部的紅色點(diǎn)為加載的數(shù)據(jù)文件。通過圖5(a)、圖5(c)、圖5(e)對比可知，當(dāng)采用默認(rèn)參數(shù)網(wǎng)格化時(節(jié)點(diǎn)數(shù)為100*78個)，多重二次法和方位加權(quán)多重二次算法整體上與克里格插值保持一致，平滑性較好，無不規(guī)則突變點(diǎn)。由于克里格法是一種無偏估計，故穩(wěn)定性較好。當(dāng)加密網(wǎng)格后(節(jié)點(diǎn)數(shù)變?yōu)?00*100個)，網(wǎng)格化效果未發(fā)生較大改變。而其他兩種方法加密后的網(wǎng)格化效果局部較節(jié)點(diǎn)數(shù)少時相比略有不光滑，但整體性特征未變。

綜上所述，當(dāng)數(shù)據(jù)點(diǎn)較少時，三種方法的網(wǎng)格化效果基本一致，而在使用多重二次法和方位加權(quán)多重插值法進(jìn)行網(wǎng)格化處理時，也要考慮節(jié)點(diǎn)數(shù)對局部特征的影響。

3.2 算例二

圖6、圖7、圖8分別是野外實(shí)測數(shù)據(jù)經(jīng)過克里格插值、多重二次法及方位加權(quán)多重二次插值所成的電阻率分布等值線圖，內(nèi)部紅色點(diǎn)顯示了測點(diǎn)的分布情況，共6 571個數(shù)據(jù)點(diǎn)。

通過對比圖6～圖8可以發(fā)現(xiàn)，三種方法都能較好地擬合地下的總體特征。圖6平滑程度較好，但是在網(wǎng)格化后，判斷等值線極值圈附近時存在有負(fù)值的情況。而圖7、圖8網(wǎng)格化效果沒有太大差別，兩者較圖6在局部特征上不光滑，但是沒有出現(xiàn)等值線負(fù)值的情況，較好地解決了等值線極值圈的判斷及插值。而在本次實(shí)際數(shù)據(jù)點(diǎn)較多時，三種方法的網(wǎng)格化效率有所差別。方位加權(quán)算法用時5 s，略快于克里格法插值的8 s，而不加權(quán)的多重二次法由于將所有的點(diǎn)都參與計算，用時多達(dá)20 s左右，計算效率遠(yuǎn)遠(yuǎn)低于其他兩種方法。

綜上所述，當(dāng)數(shù)據(jù)較多時，采用方位加權(quán)多重二次法，在提高計算效率的同時，保證了較高的節(jié)點(diǎn)插值精度，總體網(wǎng)格化質(zhì)量較好，與實(shí)際情況相符，可以推廣應(yīng)用。

4 結(jié)語

筆者采用方位加權(quán)多重二次插值算法，實(shí)現(xiàn)可視化界面后，對數(shù)據(jù)進(jìn)行網(wǎng)格化，求出網(wǎng)格節(jié)點(diǎn)的屬性值。由于等值線圖可以直觀顯示地質(zhì)體的物理信息，通過調(diào)用 Surfer 的Shapes 對象生成并顯示該實(shí)測數(shù)據(jù)分布區(qū)域的等值線圖，與其內(nèi)部插值方法進(jìn)行對比。利用MFC和Surfer的結(jié)合，一方面實(shí)現(xiàn)了交互式設(shè)計，便于用戶批處理數(shù)據(jù)，提高了數(shù)據(jù)分析和圖形繪制的效率，另一方面可以快速簡單的完成等值線圖繪制，成圖效果較為理想，保證了不規(guī)則數(shù)據(jù)的網(wǎng)格化效果和實(shí)際情況的要求。