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

(1．天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072;2．中國(guó)水電顧問(wèn)集團(tuán)成都勘測(cè)設(shè)計(jì)研究院，成都 610072)

1 前言

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

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

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

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

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

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

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

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

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

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

3．1 NURBS 曲面技術(shù)

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

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

式(1)中，wij為相應(yīng)于控制點(diǎn)Pij的權(quán)因子;k、l為階數(shù);Nki(u)(0≤i≤m)和Nlj(v)(0≤j≤n)分別是定義在u、v向節(jié)點(diǎn)矢量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。在實(shí)際工程中，一般 k、l取3基本上就可以滿(mǎn)足要求。

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

以圖2中給定的數(shù)據(jù)點(diǎn)集合為例，設(shè)該地層面的原始數(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技術(shù)的地質(zhì)曲面插值—逼近擬合方法具體實(shí)現(xiàn)如下:

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

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

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

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

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

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

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

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

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

2)原始數(shù)據(jù)的精度分析。將所有原始數(shù)據(jù)與構(gòu)造曲面進(jìn)行偏差分析，以原始數(shù)據(jù)點(diǎn)與曲面間距離為變量，求其標(biāo)準(zhǔn)偏差σ(m):

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

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

4 結(jié)語(yǔ)

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

［1］黃健全，羅明高，胡雪濤．實(shí)用計(jì)算機(jī)地質(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ì)．基于插值與逼近的復(fù)雜曲面擬合［J］．工程圖學(xué)學(xué)報(bào)，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ù)插值方法的應(yīng)用［J］．地球信息科學(xué)，2004，6(2):50－52．

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

［7］李曉軍，王長(zhǎng)虹，朱合華．Kriging插值方法在地層模型生成中的應(yīng)用［J］．巖土力學(xué)，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］秦緒佳，王 青，鮑虎軍．基于散亂點(diǎn)的增量式曲面逼近［J］．計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2006，18(9):1408－1413．

［11］徐 華，武 強(qiáng)．復(fù)雜地質(zhì)體中多值面的網(wǎng)格生成算法［J］．計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2002，14(7):609－612．

［12］鐘登華，李明超．水利水電工程地質(zhì)三維建模與分析理論及實(shí)踐［M］．北京:中國(guó)水利水電出版社，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］施法中．計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理 B樣條(CAGD＆NURBS)［M］．北京:高等教育出版社，2001．