中国指挥与控制学会会刊 
可视域分析在雷达、通信基站信号覆盖范围计算,以及观察点设置等领域具有广泛的应用[1-2]。通过地理信息系统(GIS)将地图、数据库和空间分析三者有机结合,越来越多的应用开始采用GIS对可视范围进行直观处理和显示,在可视域分析中得到广泛应用[3]。
在雷达探测应用领域,文献[4]利用DEM高程数据计算得到了对空雷达电磁波存在的遮挡盲区,文献[5]利用美国航天飞机测绘任务(SRTM)数字高程数据构建我国气象雷达网的威力评估系统,但上述研究均对可视域分析建模采用简化处理,即没有考虑障碍物后的可视区域,如图1 中的遮挡盲区后面的“可视部分”。在对海面可视域分析时,一般当辐射源架设在沿海陆地较高位置,此时海上岛屿等遮蔽物后面仍存在可视区域,须对这种情况进行定量分析建模。
在可视域算法方面,文献[1]将可视域判断算法分为比较LOS视线仰角、视线高度和视线斜率等几种,这几种算法本质上是一致的,直接计算时间复杂性为o(n2),n为数据分辨力。采用逐点计算方法的商用软件(如超图、MapGIS等软件)带有通视分析的功能,当数据分辨力较高且在大面积计算时,单点通视的计算时间一般为分钟量级,多点通视分析达到小时甚至天的时间量级,难以满足实时性要求。
为提高现有算法的计算速度,文献[6]提出一种基于环形扫描环的并行可视域分析算法,但需要较大存储量。文献[7]对比逐点计算(R3算法)、复用向内最近点法(R2算法)和复用向外逐点(XDraw)等三类方法,指出Xdraw算法是具有高性能低计算复杂度的算法。文献[8]从减少数据库访问的角度,针对Xdraw算法进行了改进可提高算法的实现速度,但并没有解决算法复杂性问题。
笔者注意到,上述相关算法均采用恒定的计算步长,在实际应用中,由于盲区面积区域和可视面积较为悬殊,如对海观测,盲区面积较小。如果对盲区和可视区域采用不同的计算步长,则可以在保持一定的计算精度基础上,大幅度减少计算量,从而提高算法计算效率。由于Xdraw算法相对其他算法具有实时性优势[7],所以本文针对高性能的Xdraw算法实现了可变步长算法。
1 基本概念
可视域分析依赖于两点的可视判断,即通视判断。地理上的两点可视是指两点之间不存在地物遮挡,在视线上可达到对方,称为通视。如果取视点的最大探测距离作为视线半径,对该长度半径的圆内的所有点进行可视域分析,得到可视点的集合即为该观察点的可视域。
在直角坐标系中,通常用两者连线所决定的视线高度来判断是否可视。对于观察点O和目标点A,其视线上任一点B的视线高度用下式计算:
h视线= (hr-ht)+ht
式中,hr和ht分别为观察点和目标点的高度,rOA和rOB分别为O点和A、B点的地面距离。则A点对O点可视表述为:当连线OA上任一点B的实际高度h小于其视线高度h视线时,A点和O点可视。
利用数字高程数据进行可视域分析时,首先需要确定数字高程模型,不同的模型决定了算法不同。数字高程模型是对地球表面地形的数字描述和模拟,最常用的是不规则三角形网(Triangulated Irregular Network,TIN)和规则网格模型(Regular Square Grid,RSG)[9]。本文针对DEM数据特点,采用应用较广的RSG模型。
2 数字高程数据处理
2000年2月11日,美国发射的“奋进”号航天飞机上搭载的SRTM系统获取了全球北纬60°至南纬60°之间总面积超过1.19亿km2的雷达影像数据,经过两年多的数据处理,制成数字地形高程模型(DEM),即SRTM地形产品数据。针对航天飞机机载雷达信号可能受到干扰,发生镜面反射、雷达阴影、回波滞后等问题导致的SRTM高程数据在水域、高山和峡谷地区存在小块的数据空缺点和空白区,Jarvis et al.(2008)整理并提供了更完整的SRTM数据[10]。最新V4.1版本由国际热带农业中心利用新的插值算法得到,更好地填补了SRTM数据的空洞。该数据可从中国科学院计算机网络信息中心国际科学数据镜像网站免费下载。
SRTM高程数据有SRTM1和SRTM3两种,对应地形分辨率为30m和90m,均为栅格数据。目前仅公开SRTM3数据,每5°经纬度存储为1个数据文件,文件名为srtm-XX-YY.zip,XX表示列数(01-72),YY表示行数(01-24),包括img、tif和ascii三种格式。这里采用tif格式数据,用16位有符号整型数据表示每个地理点的高度值。我周边沿海数据文件包括srtm-58-08、srtm-58-09、srtm-59-08、srtm-59-09、srtm-60-04-srtm-60-09、 srtm-61-04、srtm-61-09。为得到连续的地形高度数据,我们利用GIS软件的镶嵌功能把上述数据块进行拼接,得到我国海域连续数字地形数据文件。由于SRTM主要关注陆地地形数据,数据默认将海面高度设置为无效数据,在16位图文件中用-9999表示。为实现对海面可视域分析,我们通过编程对数据文件的海面高度值改为0。
对进行大范围内可视域分析时,需要考虑地球面弯曲引起的地面相对观测点的高度变化,对远距离观察点进行修正[11]。如图2 所示,距离观测点O较近的B点可认为和观测点在同一海拔上,而距离较远的A点,由于地球曲率的影响,其高度进行如下修正:
hΔ=
式中,KRe为等效地球半径,K为折射系数,标准大气下约为4/3,Re为地球平均半径,取6370km。
3 变步长的XDraw算法
可变步长的快速算法的基本思想是:对于盲区内和盲区外采用不同步长。对海分析由于盲区面积少,所以这里盲区内采用高分辨率步长,盲区外采用低分辨率步长。
鉴于XDraw算法是比较准确且CPU时间消耗少的算法,这里以XDraw算法为例进行介绍。XDraw算法将不同方位上相同距离的点组成一个观察环,并向外复用结果。当要计算i环上网格点的可视性时,利用i-1环的相邻网格点上记录的最低视线来判断该点的可视性,当高于内环视线则为可视,并记该点的上空最低视线为到该点的可视视线;否则不可视,记该点的上空最低视线为内环网格点的上空最低视线。我们根据XDraw算法复用思想,比较观察点对目标点斜率值(用绝对值表示)。
当观察点高度高于目标点时(对海观察往往如此),两者连线上不同距离采样点高度上任一点相对观察点的斜率大于前一点的斜率,则该点被前一点遮挡;当该点斜率小于前一观察点时,并不能确定该点可视,还需要跟该点与观察点连线上之前所有点的斜率作比较,当该点斜率均小于其他点斜率时,则该点可视。如图3 所示,观察点O(高度为h)对海波束向下俯视,最小可视斜率(绝对值)为h/R。沿着观察方向等距离间隔取点,其高度值分别为A、B、C、D,…,相应视线为LOS1、LOS2、LOS3、LOS4,…,当任一点视线斜率小于最小可视斜率(如视线LOS4)时,该点以远的所有点均不可视。当视线斜率大于最小可视斜率时,再与该点前面所有点中最大视线斜率相比,小于所有点的最大视线斜率时可视(如LOS3小于LOS1和LOS2的斜率,C点可视);否则不可视(如LOS2大于LOS1的斜率,B点不可视)。最后,在不同方位重复上述过程,得到所有可视点,即为给定视线半径下的可视域。
基于上述思想的计算流程如下:
1) 初始化,确定视点O(x,y)的高度h,视线半径R,方位采样步长a,距离采样步长r,计算视线最远处的斜率SR=h/R,即为最小可视斜率;
2) 取初始方位a=0°,获得距离r01 r02…r0k…R上的高度值h0k;将盲区标志初始化置0;
3) 用S表示斜率,计算目标点对观察点的斜率S0k=|h-h0k|/r0k(注意,如果目标距离视点较远,需要按照式(2)进行可视域修正),如S0k≤SR,则表示后续海面点全部不可视,转向6),该点的盲区标志置1;若S0k>SR,则进一步比较当前点斜率S0k和前一点斜率S0k-1;
4) 如S0k>S0k-1则为盲区点,再判断是否在盲区,如果盲区标志为0,表明开始进入盲区,盲区标志置1,并缩小距离和方位仿真步长。如果盲区标志为1,则表明已经进入盲区,则步长保持不变。
5) 如果S0k≤S0k-1,则和该点前面所有点最小斜率比较,如S0k<Min(S0m,m=0…k-1),该点为可视点,更新Min(S0m,m=0…k-1)。再判断是否在盲区,如果盲区标志为1,则该点为盲区结束点,恢复计算步长,否则该点在可视区,计算步长保持不变。
6) 该方位线上所有可视点连线将构成多个可能不连续的线段,将每条线段的开始点和结束点坐标分别存入数据点集合StaPointList和EndPointList中,用于可视域的表达;
7) 距离步长k遍历,重复3)-6),直到r=R;
8) a=a+a方位步长增加,重复2)-7),直到a=360°;
9) 绘制每个方位线上的所有可视线段,得到观察点在给定视线半径下的可视域。
4 仿真结果
硬件环境:双频2.6GHz的处理器和4GHz内存的单机;软件环境:基于GIS开发组件的VS2008编程环境。
误差计算:我们将变步长计算结果和定步长算法结果进行比对。参考文献[12]将算法得到的可视点数之差占可视点比例衡量误差大小。由于我们算法的可视域是用从观察点到目标点线段组成,所以取可视域面积交集占实际可视面积(以GIS软件可视域分析结果为真值)的比例为分析准确度。
设软件自带可视域分析得到的面积为S,而我们的方法得到可视域面积为S'(采用不同长度线段乘以相应的方位步长作为面积元,累积后可求得),则准确度定义如下:
K= ×100%
这里选择批山岛海域O点(121.3145°E,28.3301°N,h=435m)为传感器观察点,对360°全方位, 60km(47海里)半径区域,采用低分辨率固定步长、高分辨率固定步长、高分辨可变步长(为简单起见,这里仅对距离进行变步长采样,方位采样恒为a=0.1°,算法在盲区采用大步长,可视区采用小步长),得到可视化分析结果、运行时间和准确度如表1 所示。
表1 不同可视域分析方法的性能比较 |
| 参数 设置 | 固定步长 (rΔ=1000m) | 固定步长 (rΔ=100m) | 变步长 (rΔ=1000m,rΔ'=100m) | 逐点计算 (地形分辨率90m) |
|---|---|---|---|---|
| 仿真 结果 |  |  |  |  |
| 耗时 | 46ms | 515ms | 61ms | 551ms |
| 准确度 | 72% | 91% | 87% | 100% |
从表1 中可以得出,我们所提的变步长方法时间消耗和低分辨率固定步长相似,比高精度固定步长减小约1个数量级,而误差则增加不大(与高精度相比只有4%)。另外,所提算法可根据可视域分析区域的大小,在计算精度和时间之间折中。随着多观察点的加入,该算法的速度提升将更加明显,可满足单机版实时处理的需要。
5 结束语
针对实际系统对海可视域分析的需求以及当前算法存在的问题,本文提出了一种对海快速可视域快速分析方法。通过和GIS软件自带的可视域分析结果进行比较,本文所提方法在保持一定计算精度情况下,可得到数量级的计算速度提升。实际系统建模中,由于光波和微波高端波段的反射和绕射可忽略,所以可直接应用本文的方法。而当传感器的波长较低和应用的精度较高的场合,还需要进一步需考虑电磁波绕射和多径效应。