李明超，繆正建，劉 菲，王 剛

(1．天津大學水利工程仿真與安全國家重點實驗室，天津 300072;2．中國水電顧問集團成都勘測設計研究院，成都 610072)

1 前言

地質(zhì)勘探結果大多反映在一些離散不規(guī)則分布的數(shù)據(jù)點上，為了通過這些離散數(shù)據(jù)建立起區(qū)域性連續(xù)的整體模型，需要利用插值和逼近的曲面處理方法［1，2］。曲面插值(surface interpolation)是嚴格通過給定的數(shù)據(jù)點來構造曲面，并根據(jù)原始數(shù)據(jù)點值來插補空白區(qū)的值;這類方法不改變原始數(shù)據(jù)點值。而曲面逼近(surface approximation)則是利用相對簡單的數(shù)學曲面來近似構造復雜的地學曲面，根據(jù)一定的數(shù)學準則，使所給出的數(shù)學曲面最大限度地逼近地質(zhì)曲面;通過擬合處理的曲面，原始數(shù)據(jù)點一般有所改變，所以曲面逼近的結果往往會取得平滑的效果。曲面的插值與逼近統(tǒng)稱為曲面擬合(surface fitting)［3］。

在地質(zhì)曲面構造中運用較多的插值和逼近方法包括按近點距離加權平均法、按方位取點加權法、雙線性插值法、移動曲面插值法、二元三點插值法、Kriging插值法和三次樣條函數(shù)擬合法、趨勢面擬合法、加權最小二乘擬合法等［4～9］。插值方法能保證構造曲面嚴格通過原始數(shù)據(jù)點集合，適用于均勻分布的數(shù)據(jù)點，但外推能力和唯一性較差［1］;地質(zhì)工程師一般傾向于使用插值方法，這在二維環(huán)境下能夠有效地操作，但在三維環(huán)境下解決基于分布不均勻原始數(shù)據(jù)的復雜地質(zhì)曲面構造時就會遇到較大的困難［6，9］?；诒平枷氲那鏀M合方法對原始數(shù)據(jù)點的分布沒有任何要求，且外推能力和唯一性均較強，并能夠保證曲面對原始數(shù)據(jù)點集合具有最佳逼近效果和很好的光順性［10］;但是直接采用逼近方法應用于分布不均勻原始數(shù)據(jù)點的復雜曲面擬合也存在兩點困難［3］，一是如何確定擬合曲面的控制網(wǎng)格頂點個數(shù)保證其滿足精度要求，二是如何保證所構造的逼近曲面的形狀能夠滿足實際需要。

由于自然界地質(zhì)結構的復雜性，各種地質(zhì)曲面的擬合構造方法始終是地質(zhì)學家和計算機科學家關注的熱點問題，目前仍然缺乏完善的處理方法［11］。文章以水利水電工程為例，在分析多源地質(zhì)數(shù)據(jù)特征的基礎上，針對曲面插值和逼近方法各自的優(yōu)缺點，提出了插值—逼近相結合的復雜地質(zhì)曲面擬合方法，以滿足工程區(qū)域三維地質(zhì)建模精度高且存儲量小的要求［12］。

2 多源地質(zhì)數(shù)據(jù)

在水利水電工程中，通過多種地質(zhì)勘測方式如地質(zhì)測繪、遙感、地質(zhì)勘探等，得到了各種豐富多樣的地質(zhì)數(shù)據(jù)，如地質(zhì)點資料，遙感圖像，鉆孔、平硐及物探信息等，如圖1所示。同時，對這些原始數(shù)據(jù)加以解譯分析，可獲得相應的地層界線、斷層、褶皺等構造跡線以及平面構造地質(zhì)圖、不同位置的橫/縱剖面圖和不同高程的平切圖等。

圖1 多源地質(zhì)數(shù)據(jù)示例Fig．1 Examples of multi-source geological data

這些來源不同的數(shù)據(jù)在精度、分辨率、數(shù)量、質(zhì)量等方面都存在較大的差異，為了使所有有效數(shù)據(jù)一起成為三維地質(zhì)曲面擬合可利用的、可靠的、一致的信息，應采用統(tǒng)一的數(shù)據(jù)結構進行管理。因此，采用面向?qū)ο蠓诸惡拖到y(tǒng)的思想，將直接可用的鉆孔、平硐信息與解譯分析得到的剖面數(shù)據(jù)進行耦合處理［11］，按不同地層、斷層等對象進行自動分層和三維化，最后統(tǒng)一以點集合的數(shù)據(jù)結構進行存儲，供三維地質(zhì)曲面擬合使用。圖2為某水電工程一個地質(zhì)界面的點集合數(shù)據(jù)實例，圖中數(shù)據(jù)比較集中且分布均勻的區(qū)域是工程主要樞紐建筑物布置區(qū)域，而在其相關的周邊地區(qū)勘探分布較少且離散。

3 地質(zhì)曲面的插值—逼近擬合構造

圖2 地層界面的多源數(shù)據(jù)集合Fig．2 Data set of a horizon boundary surface

根據(jù)上述分析可知，原始的工程地質(zhì)勘探數(shù)據(jù)往往分布范圍較大且不均勻，單獨采用插值方法或者逼近方法都會遇到較大的困難，因此文章提出插值與逼近相結合的方法。在構造復雜地質(zhì)曲面的過程中，充分考慮原始數(shù)據(jù)點的分布特點和曲面構造的精度要求，對于集中且均勻分布的原始數(shù)據(jù)，采用插值方法，使曲面嚴格通過這些數(shù)據(jù)點;對于分布離散不均勻的數(shù)據(jù)，若精度要求很高則采用插值方法進行加密，一般情況下采用逼近擬合方法，使曲面在給定精度下充分逼近原始數(shù)據(jù)。在水利水電工程壩區(qū)，)地質(zhì)勘探數(shù)據(jù)集中且分布較均勻，其構造的地質(zhì)曲面精度要求很高，這部分須采用插值方法;而在其相關的周邊地區(qū)，對精度要求不高，采用逼近方法是可行的。

3．1 NURBS 曲面技術

NURBS(Non-uniform rational B-splines，非均勻有理樣條曲線)技術是ISO(International Standard Organized，國際標準化組織)1991年頒布的STEP標準中自由型曲線曲面的唯一表示方法，它是在B樣條函數(shù)基礎上發(fā)展起來的，利用非均勻節(jié)點向量表達式來構造有理B樣條，它對標準的解析圖形(如圓錐曲線、二次曲面、旋轉(zhuǎn)曲面等)和自由型曲線曲面提供了統(tǒng)一的數(shù)學描述。隨著計算機輔助幾何設計(CAGD)的發(fā)展，NURBS技術有了快速的發(fā)展，已被廣泛應用于包括三維地質(zhì)建模在內(nèi)的工程實踐中［13］。雖然如此，但由于NURBS技術的復雜性、地質(zhì)對象的錯綜復雜和所具有的不規(guī)則形態(tài)，阻礙了NURBS在三維工程地質(zhì)中更多的應用。實際上，針對復雜地質(zhì)體形態(tài)的無規(guī)律性變化，選擇NURBS技術進行水利水電工程三維地質(zhì)建模與分析，具有節(jié)省存儲空間、計算機處理簡便易行、數(shù)據(jù)庫管理方便、保證空間唯一性和幾何不變性等優(yōu)點，在地質(zhì)表示方面有很高的應用價值。

給定一組(m+1)(n+1)的網(wǎng)格控制點Pij(0≤i≤m，0≤j≤n)的NURBS曲面可定義為:

式(1)中，wij為相應于控制點Pij的權因子;k、l為階數(shù);Nki(u)(0≤i≤m)和Nlj(v)(0≤j≤n)分別是定義在u、v向節(jié)點矢量U={u0，…，um+k|ui≤ui+1，i=0，…，(m+k － 1)}和 V={v0，…，vn+l|vj≤ vj+1，j=0，…，(n+l－1)}上的 k、l階 B 樣條基函數(shù);Sij(u，v)表示擬合的曲面段，u∈［ui，ui+1］，i= (k － 1)，…，m ，v ∈ ［vj，vj+1］，j=(l－1)，…，n。在實際工程中，一般 k、l取3基本上就可以滿足要求。

3．2 基于NURBS技術的插值—逼近擬合實現(xiàn)

以圖2中給定的數(shù)據(jù)點集合為例，設該地層面的原始數(shù)據(jù)集合為 D={ps，s=1，2，…，m}，其中分布密集均勻的子集為 D1={pi，i=1，2，…，n}，剩余分布不均勻的子集為 D2={pj，j=1，2，…，r1}和 D3={pk，k=1，2，…，r2}，這里 n+r1+r2=m，如圖3(a)所示。則基于NURBS技術的地質(zhì)曲面插值—逼近擬合方法具體實現(xiàn)如下:

1)對于子集 D1，按照蒙皮法(skinning)［14］插值思想構造曲面，若曲面控制點為n，該方法的數(shù)學模型為:

式(2)中，X∈Rn，為曲面控制頂點構成的未知矢量;B∈Rn，為曲面集合D1及邊界條件構成的已知矢量;A∈Rn×Rn，為系數(shù)矩陣。

由式(2)可知，蒙皮法是一種基于插值思想的曲面構造方法，能夠保證曲面嚴格通過集合D1中的所有數(shù)據(jù)點。該方法使曲面在進行插值時順序地通過一系列截面曲線，從而很好地生成復雜的自由曲面。在NURBS蒙皮過程中，首先要求各截面曲線的拓撲性質(zhì)應該一致，即都是開曲線或都是閉曲線;其次在蒙皮插值過程中，插值的方式將影響最終的結果，根據(jù)工程地質(zhì)實際需要，一般選擇線性插值。圖3(b)為子集D1通過蒙皮插值得到的曲面S1。

2)對于子集D2和D3，則直接利用NURBS曲面技術采用反算法［15］進行擬合。為了滿足曲面的邊界約束條件(鄰接曲面邊界和區(qū)域邊界約束)，這里直接提取已構建曲面S1與子集D2、D3相鄰的邊界線加入子集D2、D3，同時根據(jù)研究區(qū)域?qū)⒃紨?shù)據(jù)轉(zhuǎn)化為一系列u、v方向的曲線矢量，利用反算法求出各曲線矢量上的控制點。進而，基于設定精度和控制點數(shù)據(jù)運用分段NURBS函數(shù)進行逼近擬合，獲得按精度充分逼近子集D2、D3的擬合曲面S2、S3，如圖3(c)所示。這里曲面控制點的權值均設置為1(可進行調(diào)整)，曲面邊界公差為0．01，內(nèi)部公差為0．1。

3)合并曲面S1和曲面S2、S3，由于構建S2、S3曲面的右邊界和左邊界數(shù)據(jù)均來自S1曲面，因此兩者能夠達到實際精度要求的縫合，最終得到完整的地質(zhì)擬合曲面S，如圖3(d)所示，該圖對完整曲面進行了渲染表達，同時顯示了相應的原始數(shù)據(jù)集，以便進行對比調(diào)整。

圖3 基于NURBS技術的地質(zhì)曲面插值—逼近擬合構造Fig．3 Interpolation-approximation fitting construction of geological surface based on NURBS technique

3．3 曲面分析與調(diào)整

在完成上述地質(zhì)曲面的構造后，需要對其進行檢查分析與調(diào)整，以滿足實際精度要求和建模需要。這主要從以下兩個方面進行:

1)地質(zhì)結構合理性和曲面幾何性檢查。即檢查所擬合的地質(zhì)結構面整體變化趨勢是否合理、在幾何結構連續(xù)性及拓撲關系上是否正確。若發(fā)現(xiàn)不合理或錯誤之處，可快速方便地對其進行局部調(diào)整(NURBS曲面局部調(diào)整不會影響其他部位)或重新構造。

2)原始數(shù)據(jù)的精度分析。將所有原始數(shù)據(jù)與構造曲面進行偏差分析，以原始數(shù)據(jù)點與曲面間距離為變量，求其標準偏差σ(m):

式(3)中，n為數(shù)據(jù)點數(shù);di為數(shù)據(jù)點i與構造曲面間的距離。

若σ小于控制誤差ε，則曲面滿足要求，否則需要對曲面偏差較大處按原始數(shù)據(jù)進行局部調(diào)整，直到達到誤差要求。對于上述曲面，其原始數(shù)據(jù)點數(shù)n=597，控制誤差 ε =0．5 m，根據(jù)式(3)求得標準偏差值σ=0．17 m，這主要是由逼近擬合曲面引起的誤差，小于控制誤差，精度已足以滿足實際需要。

4 結語

文章針對工程地質(zhì)多源數(shù)據(jù)的分布和地質(zhì)曲面的特征，充分考慮了地質(zhì)精度要求、曲面連續(xù)性和數(shù)據(jù)存儲量等方面的均衡，提出并實現(xiàn)了基于NURBS技術的復雜地質(zhì)曲面插值—逼近擬合構造方法。該方法對于工程關鍵區(qū)域集中且均勻分布的原始數(shù)據(jù)，采用NURBS蒙皮插值方法，使曲面嚴格通過這些數(shù)據(jù)點;對于周邊區(qū)域分布離散的數(shù)據(jù)，采用逼近擬合方法，使曲面在給定精度下充分逼近原始數(shù)據(jù);并可通過地質(zhì)結構合理性檢查、曲面幾何性和原始數(shù)據(jù)精度分析對構造曲面進行檢查分析和調(diào)整。實例研究表明，該方法所構造的地質(zhì)曲面能滿足地質(zhì)工程師的實際需要，并能為進一步的三維地質(zhì)建模提供基礎，效果良好。

［1］黃健全，羅明高，胡雪濤．實用計算機地質(zhì)制圖［M］．北京:地質(zhì)出版社，1998．

［2］Horsman J，Bethel W．Methods of constructing a 3D geological model from scatter data［C］//Proceedings of AVS 95，http://www － vis．lbl．gov/Publications/1995/Horsman － AVS95．ps，1995．

［3］彭芳瑜，周云飛，周 濟．基于插值與逼近的復雜曲面擬合［J］．工程圖學學報，2002，23(4):87－96．

［4］Fernández O．Obtaining a best fitting plane through 3D georeferenced data［J］．Journal of Structural Geology，2005，27(5):855－858

［5］吳春發(fā)，李 星．地質(zhì)模擬中數(shù)據(jù)插值方法的應用［J］．地球信息科學，2004，6(2):50－52．

［6］周小文，付 暉，吳昌瑜．地層特性隨機場插值方法應用研究［J］．巖土力學，2005，26(2):221－224．

［7］李曉軍，王長虹，朱合華．Kriging插值方法在地層模型生成中的應用［J］．巖土力學，2009，30(1):157－162．

［8］Wang D，Wang H T，Xi J H，et al．Three-dimensional geological modeling with discrete smooth interpolation algorithm［C］//2010 International Conference on Remote Sensing，IEEE，2010:125－129．

［9］Sirakov N M，Granado I，Muge F H．Interpolation approach for 3D smooth reconstruction of subsurface objects［J］．Computers ＆Geosciences，2002，28(8):877 －885．

［10］秦緒佳，王 青，鮑虎軍．基于散亂點的增量式曲面逼近［J］．計算機輔助設計與圖形學學報，2006，18(9):1408－1413．

［11］徐 華，武 強．復雜地質(zhì)體中多值面的網(wǎng)格生成算法［J］．計算機輔助設計與圖形學學報，2002，14(7):609－612．

［12］鐘登華，李明超．水利水電工程地質(zhì)三維建模與分析理論及實踐［M］．北京:中國水利水電出版社，2006．

［13］Zhong D H，Li M C，Song L G，et al．Enhanced NURBS modeling and visualization for large 3D geoengineering applications:an example from the Jinping first-level hydropower engineering project，China［J］．Computers＆ Geosciences，2006，32(9):1270－1282．

［14］Piegl L，Tiller W．Algorithm for approximate NURBS skinning［J］．Computer－Aided Design，1996，28(9):699－706．

［15］施法中．計算機輔助幾何設計與非均勻有理 B樣條(CAGD＆NURBS)［M］．北京:高等教育出版社，2001．