1 引言
基于网络环境的三维虚拟地球软件推动了地理信息大众化发展。但在专业GIS应用方面，由于缺乏三维空间分析功能，三维虚拟地球软件往往还不具备实用性。视域分析作为重要的空间分析方法，广泛应用在旅游中的风景评价，房地产中视线遮挡判断，以及通讯中的信号覆盖，或军事上的火力覆盖等多方面[1]。在传统三维GIS软件中，视域分析基于通视分析进行，即通过判断给定对象与一组目标点之间的通视性，从而获得指定对象的视域可见范围。主要通视分析算法有以下4种：JANUS算法，DYNTACS(Dynamic Tactical Simulation)算法，ModSAF算法和Bresenham算法[2]。以上方法都是基于平面基准，未考虑地球曲率的影响。而基于三维虚拟地球的空间视域分析，地球曲率对分析结果影响很大。例如，美军“宙斯盾”系统的核心SPY-1相控阵雷达，最大探测距离超过300千米，由于地球曲率的影响，对低空掠海10米飞行目标的探测距离却只有约31千米[3]。
对此，本文主要探讨了三维虚拟地球中基于球面基准下的视域分析方法，对JANUS算法进行改进，对全球地形数据采用金字塔分层分块数据组织基础上，顾及地球曲率条件，实现大范围视域分析，同时基于NSCGlobe平台将分析结果进行了可视化。
2 三维虚拟地球中实时视域分析方法
2.1视域分析中地形数据组织
在三维虚拟地球中，构建三维场景的地形数据具有海量特性。因此，在进行视域分析时，无法将大范围地形数据一次性载入内存进行计算分析。同时，对于视域分析，根据应用特点灵活采用不同尺度的地形，既不影响分析精度，还可以提高分析处理的效率。例如，对于地面火炮系统，采用较高分辨率的地形有利于精确打击。而对于那些与细节模型关系不大的应用中，高分辨率模型反而会影响分析处理效率。如对于离地10000米以上作战的飞机，分析其火力覆盖的视域时就可以采用较低分辨率的地形，以便火力覆盖的实时计算。
因此，根据视域分析的多样性，应该灵活采用多尺度地形进行分析计算。这是由于多尺度地形避免了信息冗余，并且易于支持任何细节层次的地形表达[3]。基于上述考虑，综合已有全球空间数据组织方法，在三维虚拟地球中进行视域分析，本文以金字塔结构分层分块数据组织方法对三维虚拟地球中海量存储的数字高程模型DEM(Digital Elevation Model)进行组织，其核心是面向全球范围，采用Ottoson提出的用椭球四叉树结构建立全球等经纬度的空间数据索引[4]，分别建立东西半球两个层次四叉树结构，每个层次四叉树结构中四个节点记录节点的空间范围，同时存储指向该范围的地形数据集指针，如图1。
图1全球范围地形金字塔结构图
基于上述数据组织方法，在进行全球范围视域分析时，根据视域分析特点，采用合适尺度的地形数据层。确定地形数据层后，根据视域覆盖，计算需要的地形块索引。然后根据索引调度数据，逐块进行视域分析，建立多尺度分块的地形数据的视域分析计算分析。
2.2全球基准下的视域分析方法
在传统基于平面基准的视域分析方法中，无需考虑地球曲率的影响。在全球基准下进行大范围视域分析，受到地球曲率影响明显。若将上述方法运用在全球大范围条件下，将导致球面下不可见的目标点判断成可见点，如图2、图3。
图2球面下A、B不能通视
图3展成平面后A、B可以通视
针对上述问题，本文对JANUS算法进行改进，在三维球面空间中进行视域分析。主要流程如下：
● 在三维虚拟地球中实时确定观察点经纬度B，L和高程Z，并转换为以地心为原点的直角坐标系统中V(x0，y0，z0);
● 根据观察点经纬度坐标(lat,lon)和分析半径R(分析半径)，计算需要分析的DEM块索引编码，计算方法如公式(1)：
地形块索引编码由Level_ID+ Block_ID两部分组成。式中，N为经验值，依显示屏幕分辨率而定，取值范围为3—6，R(Level_Zero)为0层(即最上层)地形块半径,tilesize为第Level_ID层地形块经纬度跨度。
● 根据地形块索引编码，将地形块逐块载入内存进行视域分析;
对每个DEM块进行视域分析时，遍历DEM块中各地形点，依次作为目标点T(x，y，z)，对观察点V与目标点T进行通视分析。本文在球面基准下进行分析，LOS(Line Of Sight)为V到T的视线段，判断V与T之间的可视性与否，即判断LOS上各点在地球表面投影点的高程是否低于LOS上对应点高程。对LOS上每一个点都进行通视分析显然不现实，为了同时兼顾分析的效率与精度，将LOS等分成若干段后取等分点进行分析。如图4所示，对第i个等分点进行分析，由于地心O(0,0,0)、第i个等分点Pi(xi,yi,zi)、OPi与地球球面交点Vi(Vix,Viy,Viz),OPi与地形表面交点M共线，M即为Pi在地球表面投影，Vi即为Pi、M在地球球面投影，地球半径为R，M高程为hm，Pi高程为hpi,比较hm与hpi大小即可，计算步骤如下：
图4球面基准下通视分析图
● 计算LOS的空间长度d;
● 用d除以所选择DEM层的空间分辨率，计算结果加1，取整得到整数k，由此使得k必大于0;
● 将LOS等分成k段;
● 依次遍历k个等分点，进行通视分析;
在第i个等分点计算中，先求出该等分点Pi(xi，yi，zi)坐标，按公式(2)计算：
从而可以根据空间距离公式计算出OPi长度，而i等分点高程为hpi，R为地球半径，计算得到：
同时，根据共线关系式(4)：
得到Vi点坐标：
利用公式(5)求出Vi点经纬度地理坐标。
利用Vi点经纬度地理坐标，根据NSCGlobe中等经纬度的全球空间数据索引，检索到相应地形数据文件，通过该DEM格网点高程值双线性内插得到M高程hm。如图5所示，灰色点为DEM格网点，h1, h2, h3, h4分别为其高程值
图5双线性内插计算hm
约定相邻格网点横向、纵向距离都为1，黑色点为M点在球面投影Vi点,根据经纬度确定Vi点与临近格网点横向、纵向距离，如其与左下角DEM点的横向距离为v，纵向距离为u，则M点高程由公式(7)计算:
hm = (1-u)*(1-v)* h1 + v*(1-u)* h2 + u*(1-v)* h3 + u*v* h4 (7)
最后，比较hm与hpi大小：若hm>hpi,则V与T之间不可通视;否则继续对下一个等分点分析，直到所有等分点都分析完毕，若所有等分点都可通视，则V与T之间可以通视。
3 实验
基于上述方法，本文采用全球30米分辨率中国地区STRM地形数据，构建8层金字塔结构地形数据集，进行了多尺度地形视域分析。实验中计算多尺度地形分析可见点百分比K与单个尺度可见点百分比进行比较，K由公式(8)计算得到。分别基于不同的视域分析半径进行实验。
K= 可视通视点数/总地形分析点数 (8)
从实验结果表1可以看出，与仅仅采用单个尺度地形数据进行视域分析相比，在保证分析精度前提下，提高了运行效率。
为了将分析结果在三维虚拟地球中进行可视化表达，实验中以地形块为单元，逐点进行分析，若地形点可见，则像素设置为青色，否则设置为红色。最后将全部视域分析结果组织成辐射状的射线数据集。在此基础上，基于三维虚拟地球平台软件NSC Globe，将分析结果进行动态加载，从而实现分析结果实时可视化。实验中对电视塔信号覆盖进行模拟，电视塔高100米，在有效半径3000米条件下，其信号可视域覆盖范围如图5。
图6电视塔信号可视域覆盖分析效果图
4 结束语
本文主要探讨了三维虚拟地球中基于多尺度地形数据的大范围视域分析方法。下一步研究将针对应用中专业模型的特点，例如通信基站覆盖中电磁波传播模型，将应用模型与三维视域分析集成，构建基于三维虚拟地球的专业应用分析。
参考文献：
[1] Leila De Floriani,Paola Magillo.Intervisibility on Terrains[EB/OL]. ftp://ftp.disi.unige.it/person/MagilloP/TEMP/gisbook99.pdf.
[2] 王智杰,邱晓刚,李革.RSG地形通视性快速算法设计.计算机仿真, 2004,21(12):92-95.
[3] 刘占荣. 宙斯盾作战系统结构分析[ J ].情报指挥控制系统与仿真技术, 2004, 26 (1) : 22 - 31.
[4] Ottoson.P,Hauska.H.2002,Ellipsoidal quadtrees for indexing of global geographical data,int. J.geographical information science,vol.16,no.3,pp.223-226.
[5] Leila De Floriani,Paola MAgillo.Representing the visibility structure of a polyhedral terrain through a horizon map[J].International Journal of Geographical Information Systems,1996.
作者简介
董鉥涛(1987年—)，男，汉族，湖北十堰人。硕士研究生，主要从事三维地理信息系统方向研究。
侯小波(1986年—)，男，汉族，河南南阳人。硕士研究生，主要从事三维地理信息系统方向研究。
声明：中国勘测联合网登载此文出于传递更多信息之目的，并不意味着赞同其观点或证实其描述，文章内容仅供参考。
