基于距离场映射的地形可视域分析方法

王晨昊,汤晓安,马伯宁,冷志光

1.国防科学技术大学,湖南长沙410073;2.中国人民解放军91635部队,北京102249

1 引 言

地形可视性也称为地形通视性(visibility),是指从一个或多个位置所能看到的地形范围或与其他地形点之间的可见度。可视性分析是地形分析的重要组成部分,很多与地形有关的问题都涉及地形可视性的计算问题。例如,观察站、雷达位置、广播电视或电话发射塔的位置、路径规划、航海导航、军事上的阵地布设、道路和建筑物的景观设计、日照分析等[1]。

点点通视(inter-visibility)分析和可视域(viewshed)分析是地形可视性分析的两个基本内容[2]。而可视域的计算从本质上讲,是两点之间的通视计算在面域上的实现,也就是说,给定视点的可视域就是通过视点的无数视线上的可见部分的集合[1]。目前已有的可视域分析算法全部是基于CPU实现的,算法复杂度较高,难以应付大数据量地形的在线分析;此外,这些算法几乎都是针对固定格式地形数据提出的,通用性较差。

针对以上问题,本文提出一种基于距离场映射的可视域分析方法。这种方法面向单观察点地形可视域计算问题,可以在三维地形渲染输出的同时获取地形的可视范围。与现有方法相比,这种方法具有以下特点:

(1)算法通过 GPU实现,计算速度更快。

(2)可以通过调整距离图分辨率来控制分析结果的精度,在处理速度和结果精度的平衡上具有较高灵活度。

(3)能同时适用于规则格网(RSG)和非规则格网(TIN)地形模型,也能用于地面物体的可视性分析。

2 相关研究

现有的地形可视域计算方法大部分是针对规则格网DEM提出的[3]。文献[4]提出一种基于视线(line of sight,LOS)的方法,跟踪从视点 O到任何其他点P的视线,并沿着O与P之间的射线移动,在到达 P之前检测到视线与地形相交或到达P,就终止移动。基本的LOS算法存在冗余计算,为此研究者对LOS方法做了许多改进,以期提高计算效率。相关成果有文献[5]提出的基于并行计算系统的可视性分析模型、文献[6]利用相邻视点之间的内在一致性提出的改进并行算法、文献[7]提出的平面扫描算法、文献[8]提出的利用辅助格网的改进算法、文献[9]提出的增量地形可视性算法以及文献[10]提出的基于斜率比较的可视性算法等。这些改进的LOS算法虽然一定程度上提高了执行速度,但由于必须内插计算与视线相交的DEM网格单元,计算量仍然比较大。

不同于LOS方法,文献[11]提出一种参考面法,通过视点和目标点之间的空间关系所形成的参考面来判断观测点与所有目标点是否可视。虽然这种方法比任何基于视线的算法都快,但仍然存在实现复杂度高、通用性差的问题。文献[12]则在分析可视性问题共性的基础上构建了一种DEM地形可视性统一分析模型,这个模型使得具体的可视性问题具有模式化的解决方法,但并未从速度、精度等方面对可视性算法进行改进。

面向TIN的可视域计算方法目前主要有前后顺序法、分而治之法、扫描线法、光线投射法等[13]。对比面向规则格网DEM的可视性分析方法和面向TIN的方法,各种算法往往只针对一种形式的地形模型,而不能通用。另外,对于包含建筑物、植被等地物的地形区域,现有分析方法的适用性也受到挑战。

本文方法受到文献[14]提出的阴影图方法的启发。阴影图方法用于真实感图形绘制中物体阴影的实时生成。由于研究对象以及应用侧重点不同,本文方法与阴影图方法的不同之处在于:①提出一种全新的地形可视域分析方法,根据应用中是否考虑目标高度,给出了两种不同的可视性判断算法;②针对地形距离图的特点,提出一种偏移采样法来提高可视性计算的精度。

3 本文方法

本文通过距离场映射的方法进行地形的可视性分析。主要包含两个步骤:第一步,在观察点以透视投影渲染地形,生成距离图;第二步,在顶视图以平行投影渲染地形,将距离图映射到地形表面片元,并进行片元可视性的判断。算法的执行包含对地形的两次渲染,流程如图1所示。

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

3.1 生成距离图

三维地形一般以表面模型的形式进行渲染,经过深度测试(depth test),显示设备可以获得每条投影射线上与视点最近的片元,通常情况下将其颜色值写入颜色缓存作为一个像素的显示输出,如图2(a)所示。如果不将片元颜色作为输出,而是将片元到观察点的距离值输出到颜色缓存,那么就可以获得一个以观察点为基准点的地表距离场在屏幕上的二维投影,称为距离图,记为fM。距离图记录了从O出发的视线与地形表面的第一个交点到O的距离值,如图2(b)所示。

图2 地形渲染输出Fig.2 The rendering outputs

当关注的地形相对于视点的角度较大,甚至是全向时,仅仅利用一个二维纹理将难以满足要求。这种情况下利用立方体贴图[15](cube map)来存储全向的距离场信息,如图3所示。立方体贴图包含六个面,每个面覆盖90°的视场角。在支持Direct3D 10的显卡上,利用几何着色器(geometry shader)只需通过对场景的一次渲染就可以获得立方体贴图。立方体贴图通过一个三维矢量来访问。

图3 立方体距离图Fig.3 Distance cube map

3.2 计算可视域

通过在观察点O处渲染地形,获得了距离图(立方体距离图)fM。接下来通过对地形的第二次渲染来获取地形可视域分析结果。为了能覆盖关注的地形区域,通过自顶向下的平行投影将地形渲染到目标缓冲区。经过栅格化处理之后,地形的几何数据被转换为片元,每个片元都对应于目标缓存的一个像素。设世界坐标系下观察点O到片元点 P的矢量为vP,片元法矢量为 nP。在像素着色器中通过以下算法进行片元的可视性判断。

算法1:

第一步,计算O到 P的欧氏距离d、片元的距离场映射 dM以及片元方向参数λ

第二步,判断片元的可视性。如果d=dM,且λ<0,那么相对于观察点,片元可见;否则,片元不可见,原理如图4所示。其中λ<0是方向约束条件,即只有面对视线的片元才可能是可见的。

图4 计算片元的可视性Fig.4 Estimates the visibility of the fragments

图5给出了算法1的归一化的中间数据及可视域计算结果,输出结果中的红色部分为不可视区域。所使用的地形数据为北京密云地区1 m分辨率DEM。图5(b)、(c)中的某些区域存在锯齿,对结果精度会产生一定影响,后面将给出偏移采样法来处理。

图5 算法1中间数据及输出结果Fig.5 Intermediate data and result of scheme1

算法2:

第一步,计算

第二步,判断片元处目标的可见性。如果 d

图6 计算片元处目标的可见性Fig.6 Estimates the visibility of the targets above the fragments

3.3 减小锯齿走样的影响

距离图的边缘存在锯齿走样,如图8(a)左上角所示。如果通过普通的图像采样方式来访问距离图,距离图的边缘锯齿会使距离场映射结果产生条带状分布的锯齿,如图8(a)中线划(方向参数λ=0,视线与地形表面相切)经过的区域。这是由于用接近与地形表面相切的矢量vP(图7中地形剖面的粗线部分)采样距离图时错误地获取了远处 PM点的dM值造成的。笔者通过一种偏移采样的方法来减小锯齿走样对可视域计算影响图8(b))。

设距离图的视场角为θ,每个面的纹理分辨率为m×m,点 P在距离图的视图坐标系下的视场深度为DP。则纹理空间中一个像素的宽度对应物空间宽度为[16]

图7 偏移采样法原理图Fig.7 Displaced sampling approach

图8 使用偏移采样法的前后对比Fig.8 Comparisons of outputs from sampling without and with displacement

在算法1和算法2中计算 dM时,如果检测到距离图水平边缘,则令查询矢量vP在垂直于边缘的方向向下偏移τ,那么在距离图中就产生一个像素的偏移量。算法原理如图7所示,称之为偏移采样法,即算法3。

算法3:

第一步,计算观察点到片元的方向矢量vP和向下偏移一个距离图像素的矢量vτ的距离场映射值

第二步,如果满足以下条件,则vP的距离场映射取为d′M

式中,δ是距离图边缘检测的阈值。利用偏移采样法,有效地减小了轮廓锯齿的影响,输出结果如图8(b)所示。

4 试验与讨论

试验在一台装有2 GHz CPU、1GB内存、Geforce 9500GT显卡的 PC上执行。算法通过HLSL语言编写的GPU程序实现。

4.1 本文方法的结果

图5已经给出了算法1的输出结果,这里给出了算法2应用于全向可视域分析的结果,见图9。试验数据为北京密云地区 1 m分辨率DEM数据。

图9 算法2的输出结果Fig.9 Results of scheme 2

由于在生成距离图时地形通过投影转化为二维栅格数据,可视域的计算与地形的几何数据形式无关,因此本文方法能同时适用于RSG以及TIN地形模型。另外,所提方法还可以用于处理含有地物的地形。地物可以记入DEM,也可以是独立的几何模型,只需在生成距离图以及计算可视域的时候将地物与地形同时进行渲染即可。分析结果如图10所示,图中白点为观察点。

图10 包含建筑物场景的分析结果Fig.10 Results of scene which contains buildings

可以通过调整距离图的分辨率(m×m)来控制分析结果的精度。由式(1)可推出,在观察半径为D的情况下,如果希望分析结果的空间分辨单元小于a,那么距离图的分辨率应满足

所提方法的处理速度与地形数据量、距离图分辨率的关系曲线如图11所示。试验表明,方法执行速度非常快,对包含13万个三角形的DEM进行半DEM网格精度的全向可视域分析,处理时间约为0.01秒。

图11 本文方法处理速度曲线Fig.11 Graphs of the efficiency of our approach

4.2 对比试验

相同试验条件下,将不同m值时本文方法的结果与LOS方法[3]的结果进行了对比。以高精度LOS方法的输出结果为参考图像,图12第一行为本文方法的输出结果,第二行为本文方法结果与参考图像的差值。本文方法的误差主要来自距离图采样的锯齿走样问题,从试验结果可以看出,通过提高距离图分辨率可以减小误差。

图12 本文方法与LOS方法分析结果的对比Fig.12 Comparison of the result of our approach with LOS approach

将本文方法的处理速度与LOS方法[3]以及基于参考面[8]的方法作了对比。本文方法是GPU算法,后两种算法都是CPU算法,进行对比的时间均为读入DEM后直至提交结果显示的时间。在完全相同的试验条件下,对比结果如表1所示。

表1 处理速度对比Tab.1 Comparison of the efficiency

以上试验中,本文方法在不同DEM数据量下都明显快于其他两种算法。与参考面方法相似,本文方法的处理速度也与可视域面积无关。

5 结论与进一步工作

本文给出一种基于距离场映射的地形可视域分析方法。首先,在观察点位置O绘制地形,将O到地形表面可见点的距离作为输出,保存到距离图;然后,在顶视图第二次渲染地形,将距离图映射到地形表面,并通过与各片元到O的距离值进行对比来实现地形的可视性计算。本文方法基于GPU程序实现,能以很快的速度获得高精度计算结果;能适用于多种形式的三维地形模型以及建筑物的可视性分析。

在进一步工作中,将结合军事、建筑规划、景观分析等领域的具体需求,对本文方法进行扩展和深入研究。

[1] ZHOU Qiming,LIU Xuejun.Digital Terrain Analysis [M].Beijing:Science Press,2006:181-199.(周启鸣,刘学军.数字地形分析[M].北京:科学出版社,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.(朱庆,林晖.数码城市地理信息系统[M].武汉:武汉大学出版社,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.(应申,李霖,梅洋,等.增量法地形可视计算与分析[J].测绘学报,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].测绘学院学报,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地形可视性分析的统一模型构建与应用[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.(应申.空间可视分析的关键技术与应用研究[D].武汉:武汉大学,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.