利用MFC和Surfer實現方位加權多重二次曲面插值

胥凱雄, 劉海飛, 柳建新

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

0 引言

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

1 方法原理

1.1 按方位取點

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

1.2 實現過程

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

(1)

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

CT=|c1c2…cn|

ZT=|z1z2…zn|

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

AC=Z

(2)

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

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

(3)

通過求和可以計算出待插點的屬性值zp。這樣以來，平面區域內所有網格節點的屬性值即構建的圓錐頂點便可以得到，由此所有實測數據的網格化處理完成。

2 可視化軟件設計

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

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

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

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

圖4 區域網格化可視圖Fig.4 Regional gridded view

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

2.1 數據導入

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

2.2 離散點顯示

讀取數據信息結束后，為了便于用戶了解區域網格測點的分布情況，可以將二維坐標點通過界面可視化，方便選擇合適的參數統計和排布區域資料。數據讀取完成后，利用像素點在界面上進行可視化操作。換用FillSolidRect函數，在指定矩形進行填充，便于顯示(圖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 二維網格化

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

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

2.4 網格文件輸出

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

DSAA

xnyn(x、y方向網格點的個數xn、yn)

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

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

zminzmax(z方向網格點的最小值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 克里格網格化后的等值線圖Fig.6 Contour map obtained by krigin interpolation

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

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

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

3 網格化結果

3.1 算例一

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

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

3.2 算例二

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

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

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

4 結語

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