一种基于SRTM的雷达可视域图绘制方法

赵 阳 李 伟 杨 申

(西安电子工程研究所 西安 710100)

0 引言

雷达在民用和军事方面的应用已经十分的成熟和广泛,其精确的定位能力使人们可以轻松掌握相距千里之外的目标信息。但是由于雷达依靠电磁波的回波来进行定位,而电磁波直线传播的特点导致其易受障碍物的遮挡,所以如何减小地物遮挡对雷达带来的影响至关重要[1]。在雷达选址时,除了选择地形开阔的阵地之外,还需要能够分析雷达周围一定范围内的可视域情况,进而能够避开遮挡严重的区域。

可视域是地理信息系统中的一个重要概念,是指从某一个固定位置处所能看到的区域范围,常常需要指定可视半径和观测的方位角。可视域的计算还与目标高度有关,是否可视指的是目标与观察点之间的连线上有无障碍物遮挡,常用的方法是比较目标和观察点之间连线的俯仰角与目标位置处遮蔽角的大小[2]。

本文实现了一种基于SRTM的雷达可视域图绘制方法,实现原理简单,地形数据精度高,算法通用性好,可应用于其他雷达系统。本文还将该方法的Matlab仿真结果与一款专业的地图绘制软件GlobalMapper进行了对比,验证了该方法的可行性和有效性。

1 SRTM简介

SRTM(Shuttle Radar Topography Mission)指的是航天飞机雷达地形测绘任务,是由美国航空航天局、地理空间情报局以及德国和意大利的航天机构于2000年2月开展的一项探测任务。该任务历时11天,通过发射奋进号航天飞机,使用其机载影像雷达的C波段和X波段获取了北纬60°至南纬56°之间,面积超过1.19亿平方公里的12TB雷达影像数据,覆盖全球陆地表面的80%以上。全部数据经过两年左右的处理,获取平面精度±20m,高程精度±16m的全球数字高程模型(DEM),可绘制成高精度的全球三维地形图[3]。

SRTM数据每经纬度方格提供一个文件,数据值为高程值,是纬度和经度方向上的二维采样结果。SRTM的精度有1 arc-second和3 arc-seconds两种,称作SRTM1和SRTM3,或者称作30M和90M数据[4]。SRTM1的文件里面包含3601×3601个采样点的高程数据,所以每个采样间隔经纬度增加1/3600°;SRTM3的文件里面包含1201×1201个采样点的高程数据,每个采样间隔经纬度增加1/1200°。本文仿真中使用的是SRTM3格式数据。

例如文件名为N27E086.hgt的SRTM3数据,采样的起始纬度为北纬27°,之后每个采样点的纬度增加1/1200°,终止纬度为北纬28°;采样的起始经度为东经86°,之后每个采样点的经度增加1/1200°,终止经度为东经87°。SRTM数据格式如图1所示。

图1 SRTM数据格式示例

2 雷达可视域图绘制流程

本文所述的基于SRTM的雷达可视域图绘制流程如图2所示。

图2 基于SRTM的雷达可视域图绘制流程

从图2中可以看出,在设置好输入参数以后就可以从SRTM文件中获取雷达位置所在区域的高程数据,根据前述可知获取到的SRTM数据范围是1°×1°的经纬度方块。然后从起始方位角开始以一定的步长(步长一般是1°)扫描整个扇区,循环进行某一个方位上的可视性计算,循环结束的条件是到达扇区的边界。在得到所有方位上的可视性后,绘制可视域图。

各步骤的详细设计在下文中进行描述。

2.1 设置输入参数

输入参数包含雷达参数和通视参数。雷达参数包含雷达所在位置的经度L0、纬度B0和高程H0,以及雷达天线的高度ha;通视参数包含可视半径Rmax、起始方位角θ0、扫描的扇区宽度θl、扫描的步长θs、选定点相对雷达的俯仰角φ、目标相对雷达的高度hr。

2.2 获取雷达位置所在区域的SRTM高程数据

根据雷达所处位置的经纬度可以解析出雷达位置所在区域的SRTM文件名,再根据文件名从SRTM地图数据文件夹中读取相应的SRTM数据,将获取到的数据的经度和纬度分别保存在一维向量中,高程值保存在二维矩阵中。

以我国的地形数据为例(北纬、东经),经度L0、纬度B0和文件名filename的关系为

filename=Nfloor(B0)Efloor(L0).hgt

(1)

其中,N代表北纬,E代表东经,floor()函数是向下取整函数,.hgt是文件后缀名。

例如位于东经108.332571°、北纬39.823789°位置的雷达,其SRTM数据所对应的文件名为N39E108.hgt。

高程数据在SRTM文件中的存储格式是16bit的signed integer,每两个字节以空格分隔,分别为高8位数据和低8位数据。因此,需要提取第一个字节乘以256,再与第二个字节相加,得到原始高程数据。

以上原始高程数据中还存在一些特殊数据,例如32768、65535等,属于异常值,对于这些数据需要进行异常值修正,修正后的矩阵才能用于后续计算。本文采用临近点的高程值来修正异常值,具体做法如下:

1)找出二维矩阵中的异常点坐标;

2)对于非第一列数据的异常值取其左侧相邻点赋值;

3)对于第一列数据的异常值,若非第一行则取其上一行相邻点赋值;

4)对于第一列第一行数据的异常值,取其相邻三个数据的平均值进行赋值。

2.3 某一个方位上的可视性计算

通常研究可视域方位角并非只有一个角度,而是一段扇形区域。针对固定的可视半径,可以将扇形区域的问题简化为某一个方位上的问题,只不过整个流程需要做多次循环计算。如图3所示为某一个方位上的可视性计算流程。设循环计数值为i,则该方位角θi的大小为式(2)所示。

图3 某一个方位角上的可视性计算流程

(2)

其中,各符号的含义参见2.1节。

2.3.1 根据选定点的RAE坐标计算LBH坐标

将每个方位上的选定点设为可视半径处的最远点,其RAE坐标指的是选定点相对雷达的距离、方位角和俯仰角坐标,根据2.1节和式(2)可知RAE坐标为(Rmax,θ,φ)。

LBH坐标指的是选定点的经度、纬度和高程坐标,其坐标转换过程如下:

1)由选定点的RAE坐标(Rmax,θ,φ)计算选定点的雷达站心直角坐标(rx,ry,rz)

(3)

2)由雷达的经纬度坐标(L0,B0,H0)和选定点的雷达站心直角坐标(rx,ry,rz),计算选定点的地心直角坐标(X,Y,Z)

(4)

(5)

(6)

3)采用迭代法,由选定点的地心直角坐标(X,Y,Z)计算选定点在WGS84坐标系下的经纬度坐标LBH

(7)

(8)

(9)

2.3.2 拼接SRTM数据

由2.3.1节中的计算步骤可以得到选定点的LBH坐标,根据L和B的值判断选定点是否超出当前SRTM数据的范围。若超出则重复2.2节中的计算过程,获取雷达所在位置到选定点之间所有区域的SRTM高程数据,并与原SRTM数据进行矩阵拼接;否则不进行数据拼接。

2.3.3 构建高程线方程

高程线指的是从雷达所在位置到选定点之间的连线,连线上各采样点处的高程值必须能够获取。本文采用等间隔采样法,每隔相同的距离采样一个点,可以得到采样点数目num为

(10)

其中,Rinter表示采样间隔,本文使用的间隔为100m。采样点数目与方位无关,这样可以保证各个方位上的采样点数目一致。接下来的关键问题是如何获取高程线上各采样点处的高程值。

根据第1节可知,SRTM数据也是由离散的采样点构成,单个文件所在区域的采样点数目是1201×1201(SRTM3格式)。又根据2.3.2节,即使当前SRTM数据的范围已经囊括整条高程线的数据范围,但是由于两者采样间隔的不同并不能保证经纬度重合,即高程线上采样点的经纬度不一定在SRTM区域中有对应的离散值,所以需要求高程线上采样点高程的近似值。

本文采用最近距离法获取高程线上各采样点高程的近似值:假设要获取高程线上采样点j的高程值,则在SRTM高程数据矩阵中,选取与点j经纬度距离最近的点的高程值作为点j的高程值。根据此方法可以近似获得高程线上所有采样点的高程值。

注意,此时获得的高程值并没有考虑地球曲率对距离带来的影响,所以还需要借助地球曲率修正算法,对所有高程值进行曲率修正。

(11)

其中:Hcorr为曲率修正后的高程值;Hori为修正前的原始高程值;R0为地球的等效半径;Rj为采样点j到雷达的水平距离。

2.3.4 计算高程线上各个采样点的最大遮蔽角

遮蔽角的含义是雷达天线顶端与障碍物顶端的连线和地平线之间的夹角[5],如图4所示,两条虚线之间的夹角αj就是一个遮蔽角。当目标与天线顶端的连线和地平线之间的夹角(也叫目标的俯仰角)小于该遮蔽角时目标就会被遮挡,所以有遮蔽的含义。

图4 采样点j处的遮蔽角示意图

根据图4所示,可以得出高程线上第j个采样点的遮蔽角αj的计算公式为

(12)

其中:Hj表示该采样点处曲率修正后的高程值;H0表示雷达的高程;ha表示雷达天线的高度;Rj表示该采样点到雷达的水平距离。

而某个采样点位置的最大遮蔽角,指的是该采样点以及其前面所有采样点的遮蔽角的最大值。如图5所示,αA是采样点A处的遮蔽角,αB是采样点B处的遮蔽角,虽然HB>HA,但是有αB<αA。可以看到在这种情况下,障碍物B已经被障碍物A遮挡,所以B点处的最大遮蔽角是更大的αA。

图5 最大遮蔽角示意图

根据以上描述可以获得最大遮蔽角的计算方法:使用一个全局变量保存最大遮蔽角记为αmax,然后循环遍历每个采样点并计算其遮蔽角αj,将每次计算出的αj与αmax作比较并更新αmax为较大值,则该采样点的最大遮蔽角αMj就是此时的αmax。按照此方法,图5中A点的最大遮蔽角为αA,B点的最大遮蔽角也是αA。

可以看到随着采样点的增多最大遮蔽角在不断更新为更大的遮蔽角,这样就可以保证每个采样点处的最大遮蔽角都是其前面所有采样点遮蔽角的最大值。最终的计算结果为一个向量,保存着每个采样点处的最大遮蔽角。

2.3.5 计算目标高度下各个采样点处的可视性

有了每个采样点的最大遮蔽角,就可以计算出目标高度下每个采样点处的可视性。根据2.3.4节中遮蔽角的含义可得,目标高度下在某一采样点处,当雷达与目标之间连线的俯仰角大于该采样点的最大遮蔽角时,该采样点处为可视,反之不可视。第j个采样点处的可视性vj表示为

(13)

其中:αTj为目标在相对雷达hr高度下该采样点处的俯仰角;αMj为该采样点处的最大遮蔽角。根据几何关系有式(14)为

(14)

其中:H0表示雷达的高程;ha表示雷达天线的高度;Rj表示该采样点到雷达的水平距离。示意图如图6所示。

图6 目标相对雷达hr高度下采样点j处的可视性

2.4 绘制可视域图

在扇区内各个方位上重复2.3节的计算过程,就能得到hr高度下各方位上高程线各个采样点处的可视性,据此就可以绘制出可视域图。

本文使用Matlab的polarplot()函数绘制可视域图。polarplot()函数的作用是绘制极坐标图,传入的第一个参数是角度向量,第二个参数是待绘制参数的向量。

绘制方法如下:

2)遍历各个采样点,采样点范围从1到num,其中num可由式(10)得到;

对于每个采样点j,设其可视距离向量为rangeVector,向量长度与angleVector相同。求rangeVector的步骤如下:

①遍历angleVector中的各个方位角θi

ⅰ.取方位角θi上,高程线上采样点j处的可视性vij;

ⅱ.取方位角θi上,高程线上采样点j距雷达的水平距离Rij;

ⅲ.可视距离向量的第i个元素rangeVectori=vijRij,将rangeVectori添加到rangeVector中。

②遍历完angleVector后就会得到采样点j处的可视距离向量rangeVector,使用polarplot(angleVector,rangeVector)函数绘制采样点j处的单个圆;

3)遍历完所有采样点后,此时就绘制了num个圆。可视域图绘制结束。

图7为可视域图绘制原理示意图。雷达位于点O,可视半径设为600m,在5个方位上全部都每隔100m采样6个点,形成射线OA到OE。

图7 可视域图绘制原理示意图

步骤1)中的angleVector的取值为射线OA到OE的5个方位角的大小。

步骤2)对6个采样点进行遍历,每个采样点都有自己的rangeVector,其长度都为5(图中离散的弧线)。假设图中所有的点都可视,则采样点1的rangeVector值为{100,100,100,100,100},采样点2的rangeVector值为{200,200,200,200,200},其他采样点类似。对于每个采样点都绘制出离散的圆弧。

步骤3)可以绘制出整个图7。

3 仿真结果与分析

本文使用Matlab绘制出可视域图之后,还与一款专业的地图绘制软件GlobalMapper做了比较。

3.1 GlobalMapper软件简介

GlobalMapper软件是GlobalMapper软件公司开发的一款地图绘制软件,支持多种数据源的输入和读取(包含SRTM),可以处理栅格和矢量数据。软件预置了非常多的常用坐标系和转换参数,提供了各种坐标转换方法,支持数百种地图投影。此外还拥有强大的地图编辑、转换、打印等功能[6]。

本文所述的方法需要参考GlobalMapper软件中的视图分析工具,其界面如图8所示。

图8中:“发射点高度(地面之上)”指的是雷达天线的高度ha;“接收点高度(地面之上)”指的是目标相对雷达的高度hr;“视角”参数中的“起始角度”指的是起始方位角θ0,“扫描的角度”指的是扫描的扇区宽度θl,扫描的步长θs默认为1°;选定点相对雷达的俯仰角φ默认为0;“可视半径”含义相同指的是Rmax;而雷达所在位置的经纬高参数L0、B0和H0在鼠标点击地图时可以由软件自动获取。

3.2 仿真结果对比分析

本文通过两个典型的示例,分别在雷达仰视和下视的情况下对比了Matlab的仿真结果和GlobalMapper软件的视图分析结果。

1)雷达仰视的情况

参数如下:

L0=78.943321°
B0=13.165818°

H0=455m
ha=7m
Rmax=50km
θ0=0°
θl=360°
θs=1°
φ=0°
hr=500m

仿真结果如图9所示。

图9中的右侧部分为GlobalMapper软件绘制出的可视域图,灰色阴影区域表示可视范围,其半径表示可视距离。可以看到左右两个结果的可视域轮廓基本一致,在方位90°～120°左右可视距离能够达到50km,视线良好,其余方位的可视距离不够50km。

2)雷达下视的情况

参数如下:

L0=79.270324°
B0=13.799420°
H0=940m
ha=7m
Rmax=50km

θ0=0°
θl=360°
θs=1°
φ=0°
hr=-500m

仿真结果如图10所示。

图10 雷达下视情况的仿真结果

可以看到此参数下可视域不连续。这是因为假设目标能够低空贴地飞行,其障碍物遮挡情况将变得复杂。

综上,本文所述的可视域图绘制方法能够很好地贴合于专业绘图软件GlobalMapper的绘制结果,证明了该算法的可实现性和有效性。

4 结束语

在雷达选址时可能遇到地物遮挡,因此需要绘制雷达站周围一定范围内的可视域图,针对于此,本文首先通过解析输入参数、获取SRTM高程数据和坐标转换建立各个方位上的高程线方程;然后计算出高程线上各个采样点处的最大遮蔽角;再计算出目标高度下高程线上各个采样点处的可视性;最后,使用Matlab软件绘制出了可视域图。本文还将Matlab仿真结果与地图绘制软件GlobalMapper进行了对比,验证了该方法的可行性和有效性。