基于距離場映射的地形可視域分析方法

王晨昊,湯曉安,馬伯寧,冷志光

1.國防科學(xué)技術(shù)大學(xué),湖南長沙410073;2.中國人民解放軍91635部隊,北京102249

1 引 言

地形可視性也稱為地形通視性(visibility),是指從一個或多個位置所能看到的地形范圍或與其他地形點之間的可見度?？梢曅苑治鍪堑匦畏治龅闹匾M成部分,很多與地形有關(guān)的問題都涉及地形可視性的計算問題。例如,觀察站、雷達位置、廣播電視或電話發(fā)射塔的位置、路徑規(guī)劃、航海導(dǎo)航、軍事上的陣地布設(shè)、道路和建筑物的景觀設(shè)計、日照分析等[1]。

點點通視(inter-visibility)分析和可視域(viewshed)分析是地形可視性分析的兩個基本內(nèi)容[2]。而可視域的計算從本質(zhì)上講,是兩點之間的通視計算在面域上的實現(xiàn),也就是說,給定視點的可視域就是通過視點的無數(shù)視線上的可見部分的集合[1]。目前已有的可視域分析算法全部是基于CPU實現(xiàn)的,算法復(fù)雜度較高,難以應(yīng)付大數(shù)據(jù)量地形的在線分析;此外,這些算法幾乎都是針對固定格式地形數(shù)據(jù)提出的,通用性較差。

針對以上問題,本文提出一種基于距離場映射的可視域分析方法。這種方法面向單觀察點地形可視域計算問題,可以在三維地形渲染輸出的同時獲取地形的可視范圍。與現(xiàn)有方法相比,這種方法具有以下特點:

(1)算法通過 GPU實現(xiàn),計算速度更快。

(2)可以通過調(diào)整距離圖分辨率來控制分析結(jié)果的精度,在處理速度和結(jié)果精度的平衡上具有較高靈活度。

(3)能同時適用于規(guī)則格網(wǎng)(RSG)和非規(guī)則格網(wǎng)(TIN)地形模型,也能用于地面物體的可視性分析。

2 相關(guān)研究

現(xiàn)有的地形可視域計算方法大部分是針對規(guī)則格網(wǎng)DEM提出的[3]。文獻[4]提出一種基于視線(line of sight,LOS)的方法,跟蹤從視點 O到任何其他點P的視線,并沿著O與P之間的射線移動,在到達 P之前檢測到視線與地形相交或到達P,就終止移動?；镜腖OS算法存在冗余計算,為此研究者對LOS方法做了許多改進,以期提高計算效率。相關(guān)成果有文獻[5]提出的基于并行計算系統(tǒng)的可視性分析模型、文獻[6]利用相鄰視點之間的內(nèi)在一致性提出的改進并行算法、文獻[7]提出的平面掃描算法、文獻[8]提出的利用輔助格網(wǎng)的改進算法、文獻[9]提出的增量地形可視性算法以及文獻[10]提出的基于斜率比較的可視性算法等。這些改進的LOS算法雖然一定程度上提高了執(zhí)行速度,但由于必須內(nèi)插計算與視線相交的DEM網(wǎng)格單元,計算量仍然比較大。

不同于LOS方法,文獻[11]提出一種參考面法,通過視點和目標(biāo)點之間的空間關(guān)系所形成的參考面來判斷觀測點與所有目標(biāo)點是否可視。雖然這種方法比任何基于視線的算法都快,但仍然存在實現(xiàn)復(fù)雜度高、通用性差的問題。文獻[12]則在分析可視性問題共性的基礎(chǔ)上構(gòu)建了一種DEM地形可視性統(tǒng)一分析模型,這個模型使得具體的可視性問題具有模式化的解決方法,但并未從速度、精度等方面對可視性算法進行改進。

面向TIN的可視域計算方法目前主要有前后順序法、分而治之法、掃描線法、光線投射法等[13]。對比面向規(guī)則格網(wǎng)DEM的可視性分析方法和面向TIN的方法,各種算法往往只針對一種形式的地形模型,而不能通用。另外,對于包含建筑物、植被等地物的地形區(qū)域,現(xiàn)有分析方法的適用性也受到挑戰(zhàn)。

本文方法受到文獻[14]提出的陰影圖方法的啟發(fā)。陰影圖方法用于真實感圖形繪制中物體陰影的實時生成。由于研究對象以及應(yīng)用側(cè)重點不同,本文方法與陰影圖方法的不同之處在于:①提出一種全新的地形可視域分析方法,根據(jù)應(yīng)用中是否考慮目標(biāo)高度,給出了兩種不同的可視性判斷算法;②針對地形距離圖的特點,提出一種偏移采樣法來提高可視性計算的精度。

3 本文方法

本文通過距離場映射的方法進行地形的可視性分析。主要包含兩個步驟:第一步,在觀察點以透視投影渲染地形,生成距離圖;第二步,在頂視圖以平行投影渲染地形,將距離圖映射到地形表面片元,并進行片元可視性的判斷。算法的執(zhí)行包含對地形的兩次渲染,流程如圖1所示。

圖1 本文方法的流程圖Fig.1 Flow chart of our method

3.1 生成距離圖

三維地形一般以表面模型的形式進行渲染,經(jīng)過深度測試(depth test),顯示設(shè)備可以獲得每條投影射線上與視點最近的片元,通常情況下將其顏色值寫入顏色緩存作為一個像素的顯示輸出,如圖2(a)所示。如果不將片元顏色作為輸出,而是將片元到觀察點的距離值輸出到顏色緩存,那么就可以獲得一個以觀察點為基準(zhǔn)點的地表距離場在屏幕上的二維投影,稱為距離圖,記為fM。距離圖記錄了從O出發(fā)的視線與地形表面的第一個交點到O的距離值,如圖2(b)所示。

圖2 地形渲染輸出Fig.2 The rendering outputs

當(dāng)關(guān)注的地形相對于視點的角度較大,甚至是全向時,僅僅利用一個二維紋理將難以滿足要求。這種情況下利用立方體貼圖[15](cube map)來存儲全向的距離場信息,如圖3所示。立方體貼圖包含六個面,每個面覆蓋90°的視場角。在支持Direct3D 10的顯卡上,利用幾何著色器(geometry shader)只需通過對場景的一次渲染就可以獲得立方體貼圖。立方體貼圖通過一個三維矢量來訪問。

圖3 立方體距離圖Fig.3 Distance cube map

3.2 計算可視域

通過在觀察點O處渲染地形,獲得了距離圖(立方體距離圖)fM。接下來通過對地形的第二次渲染來獲取地形可視域分析結(jié)果。為了能覆蓋關(guān)注的地形區(qū)域,通過自頂向下的平行投影將地形渲染到目標(biāo)緩沖區(qū)。經(jīng)過柵格化處理之后,地形的幾何數(shù)據(jù)被轉(zhuǎn)換為片元,每個片元都對應(yīng)于目標(biāo)緩存的一個像素。設(shè)世界坐標(biāo)系下觀察點O到片元點 P的矢量為vP,片元法矢量為 nP。在像素著色器中通過以下算法進行片元的可視性判斷。

算法1:

第一步,計算O到 P的歐氏距離d、片元的距離場映射 dM以及片元方向參數(shù)λ

第二步,判斷片元的可視性。如果d=dM,且λ<0,那么相對于觀察點,片元可見;否則,片元不可見,原理如圖4所示。其中λ<0是方向約束條件,即只有面對視線的片元才可能是可見的。

圖4 計算片元的可視性Fig.4 Estimates the visibility of the fragments

圖5給出了算法1的歸一化的中間數(shù)據(jù)及可視域計算結(jié)果,輸出結(jié)果中的紅色部分為不可視區(qū)域。所使用的地形數(shù)據(jù)為北京密云地區(qū)1 m分辨率DEM。圖5(b)、(c)中的某些區(qū)域存在鋸齒,對結(jié)果精度會產(chǎn)生一定影響,后面將給出偏移采樣法來處理。

圖5 算法1中間數(shù)據(jù)及輸出結(jié)果Fig.5 Intermediate data and result of scheme1

算法2:

第一步,計算

第二步,判斷片元處目標(biāo)的可見性。如果 d

圖6 計算片元處目標(biāo)的可見性Fig.6 Estimates the visibility of the targets above the fragments

3.3 減小鋸齒走樣的影響

距離圖的邊緣存在鋸齒走樣,如圖8(a)左上角所示。如果通過普通的圖像采樣方式來訪問距離圖,距離圖的邊緣鋸齒會使距離場映射結(jié)果產(chǎn)生條帶狀分布的鋸齒,如圖8(a)中線劃(方向參數(shù)λ=0,視線與地形表面相切)經(jīng)過的區(qū)域。這是由于用接近與地形表面相切的矢量vP(圖7中地形剖面的粗線部分)采樣距離圖時錯誤地獲取了遠處 PM點的dM值造成的。筆者通過一種偏移采樣的方法來減小鋸齒走樣對可視域計算影響圖8(b))。

設(shè)距離圖的視場角為θ,每個面的紋理分辨率為m×m,點 P在距離圖的視圖坐標(biāo)系下的視場深度為DP。則紋理空間中一個像素的寬度對應(yīng)物空間寬度為[16]

圖7 偏移采樣法原理圖Fig.7 Displaced sampling approach

圖8 使用偏移采樣法的前后對比Fig.8 Comparisons of outputs from sampling without and with displacement

在算法1和算法2中計算 dM時,如果檢測到距離圖水平邊緣,則令查詢矢量vP在垂直于邊緣的方向向下偏移τ,那么在距離圖中就產(chǎn)生一個像素的偏移量。算法原理如圖7所示,稱之為偏移采樣法,即算法3。

算法3:

第一步,計算觀察點到片元的方向矢量vP和向下偏移一個距離圖像素的矢量vτ的距離場映射值

第二步,如果滿足以下條件,則vP的距離場映射取為d′M

式中,δ是距離圖邊緣檢測的閾值。利用偏移采樣法,有效地減小了輪廓鋸齒的影響,輸出結(jié)果如圖8(b)所示。

4 試驗與討論

試驗在一臺裝有2 GHz CPU、1GB內(nèi)存、Geforce 9500GT顯卡的 PC上執(zhí)行。算法通過HLSL語言編寫的GPU程序?qū)崿F(xiàn)。

4.1 本文方法的結(jié)果

圖5已經(jīng)給出了算法1的輸出結(jié)果,這里給出了算法2應(yīng)用于全向可視域分析的結(jié)果,見圖9。試驗數(shù)據(jù)為北京密云地區(qū) 1 m分辨率DEM數(shù)據(jù)。

圖9 算法2的輸出結(jié)果Fig.9 Results of scheme 2

由于在生成距離圖時地形通過投影轉(zhuǎn)化為二維柵格數(shù)據(jù),可視域的計算與地形的幾何數(shù)據(jù)形式無關(guān),因此本文方法能同時適用于RSG以及TIN地形模型。另外,所提方法還可以用于處理含有地物的地形。地物可以記入DEM,也可以是獨立的幾何模型,只需在生成距離圖以及計算可視域的時候?qū)⒌匚锱c地形同時進行渲染即可。分析結(jié)果如圖10所示,圖中白點為觀察點。

圖10 包含建筑物場景的分析結(jié)果Fig.10 Results of scene which contains buildings

可以通過調(diào)整距離圖的分辨率(m×m)來控制分析結(jié)果的精度。由式(1)可推出,在觀察半徑為D的情況下,如果希望分析結(jié)果的空間分辨單元小于a,那么距離圖的分辨率應(yīng)滿足

所提方法的處理速度與地形數(shù)據(jù)量、距離圖分辨率的關(guān)系曲線如圖11所示。試驗表明,方法執(zhí)行速度非常快,對包含13萬個三角形的DEM進行半DEM網(wǎng)格精度的全向可視域分析,處理時間約為0.01秒。

圖11 本文方法處理速度曲線Fig.11 Graphs of the efficiency of our approach

4.2 對比試驗

相同試驗條件下,將不同m值時本文方法的結(jié)果與LOS方法[3]的結(jié)果進行了對比。以高精度LOS方法的輸出結(jié)果為參考圖像,圖12第一行為本文方法的輸出結(jié)果,第二行為本文方法結(jié)果與參考圖像的差值。本文方法的誤差主要來自距離圖采樣的鋸齒走樣問題,從試驗結(jié)果可以看出,通過提高距離圖分辨率可以減小誤差。

圖12 本文方法與LOS方法分析結(jié)果的對比Fig.12 Comparison of the result of our approach with LOS approach

將本文方法的處理速度與LOS方法[3]以及基于參考面[8]的方法作了對比。本文方法是GPU算法,后兩種算法都是CPU算法,進行對比的時間均為讀入DEM后直至提交結(jié)果顯示的時間。在完全相同的試驗條件下,對比結(jié)果如表1所示。

表1 處理速度對比Tab.1 Comparison of the efficiency

以上試驗中,本文方法在不同DEM數(shù)據(jù)量下都明顯快于其他兩種算法。與參考面方法相似,本文方法的處理速度也與可視域面積無關(guān)。

5 結(jié)論與進一步工作

本文給出一種基于距離場映射的地形可視域分析方法。首先,在觀察點位置O繪制地形,將O到地形表面可見點的距離作為輸出,保存到距離圖;然后,在頂視圖第二次渲染地形,將距離圖映射到地形表面,并通過與各片元到O的距離值進行對比來實現(xiàn)地形的可視性計算。本文方法基于GPU程序?qū)崿F(xiàn),能以很快的速度獲得高精度計算結(jié)果;能適用于多種形式的三維地形模型以及建筑物的可視性分析。

在進一步工作中,將結(jié)合軍事、建筑規(guī)劃、景觀分析等領(lǐng)域的具體需求,對本文方法進行擴展和深入研究。

[1] ZHOU Qiming,LIU Xuejun.Digital Terrain Analysis [M].Beijing:Science Press,2006:181-199.(周啟鳴,劉學(xué)軍.數(shù)字地形分析[M].北京:科學(xué)出版社,2006: 181-199.)

[2] DE FLORIANI L,MARZANO P,PUPPO E.Line-ofsight Communication on Terrain Models[J].International Journal of Geographical Information Systems,1994,8(4): 329-342.

[3] ZHU Qing,LIN Hui.Digital City Geographic Information System[M].Wuhan:Wuhan University Press,2004.(朱慶,林暉.數(shù)碼城市地理信息系統(tǒng)[M].武漢:武漢大學(xué)出版社,2004.)

[4] SHAPIRA A.Visibility and Terrain Labeling[D].Troy: Rensselaer Polytechnic Institute,1990.

[5] MILLS K,FOX G,HEIMBACH R.Implementing an Inter-visibility Analysis Model on a Parallel Computing System[J].Computers and Geosciences,1992,18(8): 1047-1054.

[6] TENG Y A,DE MENTHON D,DAVIS L S.Region-toregion Visibility Analysis Using Data Parallel Machines[J]. Concurrency:Practice and Experience,1993,5(5): 379-406.

[7] VAN KREVELD M.Variations on Sweep Algorithms: Efficient Computation of Extended Viewsheds and Class Intervals[C]∥Proceedings of Symposium on Spatial Data Handling.Delft:Taylor&Francis,1996,13A:15-27.

[8] FRANKLIN W R,RAY C K.Higher Isn’t Necessarily Better:Visibility Algorithms and Experiments[C]∥Proceedings of the 6th International Symposium on Spatial Data Handling,Advances in GIS Research.[S.l.]:Taylor &Francis,1994,(1-2):751-770.

[9] YING Shen,LI Lin,MEI Yang,et al.Characteristics Analysis ofIncrementalTerrain Visibility[J]. Acta Geodaetica et Cartographica,2007,36(2):192-197.(應(yīng)申,李霖,梅洋,等.增量法地形可視計算與分析[J].測繪學(xué)報,2007,36(2):192-197.)

[10] ZHOU Yang,TAN Bing,XU Qin.A Spatial Analysis Algorithm Based on 3D Topographic Map[J].Journal of Institute of Surveying and Mapping,2001,18(1):39-43. (周楊,譚兵,徐青.基于3維地形圖的空間分析算法[J].測繪學(xué)院學(xué)報,2001,18(1):39-43.)

[11] WU Yanlan.An Algorithm for Computing View Sheds Based on Reference Planes[J].Journal of Geomatics, 2001(1):19-21.(吳艷蘭.基于參考面的可視域算法[J].測繪信息與工程,2001(1):19-21.)

[12] YE Wei,TAO Yang.Construction and Application of Uniform Analysis Model for Terrain Visibility Based on DEMs[J].Geomatics World,2009,2(1):19-24.(葉蔚,陶旸.DEM地形可視性分析的統(tǒng)一模型構(gòu)建與應(yīng)用[J].地理信息世界,2009,2(1):19-24.)

[13] YING Shen.A Research on Key Techniques and Applications of Spatial Visibility Analysis[D].Wuhan:Wuhan University,2005.(應(yīng)申.空間可視分析的關(guān)鍵技術(shù)與應(yīng)用研究[D].武漢:武漢大學(xué),2005.)

[14] WILL IAMS L.Casting Curved Shadows on Curved Surface[C]∥Proceedings of the 5th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH′78.Atlanta:[s.n.],1978:270-274.

[15] GREENE N.Environment Mapping and Other Applications of World Projections[J].IEEE Computer Graphics and Applications,1986,6(11):21-29.

[16] COHEN J,OLANO M,MANOCHA D.Appearance Preserving Simplification[C]∥Proceedings of the 25th AnnualConference on Computer Graphics and Interactive Techniques,SIGGRAPH’98.New York:[s.n.],1998: 115-122.