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

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

1 引 言

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

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

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

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

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

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

2 相關(guān)研究

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

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

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

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

3 本文方法

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

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

3.1 生成距離圖

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

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

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

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

3.2 計(jì)算可視域

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

算法1:

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

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

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

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

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

算法2:

第一步,計(jì)算

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

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

3.3 減小鋸齒走樣的影響

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

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

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

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

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

算法3:

第一步,計(jì)算觀察點(diǎn)到片元的方向矢量vP和向下偏移一個(gè)距離圖像素的矢量vτ的距離場(chǎng)映射值

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

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

4 試驗(yàn)與討論

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

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

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

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

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

圖10 包含建筑物場(chǎng)景的分析結(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所示。試驗(yàn)表明,方法執(zhí)行速度非常快,對(duì)包含13萬個(gè)三角形的DEM進(jìn)行半DEM網(wǎng)格精度的全向可視域分析,處理時(shí)間約為0.01秒。

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

4.2 對(duì)比試驗(yàn)

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

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

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

表1 處理速度對(duì)比Tab.1 Comparison of the efficiency

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

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

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

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

[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ì)算與分析[J].測(cè)繪學(xué)報(bào),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].測(cè)繪學(xué)院學(xué)報(bào),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].測(cè)繪信息與工程,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.