利用MFC和Surfer實現(xiàn)方位加權(quán)多重二次曲面插值

胥凱雄, 劉海飛, 柳建新

(1.中南大學 地球科學與信息物理學院，長沙 410083;2.有色資源與地質(zhì)災害探查湖南省重點實驗室，長沙 410083)

0 引言

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

1 方法原理

1.1 按方位取點

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

1.2 實現(xiàn)過程

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

(1)

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

CT=|c1c2…cn|

ZT=|z1z2…zn|

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

AC=Z

(2)

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

圖1 方位取點示意圖Fig.1 The sketch of searching data points on orientation

(3)

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

2 可視化軟件設計

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

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

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

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

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

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

2.1 數(shù)據(jù)導入

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

2.2 離散點顯示

讀取數(shù)據(jù)信息結(jié)束后，為了便于用戶了解區(qū)域網(wǎng)格測點的分布情況，可以將二維坐標點通過界面可視化，方便選擇合適的參數(shù)統(tǒng)計和排布區(qū)域資料。數(shù)據(jù)讀取完成后，利用像素點在界面上進行可視化操作。換用FillSolidRect函數(shù)，在指定矩形進行填充，便于顯示(圖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實現(xiàn)Onpaint對窗口的不斷重繪，轉(zhuǎn)換實際點坐標到畫布的像素點，由此畫出二維坐標點以及網(wǎng)格區(qū)域，從而完成區(qū)域網(wǎng)格化的界面顯示。為了方便用戶識讀，利用DrawText函數(shù)的指定矩形輸出功能顯示出對應的刻度值。然后生成xn個x方向坐標，yn個y方向坐標，根據(jù)兩個方向的間距，求出區(qū)域內(nèi)各個節(jié)點的坐標即待插點的二維坐標。節(jié)點坐標的最值確定出已知屬性值的最大矩形區(qū)域，定義數(shù)組存放八個方位的屬性值以及屬性值的個數(shù)，按八個方位搜索已知屬性點，計算出這些待插點的z值。

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

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

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

DSAA

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

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

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

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

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

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

grid row……grid row n (z的值在點(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)格點的屬性值z后，利用VC編程，按行寫入上述信息，即完成了*grd格式的生成。

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

3.1 算例一

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

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

3.2 算例二

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

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

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

4 結(jié)語

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