李明超,繆正建,劉 菲,王 剛
(1.天津大學(xué)水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室,天津 300072;2.中國水電顧問集團(tuán)成都勘測設(shè)計(jì)研究院,成都 610072)
地質(zhì)勘探結(jié)果大多反映在一些離散不規(guī)則分布的數(shù)據(jù)點(diǎn)上,為了通過這些離散數(shù)據(jù)建立起區(qū)域性連續(xù)的整體模型,需要利用插值和逼近的曲面處理方法[1,2]。曲面插值(surface interpolation)是嚴(yán)格通過給定的數(shù)據(jù)點(diǎn)來構(gòu)造曲面,并根據(jù)原始數(shù)據(jù)點(diǎn)值來插補(bǔ)空白區(qū)的值;這類方法不改變原始數(shù)據(jù)點(diǎn)值。而曲面逼近(surface approximation)則是利用相對(duì)簡單的數(shù)學(xué)曲面來近似構(gòu)造復(fù)雜的地學(xué)曲面,根據(jù)一定的數(shù)學(xué)準(zhǔn)則,使所給出的數(shù)學(xué)曲面最大限度地逼近地質(zhì)曲面;通過擬合處理的曲面,原始數(shù)據(jù)點(diǎn)一般有所改變,所以曲面逼近的結(jié)果往往會(huì)取得平滑的效果。曲面的插值與逼近統(tǒng)稱為曲面擬合(surface fitting)[3]。
在地質(zhì)曲面構(gòu)造中運(yùn)用較多的插值和逼近方法包括按近點(diǎn)距離加權(quán)平均法、按方位取點(diǎn)加權(quán)法、雙線性插值法、移動(dòng)曲面插值法、二元三點(diǎn)插值法、Kriging插值法和三次樣條函數(shù)擬合法、趨勢面擬合法、加權(quán)最小二乘擬合法等[4~9]。插值方法能保證構(gòu)造曲面嚴(yán)格通過原始數(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)的分布沒有任何要求,且外推能力和唯一性均較強(qiáng),并能夠保證曲面對(duì)原始數(shù)據(jù)點(diǎn)集合具有最佳逼近效果和很好的光順性[10];但是直接采用逼近方法應(yīng)用于分布不均勻原始數(shù)據(jù)點(diǎn)的復(fù)雜曲面擬合也存在兩點(diǎn)困難[3],一是如何確定擬合曲面的控制網(wǎng)格頂點(diǎn)個(gè)數(shù)保證其滿足精度要求,二是如何保證所構(gòu)造的逼近曲面的形狀能夠滿足實(shí)際需要。
由于自然界地質(zhì)結(jié)構(gòu)的復(fù)雜性,各種地質(zhì)曲面的擬合構(gòu)造方法始終是地質(zhì)學(xué)家和計(jì)算機(jī)科學(xué)家關(guān)注的熱點(diǎn)問題,目前仍然缺乏完善的處理方法[11]。文章以水利水電工程為例,在分析多源地質(zhì)數(shù)據(jù)特征的基礎(chǔ)上,針對(duì)曲面插值和逼近方法各自的優(yōu)缺點(diǎn),提出了插值—逼近相結(jié)合的復(fù)雜地質(zhì)曲面擬合方法,以滿足工程區(qū)域三維地質(zhì)建模精度高且存儲(chǔ)量小的要求[12]。
在水利水電工程中,通過多種地質(zhì)勘測方式如地質(zhì)測繪、遙感、地質(zhì)勘探等,得到了各種豐富多樣的地質(zhì)數(shù)據(jù),如地質(zhì)點(diǎn)資料,遙感圖像,鉆孔、平硐及物探信息等,如圖1所示。同時(shí),對(duì)這些原始數(shù)據(jù)加以解譯分析,可獲得相應(yīng)的地層界線、斷層、褶皺等構(gòu)造跡線以及平面構(gòu)造地質(zhì)圖、不同位置的橫/縱剖面圖和不同高程的平切圖等。
圖1 多源地質(zhì)數(shù)據(jù)示例Fig.1 Examples of multi-source geological data
這些來源不同的數(shù)據(jù)在精度、分辨率、數(shù)量、質(zhì)量等方面都存在較大的差異,為了使所有有效數(shù)據(jù)一起成為三維地質(zhì)曲面擬合可利用的、可靠的、一致的信息,應(yīng)采用統(tǒng)一的數(shù)據(jù)結(jié)構(gòu)進(jìn)行管理。因此,采用面向?qū)ο蠓诸惡拖到y(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ū)勘探分布較少且離散。
圖2 地層界面的多源數(shù)據(jù)集合Fig.2 Data set of a horizon boundary surface
根據(jù)上述分析可知,原始的工程地質(zhì)勘探數(shù)據(jù)往往分布范圍較大且不均勻,單獨(dú)采用插值方法或者逼近方法都會(huì)遇到較大的困難,因此文章提出插值與逼近相結(jié)合的方法。在構(gòu)造復(fù)雜地質(zhì)曲面的過程中,充分考慮原始數(shù)據(jù)點(diǎn)的分布特點(diǎn)和曲面構(gòu)造的精度要求,對(duì)于集中且均勻分布的原始數(shù)據(jù),采用插值方法,使曲面嚴(yán)格通過這些數(shù)據(jù)點(diǎn);對(duì)于分布離散不均勻的數(shù)據(jù),若精度要求很高則采用插值方法進(jìn)行加密,一般情況下采用逼近擬合方法,使曲面在給定精度下充分逼近原始數(shù)據(jù)。在水利水電工程壩區(qū),)地質(zhì)勘探數(shù)據(jù)集中且分布較均勻,其構(gòu)造的地質(zhì)曲面精度要求很高,這部分須采用插值方法;而在其相關(guān)的周邊地區(qū),對(duì)精度要求不高,采用逼近方法是可行的。
NURBS(Non-uniform rational B-splines,非均勻有理樣條曲線)技術(shù)是ISO(International Standard Organized,國際標(biāo)準(zhǔn)化組織)1991年頒布的STEP標(biāo)準(zhǔn)中自由型曲線曲面的唯一表示方法,它是在B樣條函數(shù)基礎(chǔ)上發(fā)展起來的,利用非均勻節(jié)點(diǎn)向量表達(dá)式來構(gòu)造有理B樣條,它對(duì)標(biāo)準(zhǔn)的解析圖形(如圓錐曲線、二次曲面、旋轉(zhuǎ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)的無規(guī)律性變化,選擇NURBS技術(shù)進(jìn)行水利水電工程三維地質(zhì)建模與分析,具有節(jié)省存儲(chǔ)空間、計(jì)算機(jī)處理簡便易行、數(shù)據(jù)庫管理方便、保證空間唯一性和幾何不變性等優(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基本上就可以滿足要求。
以圖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)格通過集合D1中的所有數(shù)據(jù)點(diǎn)。該方法使曲面在進(jìn)行插值時(shí)順序地通過一系列截面曲線,從而很好地生成復(fù)雜的自由曲面。在NURBS蒙皮過程中,首先要求各截面曲線的拓?fù)湫再|(zhì)應(yīng)該一致,即都是開曲線或都是閉曲線;其次在蒙皮插值過程中,插值的方式將影響最終的結(jié)果,根據(jù)工程地質(zhì)實(shí)際需要,一般選擇線性插值。圖3(b)為子集D1通過蒙皮插值得到的曲面S1。
2)對(duì)于子集D2和D3,則直接利用NURBS曲面技術(shù)采用反算法[15]進(jìn)行擬合。為了滿足曲面的邊界約束條件(鄰接曲面邊界和區(qū)域邊界約束),這里直接提取已構(gòu)建曲面S1與子集D2、D3相鄰的邊界線加入子集D2、D3,同時(shí)根據(jù)研究區(qū)域?qū)⒃紨?shù)據(jù)轉(zhuǎn)化為一系列u、v方向的曲線矢量,利用反算法求出各曲線矢量上的控制點(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ù)均來自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
在完成上述地質(zhì)曲面的構(gòu)造后,需要對(duì)其進(jìn)行檢查分析與調(diào)整,以滿足實(shí)際精度要求和建模需要。這主要從以下兩個(gè)方面進(jìn)行:
1)地質(zhì)結(jié)構(gòu)合理性和曲面幾何性檢查。即檢查所擬合的地質(zhì)結(jié)構(gòu)面整體變化趨勢是否合理、在幾何結(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)造曲面間的距離。
若σ小于控制誤差ε,則曲面滿足要求,否則需要對(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,這主要是由逼近擬合曲面引起的誤差,小于控制誤差,精度已足以滿足實(shí)際需要。
文章針對(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)格通過這些數(shù)據(jù)點(diǎn);對(duì)于周邊區(qū)域分布離散的數(shù)據(jù),采用逼近擬合方法,使曲面在給定精度下充分逼近原始數(shù)據(jù);并可通過地質(zhì)結(jié)構(gòu)合理性檢查、曲面幾何性和原始數(shù)據(jù)精度分析對(duì)構(gòu)造曲面進(jìn)行檢查分析和調(diào)整。實(shí)例研究表明,該方法所構(gòu)造的地質(zhì)曲面能滿足地質(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ī)場插值方法應(yīng)用研究[J].巖土力學(xué),2005,26(2):221-224.
[7]李曉軍,王長虹,朱合華.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].北京:中國水利水電出版社,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.