李勇 段曉峰 王宏遠(yuǎn)

蘭州交通大學(xué) 土木工程學(xué)院，蘭州 730070

鐵路線路地形具有長帶狀特征，現(xiàn)有的BIM 和GIS（Geographic Information Science）軟件使用過程中地形模型數(shù)據(jù)量過大且冗余過多，導(dǎo)致加載緩慢或長時(shí)間無法響應(yīng)。若大范圍刪減，又可能因刪減不當(dāng)導(dǎo)致數(shù)據(jù)點(diǎn)無法較好地表征實(shí)際地形特征；同時(shí)，構(gòu)網(wǎng)方式單一化，無法滿足工程多式樣的地形顯示。

在地形建模和輕量化研究方面，王倩［1］研究了三角形網(wǎng)格地形建模的基本方法和理論，包括規(guī)則數(shù)據(jù)及不規(guī)則數(shù)據(jù)三角網(wǎng)的生成算法。文雄志［2］運(yùn)用層次細(xì)節(jié)技術(shù)和基于網(wǎng)格的模型簡化技術(shù)對(duì)三維地形建模進(jìn)行簡化，使其能夠更好地實(shí)現(xiàn)建模簡化技術(shù)的應(yīng)用。陸風(fēng)等［3］對(duì)基于地形線的Memoryless 改進(jìn)算法的簡化效果進(jìn)行了驗(yàn)證，顯示水下地形模型簡化后保持了原有特征和細(xì)節(jié)，但模型會(huì)產(chǎn)生跳躍現(xiàn)象。宋偉華等［4］通過Delaunay 逐點(diǎn)插入算法構(gòu)建TIN 模型，該算法能用更少的數(shù)據(jù)表達(dá)復(fù)雜地表形態(tài)，但時(shí)間復(fù)雜度差、運(yùn)行速度慢。分形理論的出現(xiàn)和發(fā)展，為描述自然界的復(fù)雜事物提供了解決思路。李鑫［5］利用數(shù)字圖像處理技術(shù)實(shí)現(xiàn)了獲取分形維數(shù)的流程，具有較好的測量穩(wěn)定性，無需占用太多計(jì)算機(jī)資源即可獲得較高精度。張欣悅等［6］提出利用有理分形插值進(jìn)行分形曲線建模，驗(yàn)證了可變參數(shù)的有理分形插值曲線在處理實(shí)際問題中具有較好的應(yīng)用價(jià)值。彭曉蕾等［7］對(duì)蘊(yùn)含地貌特征的黃土高原地貌數(shù)據(jù)進(jìn)行分形設(shè)計(jì)，將自然景觀特征與公式化的數(shù)學(xué)思維相結(jié)合，為自然景觀的美學(xué)研究提供了一個(gè)新的方向。

鑒于分形算法由極少部分地形控制點(diǎn)就可表達(dá)出原地形的特征，還能在保持地形真實(shí)感的狀態(tài)下根據(jù)實(shí)際需求動(dòng)態(tài)控制地形的細(xì)化程度，在地形構(gòu)網(wǎng)上具有優(yōu)越性，因此本文嘗試基于分形算法解決BIM 地形模型輕量化表達(dá)問題。

1 分形理論概述

1.1 分形的概念及維數(shù)

若集合F在歐氏空間中的Hausdorff 維數(shù)DH恒大于其拓?fù)渚S數(shù)Dr，即DH>Dr，則該集合是分形集，簡稱分形［8-9］。分形維數(shù)是一種非整數(shù)維數(shù)，用來解決整數(shù)維在尺度上出現(xiàn)過大或過小而無法準(zhǔn)確評(píng)價(jià)的問題。如果將分形拓展到表示自然界中的行為或現(xiàn)象，那么分形維數(shù)描述的則是零碎的局部特征構(gòu)成整體系統(tǒng)行為的相關(guān)性。

1.2 地形分形維數(shù)的計(jì)算

地形分形維數(shù)計(jì)算常用的方法為表面積-體積法和方差圖法。

1.2.1 表面積-體積法

用Hausdorff 面積SH表示地形面積，V表示地形曲面的體積。由量綱分析可得

式中：D為地形曲面的分形維數(shù)。

若用S表示地形曲面的歐式面積，則有

式中：k0為系數(shù)，又稱形狀因子；l為尺度值。

則可推導(dǎo)出直接計(jì)算分形維數(shù)的公式為

1.2.2 方差圖法

方差圖法是建立在圖像（分形數(shù)據(jù)）的高斯統(tǒng)計(jì)模型上的。給定一個(gè)分形維數(shù)D和方差δ，用分形布朗運(yùn)動(dòng)（Fractional Brownian Motion，F(xiàn)BM）模型可模擬生成一幅圖像。反之，如果將圖像或地形表面看成是分形布朗曲面，就可以得到其FBM模型的分形維數(shù)。

2 基于分形理論的數(shù)字地形模型生成算法

用于地形模擬的分形算法主要為基于迭代函數(shù)系統(tǒng)（Iterated Function System，IFS）方法和基于分形布朗運(yùn)動(dòng)（FBM）方法［9-10］。IFS 方法是對(duì)已有數(shù)據(jù)集進(jìn)行仿射變換，構(gòu)造一類插值函數(shù)，進(jìn)而通過迭代過程，得到需要的形狀。FBM 方法通過控制表征值H、δ等參數(shù)，加入隨機(jī)變量然后構(gòu)造出所需的分形圖形［11］。本文采用基于FBM的方法對(duì)地形進(jìn)行重構(gòu)。

FBM 是一種隨機(jī)過程，其增量按正態(tài)分布［12-13］。將FBM 二維曲面定義為在某概率空間上的隨機(jī)過程，且滿足：X（0，0）=0，X（x，y）是（x，y）的連續(xù)函數(shù)。

對(duì)任意變量（x，y），(h，k) ∈R2，其二維增量X(x+h，y+k) -X(x，y)服從期望為0，方差為(h2+k2)2H的正態(tài)分布，其概率滿足

基于FBM 的插值方法主要有泊松階躍法、隨機(jī)中點(diǎn)位移法（Random Midpoint Displacement，RMD）、傅立葉濾波法、逐次隨機(jī)增加法等［13］。本文采用基于FBM的改進(jìn)隨機(jī)中點(diǎn)位移法，通過對(duì)OpenRail 軟件二次開發(fā)以改進(jìn)優(yōu)化原有地形生成時(shí)的臃腫情況。

中點(diǎn)位移法是標(biāo)準(zhǔn)的分形幾何法，其在細(xì)分過程中，通過對(duì)兩個(gè)或多個(gè)頂點(diǎn)之間插值進(jìn)行地形建模，具有高速度以及為已有形狀增加細(xì)節(jié)的能力，但在實(shí)際生成地形過程中地形曲面會(huì)出現(xiàn)褶皺現(xiàn)象。為此采用一種加入與分形特征表征值H、σ相關(guān)的補(bǔ)償項(xiàng)的改進(jìn)隨機(jī)中點(diǎn)位移法。該方法既能滿足用分形布朗運(yùn)動(dòng)特征的分形表面來表達(dá)自然地形，同時(shí)又能很好解決此問題［14-16］。

設(shè)01，02，…，0n為初始柵格數(shù)據(jù)點(diǎn)，格網(wǎng)間距為d0。每一級(jí)細(xì)分包括兩步內(nèi)插過程。

1）第一步：diamond步

由最鄰近的四個(gè)格網(wǎng)點(diǎn)高程h0i（i=1，2，3，4），內(nèi)插其中心點(diǎn)h1i的高程。高程值等于四個(gè)格網(wǎng)點(diǎn)高程的平均值加上一個(gè)隨機(jī)位移Δ ～N(0，V2?)。h1i計(jì)算式為

此時(shí)，原始的正方形格網(wǎng)旋轉(zhuǎn)45°成為菱形格網(wǎng)，格網(wǎng)間距也由d0變?yōu)閐1，即

2）第二步：Square步

由第一步得到的已知菱形四個(gè)角點(diǎn)高程內(nèi)插其中心點(diǎn)高程，即

經(jīng)過第二步，格網(wǎng)間距又變?yōu)閐2，即

如此，隨機(jī)中點(diǎn)位移法的一個(gè)完整步驟就完成了，如圖1所示，接著按照所需迭代次數(shù)反復(fù)進(jìn)行迭代計(jì)算就可得到所需的分形曲面。但存在隨機(jī)變量的取值問題。通常，對(duì)Δ 的選擇要使待插點(diǎn)的高程增量滿足FBM 所需的冪指數(shù)規(guī)律。令格網(wǎng)間距的變化尺度為γn（n為迭代水平），那么隨機(jī)位移的方差為

圖1 隨機(jī)中點(diǎn)位移法計(jì)算示意

由d1和d2可知γ=2-n/2，因此演變出常用的隨機(jī)位移公式，即

式中：sc為位移大小調(diào)整因子；G為高斯隨機(jī)分量，G～N(0，1)。

式（10）一般用于虛擬分形地形的生成，不適用于自然世界中實(shí)測數(shù)據(jù)的插值。為使待插點(diǎn)的分形特征符合實(shí)際數(shù)據(jù)的分形特征，Yokaya 提出了式（11）的隨機(jī)位移增量。

本文采用文獻(xiàn)［9］提出的改進(jìn)Yokaya 隨機(jī)增量，即

3 基于分形算法的地形模型重構(gòu)技術(shù)

為了直接將分形地形與實(shí)際線路BIM 設(shè)計(jì)相結(jié)合，選用OpenRail Designer進(jìn)行二次開發(fā)。

3.1 OpenRail二次開發(fā)環(huán)境搭建

本文所使用的軟件及配置為：①操作系統(tǒng)為Windows10，內(nèi)存16G，顯卡GTX1650，顯存8G；②編譯平臺(tái)為Visul Studio2015；③開發(fā)環(huán)境為.NET4.6.2 及以上版本的框架，與OpenRail版本匹配的SDK。

3.2 基于分形算法重構(gòu)地形模型

3.2.1 原始數(shù)據(jù)重采樣

本次地形數(shù)據(jù)采用無人機(jī)航測數(shù)據(jù)，在原有數(shù)據(jù)的基礎(chǔ)上進(jìn)行重采樣，以滿足算法需要。

由分形布朗運(yùn)動(dòng)相關(guān)知識(shí)可知，要計(jì)算分形維數(shù)D，需要不斷改變其尺度。因此本文在提取分形特征值時(shí)，考慮實(shí)際地形數(shù)據(jù)特征，將采樣間隔設(shè)置為20、40、50、60 m，且為判斷不同軸采樣對(duì)分形特征值的影響，分別沿?cái)?shù)據(jù)的x軸、y軸和對(duì)角線采樣，獲得不同的采樣結(jié)果，見圖2、圖3。

圖2 高程點(diǎn)重采樣

圖3 重采樣?xùn)鸥駡D

3.2.2 分形特征提取

假設(shè)f（x）為一個(gè)隨機(jī)分形布朗運(yùn)動(dòng)，則滿足式（13）的概率分布函數(shù)［14-15］。

式中：f（x）為采樣值；Δx為采樣間隔；H為頻譜指數(shù)；F（y）為累積概率分布函數(shù)，服從標(biāo)準(zhǔn)正態(tài)分布（0，δ2）。而分形維數(shù)D=N+1-H，N為拓?fù)渚S數(shù)。對(duì)地形曲面進(jìn)行分形模擬時(shí)，H表示地形復(fù)雜情況的一種抽象和概括，H越大分形布朗運(yùn)動(dòng)變化趨于平緩，H越小變化趨于精細(xì)；δ反映地形曲面的起伏特征。

式（13）又可改寫為

最小二乘法的直線回歸函數(shù)為

由此求得系數(shù)H、C，便可推得D和δ。

為比較采樣方向?qū)Φ匦翁卣鞯挠绊?，將整體數(shù)據(jù)劃分為四個(gè)大小一致的區(qū)塊依次計(jì)算各分塊分形特征值，用來比較區(qū)塊分形特征值與整體分形特征值的差別，見圖4、圖5。

圖4 整體數(shù)據(jù)分形特征擬合

圖5 各區(qū)塊分形特征擬合

由圖4（a）可知，沿x軸采樣數(shù)據(jù)整體只呈現(xiàn)一個(gè)無標(biāo)度區(qū)，采樣點(diǎn)數(shù)據(jù)分布趨勢(shì)基本一致，擬合可得該地形沿x軸的分形特征值為H=0.432 8，δ=3.845。

由圖4（b）可知，沿y軸采樣數(shù)據(jù)整體只呈現(xiàn)一個(gè)無標(biāo)度區(qū)，無錯(cuò)誤點(diǎn)，擬合可得該地形沿y軸的分形特征值為H=0.476 5，δ=1.906。

由圖4（c）可知，沿對(duì)角線采樣數(shù)據(jù)整體只呈現(xiàn)一個(gè)無標(biāo)度區(qū)，無錯(cuò)誤點(diǎn)，擬合可得該地形沿對(duì)角線的分形特征值為H=0.659 1，δ=0.845。

由圖5（a）可知，對(duì)區(qū)塊1 進(jìn)行采樣獲得的數(shù)據(jù)以lg Δr=8 為界限呈現(xiàn)兩個(gè)無標(biāo)度區(qū)，在lg Δr=7.8 后出現(xiàn)另一個(gè)無標(biāo)度區(qū)，與區(qū)塊1 多數(shù)數(shù)據(jù)呈現(xiàn)的特征不符，因此需要剔除后面的采樣點(diǎn)，最終擬合可得該區(qū)塊地形分形特征值為H=0.500 7，δ=1.998。

由圖5（b）可知，對(duì)區(qū)塊2 進(jìn)行采樣獲得的數(shù)據(jù)整體呈現(xiàn)一個(gè)無標(biāo)度區(qū)，最終擬合可得該區(qū)塊地形分形特征值為H=0.430 3，δ=2.228。

由圖5（c）可知，對(duì)區(qū)塊3 進(jìn)行采樣獲得的數(shù)據(jù)以lg Δr=7.8 為界限呈現(xiàn)兩個(gè)無標(biāo)度區(qū)，在lg Δr=7.8后出現(xiàn)另一個(gè)無標(biāo)度區(qū)，與區(qū)塊3 多數(shù)數(shù)據(jù)呈現(xiàn)的特征不符，因此擬合時(shí)只對(duì)第一個(gè)無標(biāo)度區(qū)進(jìn)行擬合，最終擬合可得該區(qū)塊地形分形特征值為H=0.529 1，δ=1.935。

由圖5（d）可知，對(duì)區(qū)塊4 進(jìn)行采樣獲得的數(shù)據(jù)呈現(xiàn)兩個(gè)無標(biāo)度區(qū)，在lg Δr=8.2 后出現(xiàn)另一個(gè)無標(biāo)度區(qū)，與區(qū)塊4多數(shù)數(shù)據(jù)呈現(xiàn)的特征不符，因此擬合時(shí)只對(duì)第一個(gè)無標(biāo)度區(qū)進(jìn)行擬合，最終擬合可得該區(qū)塊地形分形特征值為H=0.543 6，δ=2.202。

對(duì)所有的數(shù)據(jù)結(jié)果進(jìn)行整理，結(jié)果見表1?？芍?，無人機(jī)獲取的試驗(yàn)區(qū)域地形數(shù)據(jù)，無論采用三個(gè)方向?qū)φw數(shù)據(jù)采樣，或者是對(duì)數(shù)據(jù)劃分，提取出的分形特征值基本一致，分形維數(shù)D始終保持在2.5 左右。因此在進(jìn)行分形插值構(gòu)網(wǎng)時(shí)，可以采用全局使用一個(gè)分形特征值進(jìn)行插值算法的編譯。

表1 地形分形特征值統(tǒng)計(jì)

3.2.3 曲面生成

根據(jù)分形布朗運(yùn)動(dòng)原理，基于地形數(shù)據(jù)、重采樣數(shù)據(jù)及擬合出的分形特征值，利用改進(jìn)隨機(jī)中點(diǎn)位移法，通過編制分形程序并加載到BIM 軟件OpenRail Designer，可實(shí)現(xiàn)不同特征值的地形分形模擬。

對(duì)同一地形（圖6），采用區(qū)塊分形特征值（H=0.501，δ=2.091）和整體分形特征值（H=0.523，δ=2.199）進(jìn)行插值構(gòu)網(wǎng)，迭代一次得到圖7（a）和圖7（b），視覺上兩個(gè)地形的表達(dá)差異極小。

圖6 原始地形

當(dāng)使用整體分形特征值迭代兩次以后，得到的地形［圖7（c）］顯示效果優(yōu)秀，已經(jīng)完全能夠表達(dá)原地形的起伏特征，從表2的平均坡率也能驗(yàn)證這一點(diǎn)。

圖7 三種工況的分形地形

表2 原始地形與分形地形指標(biāo)對(duì)比

由表2 可知：第一次迭代點(diǎn)數(shù)為8 320，僅為原始地形點(diǎn)數(shù)的0.831%，網(wǎng)格數(shù)也僅為原始網(wǎng)格數(shù)的0.821%；經(jīng)第二次迭代，點(diǎn)數(shù)與網(wǎng)格數(shù)分別是原始數(shù)據(jù)的2.270%和2.286%。結(jié)合對(duì)比運(yùn)算時(shí)間，采用分形算法對(duì)地形進(jìn)行輕量化是切實(shí)可靠的，可用極少地形控制點(diǎn)完美實(shí)現(xiàn)地形特征表達(dá)，在節(jié)約內(nèi)存、提高地形構(gòu)網(wǎng)效率方面也有較大優(yōu)勢(shì)。

此方法能更好地解決鐵路BIM 與GIS設(shè)計(jì)過程中長帶狀地形數(shù)據(jù)量龐大、加載緩慢的問題。

4 結(jié)論

1）通過比較區(qū)塊與整體分形特征值，無論是采用三個(gè)方向?qū)φw數(shù)據(jù)采樣，或是對(duì)數(shù)據(jù)劃分，對(duì)每個(gè)區(qū)塊采樣進(jìn)行分形特征值的擬合，提取出的分形特征值基本一致，分形維數(shù)D始終保持在2.5左右。因此，在進(jìn)行分形插值構(gòu)網(wǎng)時(shí)，可以采用全局使用一個(gè)分形特征值進(jìn)行插值算法的編譯。

2）基于分形算法的地形模型重構(gòu)可以由極少部分地形控制點(diǎn)較好地表達(dá)原地形的特征，并能根據(jù)實(shí)際需求動(dòng)態(tài)控制地形的細(xì)化程度，解決軟件構(gòu)網(wǎng)方式單一化問題，滿足工程實(shí)際中多式樣的地形顯示需求，解決了地形輕量化表達(dá)問題。