呂希奎，白嬌嬌，莊建杰，聶良濤，高崇文

(1.石家莊鐵道大學(xué) 交通運輸學(xué)院，河北 石家莊 050043；2.河北雄安軌道快線有限責(zé)任公司，河北 保定 071700；3.中國鐵路設(shè)計集團(tuán)有限公司，天津 300308)

山區(qū)普遍構(gòu)造活動強烈,工程地質(zhì)條件復(fù)雜,鐵路選線設(shè)計面臨著巨大難度。如何在地質(zhì)條件復(fù)雜的山區(qū)環(huán)境中科學(xué)合理的規(guī)劃鐵路線位,將鐵路的地質(zhì)災(zāi)害問題解決在成災(zāi)之前,已成為國內(nèi)外學(xué)者關(guān)注的焦點問題[1]。全面整合、整理和分析地質(zhì)信息數(shù)據(jù),有效利用地質(zhì)信息進(jìn)行選線,解決復(fù)雜地質(zhì)區(qū)域選線難題,研究人員嘗試對二維地質(zhì)信息進(jìn)行三維可視化建模。但復(fù)雜山區(qū)的地質(zhì)災(zāi)害和地質(zhì)構(gòu)造的復(fù)雜性、多樣性、地質(zhì)數(shù)據(jù)的稀疏性和地質(zhì)數(shù)據(jù)的海量特性,使得鐵路選線大范圍復(fù)雜地質(zhì)環(huán)境的三維建模一直未能真正實現(xiàn),成為難點問題。針對三維地質(zhì)建模,許多學(xué)者都進(jìn)行了研究。李健等[2]利用鉆孔與地質(zhì)剖面資料進(jìn)行數(shù)據(jù)轉(zhuǎn)換融合,通過插值算法加密地層矢量點,提出一種基于規(guī)則體元分裂的三維地質(zhì)建模方法。文獻(xiàn)[3-5]結(jié)合ArcGIS、SketchUp等三維建模軟件,實現(xiàn)一種半自動的復(fù)雜地質(zhì)體三維建模和可視化方法,并將其與三維地形環(huán)境中的不良地質(zhì)區(qū)域進(jìn)行疊加。陳兵等[6]基于BIM技術(shù)了建立了小范圍的橋址區(qū)三維地質(zhì)環(huán)境模型,用于更精準(zhǔn)地揭示橋址區(qū)的巖溶分布情況。文獻(xiàn)[7-8]基于鉆孔數(shù)據(jù),采用CATIA二次開發(fā)技術(shù),實現(xiàn)了完成了三維地層實體的建模。Hassen等[9]基于地表地質(zhì)和高程數(shù)據(jù),使用三維隱式方法構(gòu)建了法國諾曼底Pays d’Auge高原區(qū)域尺度三維地質(zhì)模型。文獻(xiàn)[10-12]應(yīng)用DEM和鉆孔數(shù)據(jù),結(jié)合GIS平臺和三維建模軟件,分別實現(xiàn)了杜卡拉地區(qū)、摩洛哥東北部Gareb Bouareg盆地和希臘克里特島Tymbakion盆地的三維地質(zhì)模型構(gòu)建。文獻(xiàn)[13-14]利用既有二維地質(zhì)成果,通過參數(shù)化,交互式快速建立了地質(zhì)三維模型。Olierook等[15]研究出地表地質(zhì)觀測與地球物理調(diào)查數(shù)據(jù)相融合可以建立區(qū)域三維地質(zhì)模型方法。

綜上研究成果,目前三維地質(zhì)模型的構(gòu)建主要依賴大量鉆孔資料基礎(chǔ)上進(jìn)行的三維地質(zhì)建模方法,建模的范圍也以小范圍區(qū)域為主。鐵路選線設(shè)計的分階段勘察特性,使得鐵路預(yù)可研和可研的前期規(guī)劃階段所獲取的地質(zhì)資料往往十分匱乏,更是缺少目前三維地質(zhì)建模模型時所需的海量勘探鉆孔數(shù)據(jù)[16]這使得前期選線規(guī)劃階段的三維地質(zhì)建模面臨著難題。在前期規(guī)劃設(shè)計階段,衛(wèi)星遙感是獲取地質(zhì)信息主要的應(yīng)用手段[17],因此,設(shè)計人員在預(yù)可研階段也主要通過二維柵格地質(zhì)遙感解譯圖進(jìn)行粗略的線路規(guī)劃方案設(shè)計,但在進(jìn)行地質(zhì)條件分析及方案比選時缺少直觀性,導(dǎo)致預(yù)可研階段設(shè)計的方案往往與最終設(shè)計方案存在較大的偏差。本文針對鐵路選線前期規(guī)劃設(shè)計階段三維地質(zhì)建模難題,基于空間數(shù)據(jù)三維可視化建模技術(shù),研究鐵路大范圍復(fù)雜地質(zhì)區(qū)域環(huán)境整體三維可視化方法,解決以下關(guān)鍵問題:①選線整體地質(zhì)區(qū)域環(huán)境海量三維地質(zhì)建模數(shù)據(jù)處理和矢量化地質(zhì)對象生成算法;②解決復(fù)雜山區(qū)大范圍所有地質(zhì)信息一次性整體三維建模和可視化關(guān)鍵問題;③實現(xiàn)在海量三維地質(zhì)區(qū)域環(huán)境中的選線應(yīng)用方法。實現(xiàn)充分展示整個鐵路區(qū)域的三維地質(zhì)環(huán)境,實現(xiàn)三維地質(zhì)區(qū)域環(huán)境和三維地理信息的統(tǒng)一集成和應(yīng)用,為預(yù)可研和可研階段復(fù)雜山區(qū)科學(xué)合理的路規(guī)鐵劃線位提供解決方案。

1 基于遙感地質(zhì)解譯圖的地質(zhì)對象空間數(shù)據(jù)生成方法

1.1 遙感解譯地質(zhì)對象的分類處理

遙感地質(zhì)解譯圖中存儲著多種地質(zhì)體、不良地質(zhì)和地質(zhì)現(xiàn)象等大量的地質(zhì)信息,如圖1所示。

圖1 地質(zhì)遙感解譯圖示例

這些地質(zhì)信息以柵格圖像格式在地質(zhì)解譯圖上進(jìn)行表達(dá)。解譯的地質(zhì)類型可在地質(zhì)解譯圖圖例中獲取,見圖2。

圖2 地質(zhì)解譯圖圖例

為更好地為三維數(shù)字化選線設(shè)計服務(wù),需要將柵格化的地質(zhì)信息提取并轉(zhuǎn)化為矢量化信息才能被計算機識別。對于各類地質(zhì)對象的邊界信息,利用遙感解譯軟件AOI功能得到各解譯地質(zhì)對象的二維邊界坐標(biāo)數(shù)據(jù)(x,y),為實現(xiàn)三維建模則需要進(jìn)行三維坐標(biāo)(x,y,z)的計算處理。

依據(jù)各類地質(zhì)對象在遙感解譯圖上的標(biāo)識形式,將圖1中的地質(zhì)類型劃分為點狀、線狀和面狀三類,詳細(xì)劃分見表1。

表1 地質(zhì)類型分類

1.2 點狀空間地質(zhì)對象生成

點狀范圍地質(zhì)對象可用三維點數(shù)據(jù)(x,y,z)格式表示,利用解譯的地質(zhì)對象數(shù)據(jù)二維邊界坐標(biāo)數(shù)據(jù),結(jié)合數(shù)字高程模型,以邊界點Pj(xj,yj)與內(nèi)插點P(x,y)之間距離dj平方的倒數(shù)作為權(quán)重ωj,基于Pj的高程Zj,采用加權(quán)平均法內(nèi)插計算邊界坐標(biāo)點高程Zp,從而完成三維坐標(biāo)點數(shù)據(jù)的生成,計算式為

(1)

利用此方法可確定古火山口、溫泉、熱泉等點狀地質(zhì)對象三維坐標(biāo)數(shù)據(jù),即可定位其空間位置。

1.3 線狀空間地質(zhì)對象生成

線狀地質(zhì)對象實體主要是指在幾何形態(tài)上呈線狀的地質(zhì)空間實體[18],其空間邊界由一系列有序的控制點連接而成,如斷層、斷裂等。拉格朗日插值曲線原理是從一系列控制點生成一系列點構(gòu)成的插值曲線,因此,為更真實地在三維空間中表達(dá)斷層等具有不同走向或形狀的線性地質(zhì)對象,對其控制點集進(jìn)行拉格朗日插值處理,定義如下:

線狀地質(zhì)對象三維控制點坐標(biāo)集Qi(xi,yi,zi),i=0,1,…,n,取其二維平面坐標(biāo)集構(gòu)造n個插值基函數(shù),利用插值基函數(shù)和已知節(jié)點坐標(biāo)求得插值多項式,再使用插值多項式獲得插值點坐標(biāo)。

假定二維平面坐標(biāo)集中的點滿足條件

(2)

為滿足條件式(2),構(gòu)建n次插值基函數(shù)為

i=0,1,2,…,n

(3)

則滿足條件式(2)的插值多項式為

(4)

通過求解線性方程組(4),解得Ln(xp),以此求得插值點的二維平面坐標(biāo),插值點Pi(i=0,1,2,…,n)的高程Zi通過式(1)求解獲得。最終連接空間三維型值點Pi(xi,yi,zi),即可實現(xiàn)線狀地質(zhì)對象在三維空間上的繪制,見圖3。其中,紅色點為地質(zhì)對象原有型值點,藍(lán)色點為拉格朗日插值點。繪制的線狀地質(zhì)對象,見圖4。

圖4 線狀地質(zhì)對象繪制效果

1.4 面狀空間地質(zhì)對象生成

滑坡、泥石流、巖土層界線等特殊巖土體都屬于面狀地質(zhì)實體,這些地質(zhì)對象在空間分布上呈封閉多邊形。面狀地質(zhì)實體邊界由空間型值點序列組成,采用分段三次Hermite插值法對面狀地質(zhì)對象邊界曲線進(jìn)行插值,可同時滿足邊界曲線的連續(xù)性和光滑性要求,能夠很好地描述面狀地質(zhì)對象,構(gòu)造過程如下所述。

面狀地質(zhì)對象邊界函數(shù)構(gòu)造方法。設(shè)有P1,P2,…,Pn共n個空間型值點,Pi(i=0,1,…,n)的型值點坐標(biāo)為(xi,yi),相鄰兩點之間通過三次Hermite曲線段連接,插值多項式H(x)可計算為

(5)

插值多項式求解。

Hi(x),x∈[xi-1,xi]滿足如下條件:

Hi(xi-1)=f(xi-1)=yi-1Hi(xi)=f(xi)=yi

Hi′(xi-1)=f′(xi-1)=y′i-1Hi′(xi)=f′(xi)=y′i

通過φi-1(x)、φi(x)、ψi-1(x)、ψi(x)四個三次多項式作為基函數(shù)構(gòu)造Hi(x),得

Hi(x)=yi-1φi-1(x)+yiφi(x)+y′i-1ψi-1(x)+y′iψi(x)

(6)

則由插值條件可以分別給出基函數(shù)滿足的條件為

令hi=xi-xi-1,解此方程組得

由此得到

(7)

這樣便求出滿足面狀地質(zhì)對象邊界的分段三次Hermite插值多項式。將面狀地質(zhì)對象的邊界輪廓等分為m個區(qū)間,對每個區(qū)間分別構(gòu)造三次Hermite插值多項式,利用方程(7)計算各插值點的坐標(biāo)并添加至原型值點集,然后按式(1)的加權(quán)平均法計算各型值點的高程值,最終在三維空間將該區(qū)間的點依次連接,繪制出面狀地質(zhì)對象,見圖5。

圖5 面狀地質(zhì)對象繪制效果

2 矢量化地質(zhì)對象三維模型創(chuàng)建方法

2.1 地質(zhì)對象矢量化模型創(chuàng)建

以插值處理后的地質(zhì)對象空間點集為數(shù)據(jù)源,基于CADC++接口,開發(fā)地質(zhì)遙感解譯圖矢量化地質(zhì)對象三維模型創(chuàng)建模塊,實現(xiàn)點、線和面狀矢量化地質(zhì)對象模型的自動創(chuàng)建。生成的矢量化地質(zhì)對象模型見圖6。

圖6 生成的矢量模式地質(zhì)對象模型

2.2 地質(zhì)對象矢量化顏色填充模式模型創(chuàng)建

為更充分展示各地質(zhì)對象模型,在矢量化基礎(chǔ)上實現(xiàn)顏色填充模式建模。

(8)

式中:m為地質(zhì)對象i的邊界數(shù)據(jù)點數(shù)量。

建立色系,其中線狀質(zhì)對象實測斷裂和推薦斷裂為紅色,解譯斷裂為藍(lán)色,標(biāo)注文字為黑色。

對面狀地質(zhì)對象建立RGB填充色系,顏色由淺到深,見表2。

表2 面狀地質(zhì)對象填充顏色色系

統(tǒng)計各類地質(zhì)對象所占范圍面積總和,按照面積由大到小順序排列,進(jìn)行填充顏色確定,所占比例越高,則其配色越淺,以達(dá)到更為清晰的地質(zhì)展現(xiàn)效果。

最終得到的矢量化地質(zhì)填充模式建模效果可直觀地展現(xiàn)各類地質(zhì)對象的分布情況,見圖7。圖7(a)為原始地質(zhì)遙感解譯柵格圖,7(b)為對應(yīng)的矢量化填充模式地質(zhì)建模效果。

圖7 填充模式的地質(zhì)建模效果

3 地質(zhì)信息三維空間定位算法

3.1 基于七參數(shù)的坐標(biāo)系轉(zhuǎn)換方法

為實現(xiàn)大范圍矢量化地質(zhì)模型在3DGIS中集成和展示,需要解決二者的坐標(biāo)系匹配問題,本文采用基于七參數(shù)的坐標(biāo)系轉(zhuǎn)換方法(以下簡稱七參數(shù)法)解決,其基于橢球間三維轉(zhuǎn)換,是由兩個坐標(biāo)系之間的平移值Δx0、y2平移值Δy0、z2平移值Δz0、x軸旋轉(zhuǎn)角度ωx、y軸旋轉(zhuǎn)角度ωy、z軸旋轉(zhuǎn)角度ωz,縮放尺度比m組成,數(shù)學(xué)模型為

(9)

七參數(shù)法坐標(biāo)轉(zhuǎn)換見圖8。由七參數(shù)法數(shù)學(xué)模型可知,要求解七參數(shù)最少需要知道原坐標(biāo)系與目標(biāo)坐標(biāo)系中位置相對應(yīng)的三組控制點?？蓪⒛繕?biāo)坐標(biāo)系中的三維地形坐標(biāo)賦予地質(zhì)矢量模型中對應(yīng)的控制點,計算出坐標(biāo)轉(zhuǎn)換七參數(shù),進(jìn)而計算得到地質(zhì)矢量模型中各點的三維坐標(biāo),實現(xiàn)地質(zhì)矢量模型在3DGIS環(huán)境的集成。

圖8 七參數(shù)法示意

3.2 控制點選取策略研究

七參數(shù)法中的控制點選取精確與否是影響匹配精度的主要因素,其誤差主要體現(xiàn)在空間偏移距離與俯瞰圖像匹配度兩方面。本文擬定控制點的選取范圍參數(shù)(控制點三點所圍成的三角形區(qū)域面積與正射影像圖(DOM)影像整體面積的比值r)與圖像匹配精度參數(shù)(兩個標(biāo)記區(qū)域的空間距離誤差L、原三維地理環(huán)境的標(biāo)記區(qū)域圖像(基準(zhǔn)圖)與地質(zhì)矢量模型上相對應(yīng)的標(biāo)記區(qū)域圖像(匹配圖)各灰度狀態(tài)比值的均值A(chǔ)vg),研究通過對不同范圍內(nèi)的控制點組的結(jié)果進(jìn)行對比,得到兩類參數(shù)之間的相關(guān)關(guān)系,從而確定能夠得到高匹配度結(jié)果的控制點選取策略。

3.2.1 基于灰度直方圖的圖像匹配度分析方法

本文以冰磧湖為例,將部分冰磧湖區(qū)域設(shè)定為標(biāo)記區(qū)域,其中紅色邊界線為冰磧湖邊界輪廓,黃色邊界為矢量化模型邊界。通過調(diào)整地質(zhì)矢量模型的透明度,將地形覆蓋前后的標(biāo)記區(qū)域范圍進(jìn)行對比,利用圖像處理與分離技術(shù),得到基準(zhǔn)圖與匹配圖,圖像分離過程見圖9。

圖9 標(biāo)記區(qū)域圖像分離

將圖像分離完畢后,利用計算機獲取基準(zhǔn)圖與匹配圖的灰度直方圖,將基準(zhǔn)圖H1與匹配圖H2的直方圖相似度數(shù)值d作為評價匹配精度的指標(biāo)之一。采用OpenCV計算直方圖相似度的數(shù)值方法如下:

采用相關(guān)對比標(biāo)準(zhǔn)(CV_COMP_CORREL)計算d(H1,H2),d(H1,H2)越接近1,表示兩個圖像匹配度越高,計算式為

計算流程如下:

1)裝載基準(zhǔn)圖像和匹配圖像進(jìn)行對比。

2)轉(zhuǎn)換圖像為HSV圖像格式。

3)計算所有圖像的H-S直方圖并歸一化。

4)對比基準(zhǔn)圖像與匹配圖像的直方圖。

5)計算直方圖相似度數(shù)值。

若完全匹配,則兩邊界完全重合,直方圖相似度數(shù)值為1,圖10 中則具有一定的偏差。標(biāo)記區(qū)域空間距離誤差L越小,則重合度越高。

圖10 兩個多邊形區(qū)域

3.2.2 標(biāo)記區(qū)域空間距離誤差計算方法

Hausdorff距離是一種基于輪廓的匹配方法,用來衡量兩個點集間的距離[19]。設(shè)A={a1,a2,…,an},ai(xi,yi)∈R2,i=1,2,…,n;B={b1,b2,…,bm},bj(xj,yj)∈R2,j=1,2,…,m為兩個多邊形區(qū)域有限集,輪廓邊緣由連續(xù)的點表示,見圖10。

A,B之間的雙向Hausdorff距離定義為

H(A,B)=max{h(A,B),h(B,A)}

(11)

式中:‖·‖為點集A和B的歐氏距離。

Hausdorff距離H(A,B)取h(A,B)和h(B,A)的最大值,H(A,B)越小,A,B點集之間匹配程度越高[20]。因此,Hausdorff距離匹配方法合適計算標(biāo)記區(qū)域與疊加區(qū)域的距離L值,將L取值為雙向Hausdorff距離。

取D(a,b)為a與b之間的歐幾里得距離,采用窮舉算法計算集合A和集合B兩個點集距離,算法流程見圖11。

3.2.3 控制點的選取策略

通過五組控制點組得到的控制點選取范圍參數(shù)與匹配度參數(shù)見表3。由表3可以看出,控制點三點所圍成的三角形區(qū)域面積與DOM影像整體面積的比值r,影響著圖像直方圖相似度系數(shù)d與標(biāo)記區(qū)域空間距離誤差L:r越大,L越小,d 越大,圖像匹配度越高。

表3 控制點選取范圍參數(shù)與匹配度參數(shù)

為進(jìn)一步研究兩類參數(shù)的關(guān)系,采用Matlab對表3中的控制點數(shù)據(jù)擬合,得到了在95%的置信區(qū)間下,控制點選取范圍參數(shù)與圖像匹配參數(shù)的函數(shù)關(guān)系表達(dá)式與關(guān)系圖像,見圖12和圖13。根據(jù)圖像匹配度要求,可反算出控制點的選取范圍。

圖12 控制點區(qū)域占比與圖像灰度比均值擬合

圖13 控制點區(qū)域占比與空間距離誤差擬合

r與d的函數(shù)關(guān)系為

d=0.105 2×r1.561+73.09

r與L的函數(shù)關(guān)系為

L=6 922×r-0.014 2-6 593

從圖12和圖13中可以得出,面積比值r越大,圖像直方圖相似度值d越大,標(biāo)記區(qū)域空間距離誤差L越小,二者分別呈指數(shù)增長和指數(shù)遞減關(guān)系。

3.3 地質(zhì)信息三維呈現(xiàn)方法

本文選用WGS84坐標(biāo)系(代碼為LL84),將DEM與DOM影像均以光柵文件在InfraworksBIM平臺中重建三維場景,生成項目所在地真實的三維地形環(huán)境,局部地形效果見圖14。

圖14 局部地形效果

采用最優(yōu)控制點的選取策略進(jìn)行控制點選取和空間定位,采用地形覆蓋形式,將整體矢量化地質(zhì)模型在Infraworks中集成,實現(xiàn)大范圍復(fù)雜地質(zhì)信息一次性整體三維建模和可視化,圖15為示例控制點組,圖16為實現(xiàn)方法流程。

圖15 控制點組選取

圖16 整體三維建模和可視化方法實現(xiàn)流程

整體一次性的大范圍地質(zhì)區(qū)域環(huán)境三維建模方法,能夠快速建立選線區(qū)域完整的復(fù)雜地質(zhì)環(huán)境,能夠更加清晰直觀地對地質(zhì)情況進(jìn)行展示。在建立的三維地質(zhì)區(qū)域環(huán)境中,可導(dǎo)入已有設(shè)計方案。為更直觀地展示設(shè)計方案三維效果,將Revit參數(shù)化生成的線路BIM模型沿線路中心線進(jìn)行空間定位,得到在大范圍三維地質(zhì)區(qū)域環(huán)境中設(shè)計方案BIM模型集成展示。同時可更清晰直觀地了解和查詢線路方案周邊區(qū)域的地質(zhì)信息。

4 結(jié)論

1)以衛(wèi)星遙感地質(zhì)解譯影像為主要數(shù)據(jù)源,實現(xiàn)了基于遙感地質(zhì)解譯圖的矢量化信息提取與處理方法。采用拉格朗日插值和三次Hermite曲線插值算法,實現(xiàn)了數(shù)據(jù)驅(qū)動的點狀、線狀和面狀矢量化地質(zhì)對象三維模型創(chuàng)建方法。

2)以插值處理后的地質(zhì)對象空間點集為數(shù)據(jù)源,對各類地質(zhì)對象所占范圍面積總和比進(jìn)行配色選取,實現(xiàn)了顏色填充模式的矢量化地質(zhì)環(huán)境建模效果。

3)以面積比值、雙向Hausdorff距離的空間距離誤差和灰度狀態(tài)比值的均值,建立了基于灰度直方圖的圖像匹配度分析方法和最優(yōu)控制點選取策略,實現(xiàn)了地質(zhì)信息三維空間定位算法。

4)基于七參數(shù)法以地形覆蓋形式實現(xiàn)了完整矢量化地質(zhì)模型在InfraworksBIM平臺的空間定位、集成和三維展示,實現(xiàn)了大范圍復(fù)雜地質(zhì)區(qū)域環(huán)境的三維可視化建模。

5)為進(jìn)一步深入應(yīng)用,將設(shè)計方案和鐵路線路BIM設(shè)計結(jié)果在InfraworksBIM平臺集成,實現(xiàn)了在三維環(huán)境中無損展示三維線路方案BIM模型設(shè)計效果。研究成果為預(yù)可研階段的大范圍三維地質(zhì)建模和選線應(yīng)用提供了可行的解決方案。