趙 陽 李 偉 楊 申

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

0 引言

雷達在民用和軍事方面的應用已經(jīng)十分的成熟和廣泛,其精確的定位能力使人們可以輕松掌握相距千里之外的目標信息。但是由于雷達依靠電磁波的回波來進行定位,而電磁波直線傳播的特點導致其易受障礙物的遮擋,所以如何減小地物遮擋對雷達帶來的影響至關(guān)重要[1]。在雷達選址時,除了選擇地形開闊的陣地之外,還需要能夠分析雷達周圍一定范圍內(nèi)的可視域情況,進而能夠避開遮擋嚴重的區(qū)域。

可視域是地理信息系統(tǒng)中的一個重要概念,是指從某一個固定位置處所能看到的區(qū)域范圍,常常需要指定可視半徑和觀測的方位角。可視域的計算還與目標高度有關(guān),是否可視指的是目標與觀察點之間的連線上有無障礙物遮擋,常用的方法是比較目標和觀察點之間連線的俯仰角與目標位置處遮蔽角的大小[2]。

本文實現(xiàn)了一種基于SRTM的雷達可視域圖繪制方法,實現(xiàn)原理簡單,地形數(shù)據(jù)精度高,算法通用性好,可應用于其他雷達系統(tǒng)。本文還將該方法的Matlab仿真結(jié)果與一款專業(yè)的地圖繪制軟件GlobalMapper進行了對比,驗證了該方法的可行性和有效性。

1 SRTM簡介

SRTM(Shuttle Radar Topography Mission)指的是航天飛機雷達地形測繪任務,是由美國航空航天局、地理空間情報局以及德國和意大利的航天機構(gòu)于2000年2月開展的一項探測任務。該任務歷時11天,通過發(fā)射奮進號航天飛機,使用其機載影像雷達的C波段和X波段獲取了北緯60°至南緯56°之間,面積超過1.19億平方公里的12TB雷達影像數(shù)據(jù),覆蓋全球陸地表面的80%以上。全部數(shù)據(jù)經(jīng)過兩年左右的處理,獲取平面精度±20m,高程精度±16m的全球數(shù)字高程模型(DEM),可繪制成高精度的全球三維地形圖[3]。

SRTM數(shù)據(jù)每經(jīng)緯度方格提供一個文件,數(shù)據(jù)值為高程值,是緯度和經(jīng)度方向上的二維采樣結(jié)果。SRTM的精度有1 arc-second和3 arc-seconds兩種,稱作SRTM1和SRTM3,或者稱作30M和90M數(shù)據(jù)[4]。SRTM1的文件里面包含3601×3601個采樣點的高程數(shù)據(jù),所以每個采樣間隔經(jīng)緯度增加1/3600°;SRTM3的文件里面包含1201×1201個采樣點的高程數(shù)據(jù),每個采樣間隔經(jīng)緯度增加1/1200°。本文仿真中使用的是SRTM3格式數(shù)據(jù)。

例如文件名為N27E086.hgt的SRTM3數(shù)據(jù),采樣的起始緯度為北緯27°,之后每個采樣點的緯度增加1/1200°,終止緯度為北緯28°;采樣的起始經(jīng)度為東經(jīng)86°,之后每個采樣點的經(jīng)度增加1/1200°,終止經(jīng)度為東經(jīng)87°。SRTM數(shù)據(jù)格式如圖1所示。

圖1 SRTM數(shù)據(jù)格式示例

2 雷達可視域圖繪制流程

本文所述的基于SRTM的雷達可視域圖繪制流程如圖2所示。

圖2 基于SRTM的雷達可視域圖繪制流程

從圖2中可以看出,在設(shè)置好輸入?yún)?shù)以后就可以從SRTM文件中獲取雷達位置所在區(qū)域的高程數(shù)據(jù),根據(jù)前述可知獲取到的SRTM數(shù)據(jù)范圍是1°×1°的經(jīng)緯度方塊。然后從起始方位角開始以一定的步長(步長一般是1°)掃描整個扇區(qū),循環(huán)進行某一個方位上的可視性計算,循環(huán)結(jié)束的條件是到達扇區(qū)的邊界。在得到所有方位上的可視性后,繪制可視域圖。

各步驟的詳細設(shè)計在下文中進行描述。

2.1 設(shè)置輸入?yún)?shù)

輸入?yún)?shù)包含雷達參數(shù)和通視參數(shù)。雷達參數(shù)包含雷達所在位置的經(jīng)度L0、緯度B0和高程H0,以及雷達天線的高度ha;通視參數(shù)包含可視半徑Rmax、起始方位角θ0、掃描的扇區(qū)寬度θl、掃描的步長θs、選定點相對雷達的俯仰角φ、目標相對雷達的高度hr。

2.2 獲取雷達位置所在區(qū)域的SRTM高程數(shù)據(jù)

根據(jù)雷達所處位置的經(jīng)緯度可以解析出雷達位置所在區(qū)域的SRTM文件名,再根據(jù)文件名從SRTM地圖數(shù)據(jù)文件夾中讀取相應的SRTM數(shù)據(jù),將獲取到的數(shù)據(jù)的經(jīng)度和緯度分別保存在一維向量中,高程值保存在二維矩陣中。

以我國的地形數(shù)據(jù)為例(北緯、東經(jīng)),經(jīng)度L0、緯度B0和文件名filename的關(guān)系為

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

(1)

其中,N代表北緯,E代表東經(jīng),floor()函數(shù)是向下取整函數(shù),.hgt是文件后綴名。

例如位于東經(jīng)108.332571°、北緯39.823789°位置的雷達,其SRTM數(shù)據(jù)所對應的文件名為N39E108.hgt。

高程數(shù)據(jù)在SRTM文件中的存儲格式是16bit的signed integer,每兩個字節(jié)以空格分隔,分別為高8位數(shù)據(jù)和低8位數(shù)據(jù)。因此,需要提取第一個字節(jié)乘以256,再與第二個字節(jié)相加,得到原始高程數(shù)據(jù)。

以上原始高程數(shù)據(jù)中還存在一些特殊數(shù)據(jù),例如32768、65535等,屬于異常值,對于這些數(shù)據(jù)需要進行異常值修正,修正后的矩陣才能用于后續(xù)計算。本文采用臨近點的高程值來修正異常值,具體做法如下:

1)找出二維矩陣中的異常點坐標;

2)對于非第一列數(shù)據(jù)的異常值取其左側(cè)相鄰點賦值;

3)對于第一列數(shù)據(jù)的異常值,若非第一行則取其上一行相鄰點賦值;

4)對于第一列第一行數(shù)據(jù)的異常值,取其相鄰三個數(shù)據(jù)的平均值進行賦值。

2.3 某一個方位上的可視性計算

通常研究可視域方位角并非只有一個角度,而是一段扇形區(qū)域。針對固定的可視半徑,可以將扇形區(qū)域的問題簡化為某一個方位上的問題,只不過整個流程需要做多次循環(huán)計算。如圖3所示為某一個方位上的可視性計算流程。設(shè)循環(huán)計數(shù)值為i,則該方位角θi的大小為式(2)所示。

圖3 某一個方位角上的可視性計算流程

(2)

其中,各符號的含義參見2.1節(jié)。

2.3.1 根據(jù)選定點的RAE坐標計算LBH坐標

將每個方位上的選定點設(shè)為可視半徑處的最遠點,其RAE坐標指的是選定點相對雷達的距離、方位角和俯仰角坐標,根據(jù)2.1節(jié)和式(2)可知RAE坐標為(Rmax,θ,φ)。

LBH坐標指的是選定點的經(jīng)度、緯度和高程坐標,其坐標轉(zhuǎn)換過程如下:

1)由選定點的RAE坐標(Rmax,θ,φ)計算選定點的雷達站心直角坐標(rx,ry,rz)

(3)

2)由雷達的經(jīng)緯度坐標(L0,B0,H0)和選定點的雷達站心直角坐標(rx,ry,rz),計算選定點的地心直角坐標(X,Y,Z)

(4)

(5)

(6)

3)采用迭代法,由選定點的地心直角坐標(X,Y,Z)計算選定點在WGS84坐標系下的經(jīng)緯度坐標LBH

(7)

(8)

(9)

2.3.2 拼接SRTM數(shù)據(jù)

由2.3.1節(jié)中的計算步驟可以得到選定點的LBH坐標,根據(jù)L和B的值判斷選定點是否超出當前SRTM數(shù)據(jù)的范圍。若超出則重復2.2節(jié)中的計算過程,獲取雷達所在位置到選定點之間所有區(qū)域的SRTM高程數(shù)據(jù),并與原SRTM數(shù)據(jù)進行矩陣拼接;否則不進行數(shù)據(jù)拼接。

2.3.3 構(gòu)建高程線方程

高程線指的是從雷達所在位置到選定點之間的連線,連線上各采樣點處的高程值必須能夠獲取。本文采用等間隔采樣法,每隔相同的距離采樣一個點,可以得到采樣點數(shù)目num為

(10)

其中,Rinter表示采樣間隔,本文使用的間隔為100m。采樣點數(shù)目與方位無關(guān),這樣可以保證各個方位上的采樣點數(shù)目一致。接下來的關(guān)鍵問題是如何獲取高程線上各采樣點處的高程值。

根據(jù)第1節(jié)可知,SRTM數(shù)據(jù)也是由離散的采樣點構(gòu)成,單個文件所在區(qū)域的采樣點數(shù)目是1201×1201(SRTM3格式)。又根據(jù)2.3.2節(jié),即使當前SRTM數(shù)據(jù)的范圍已經(jīng)囊括整條高程線的數(shù)據(jù)范圍,但是由于兩者采樣間隔的不同并不能保證經(jīng)緯度重合,即高程線上采樣點的經(jīng)緯度不一定在SRTM區(qū)域中有對應的離散值,所以需要求高程線上采樣點高程的近似值。

本文采用最近距離法獲取高程線上各采樣點高程的近似值:假設(shè)要獲取高程線上采樣點j的高程值,則在SRTM高程數(shù)據(jù)矩陣中,選取與點j經(jīng)緯度距離最近的點的高程值作為點j的高程值。根據(jù)此方法可以近似獲得高程線上所有采樣點的高程值。

注意,此時獲得的高程值并沒有考慮地球曲率對距離帶來的影響,所以還需要借助地球曲率修正算法,對所有高程值進行曲率修正。

(11)

其中:Hcorr為曲率修正后的高程值;Hori為修正前的原始高程值;R0為地球的等效半徑;Rj為采樣點j到雷達的水平距離。

2.3.4 計算高程線上各個采樣點的最大遮蔽角

遮蔽角的含義是雷達天線頂端與障礙物頂端的連線和地平線之間的夾角[5],如圖4所示,兩條虛線之間的夾角αj就是一個遮蔽角。當目標與天線頂端的連線和地平線之間的夾角(也叫目標的俯仰角)小于該遮蔽角時目標就會被遮擋,所以有遮蔽的含義。

圖4 采樣點j處的遮蔽角示意圖

根據(jù)圖4所示,可以得出高程線上第j個采樣點的遮蔽角αj的計算公式為

(12)

其中:Hj表示該采樣點處曲率修正后的高程值;H0表示雷達的高程;ha表示雷達天線的高度;Rj表示該采樣點到雷達的水平距離。

而某個采樣點位置的最大遮蔽角,指的是該采樣點以及其前面所有采樣點的遮蔽角的最大值。如圖5所示,αA是采樣點A處的遮蔽角,αB是采樣點B處的遮蔽角,雖然HB>HA,但是有αB<αA?？梢钥吹皆谶@種情況下,障礙物B已經(jīng)被障礙物A遮擋,所以B點處的最大遮蔽角是更大的αA。

圖5 最大遮蔽角示意圖

根據(jù)以上描述可以獲得最大遮蔽角的計算方法:使用一個全局變量保存最大遮蔽角記為αmax,然后循環(huán)遍歷每個采樣點并計算其遮蔽角αj,將每次計算出的αj與αmax作比較并更新αmax為較大值,則該采樣點的最大遮蔽角αMj就是此時的αmax。按照此方法,圖5中A點的最大遮蔽角為αA,B點的最大遮蔽角也是αA。

可以看到隨著采樣點的增多最大遮蔽角在不斷更新為更大的遮蔽角,這樣就可以保證每個采樣點處的最大遮蔽角都是其前面所有采樣點遮蔽角的最大值。最終的計算結(jié)果為一個向量,保存著每個采樣點處的最大遮蔽角。

2.3.5 計算目標高度下各個采樣點處的可視性

有了每個采樣點的最大遮蔽角,就可以計算出目標高度下每個采樣點處的可視性。根據(jù)2.3.4節(jié)中遮蔽角的含義可得,目標高度下在某一采樣點處,當雷達與目標之間連線的俯仰角大于該采樣點的最大遮蔽角時,該采樣點處為可視,反之不可視。第j個采樣點處的可視性vj表示為

(13)

其中:αTj為目標在相對雷達hr高度下該采樣點處的俯仰角;αMj為該采樣點處的最大遮蔽角。根據(jù)幾何關(guān)系有式(14)為

(14)

其中:H0表示雷達的高程;ha表示雷達天線的高度;Rj表示該采樣點到雷達的水平距離。示意圖如圖6所示。

圖6 目標相對雷達hr高度下采樣點j處的可視性

2.4 繪制可視域圖

在扇區(qū)內(nèi)各個方位上重復2.3節(jié)的計算過程,就能得到hr高度下各方位上高程線各個采樣點處的可視性,據(jù)此就可以繪制出可視域圖。

本文使用Matlab的polarplot()函數(shù)繪制可視域圖。polarplot()函數(shù)的作用是繪制極坐標圖,傳入的第一個參數(shù)是角度向量,第二個參數(shù)是待繪制參數(shù)的向量。

繪制方法如下:

2)遍歷各個采樣點,采樣點范圍從1到num,其中num可由式(10)得到;

對于每個采樣點j,設(shè)其可視距離向量為rangeVector,向量長度與angleVector相同。求rangeVector的步驟如下:

①遍歷angleVector中的各個方位角θi

ⅰ.取方位角θi上,高程線上采樣點j處的可視性vij;

ⅱ.取方位角θi上,高程線上采樣點j距雷達的水平距離Rij;

ⅲ.可視距離向量的第i個元素rangeVectori=vijRij,將rangeVectori添加到rangeVector中。

②遍歷完angleVector后就會得到采樣點j處的可視距離向量rangeVector,使用polarplot(angleVector,rangeVector)函數(shù)繪制采樣點j處的單個圓;

3)遍歷完所有采樣點后,此時就繪制了num個圓?？梢曈驁D繪制結(jié)束。

圖7為可視域圖繪制原理示意圖。雷達位于點O,可視半徑設(shè)為600m,在5個方位上全部都每隔100m采樣6個點,形成射線OA到OE。

圖7 可視域圖繪制原理示意圖

步驟1)中的angleVector的取值為射線OA到OE的5個方位角的大小。

步驟2)對6個采樣點進行遍歷,每個采樣點都有自己的rangeVector,其長度都為5(圖中離散的弧線)。假設(shè)圖中所有的點都可視,則采樣點1的rangeVector值為{100,100,100,100,100},采樣點2的rangeVector值為{200,200,200,200,200},其他采樣點類似。對于每個采樣點都繪制出離散的圓弧。

步驟3)可以繪制出整個圖7。

3 仿真結(jié)果與分析

本文使用Matlab繪制出可視域圖之后,還與一款專業(yè)的地圖繪制軟件GlobalMapper做了比較。

3.1 GlobalMapper軟件簡介

GlobalMapper軟件是GlobalMapper軟件公司開發(fā)的一款地圖繪制軟件,支持多種數(shù)據(jù)源的輸入和讀取(包含SRTM),可以處理柵格和矢量數(shù)據(jù)。軟件預置了非常多的常用坐標系和轉(zhuǎn)換參數(shù),提供了各種坐標轉(zhuǎn)換方法,支持數(shù)百種地圖投影。此外還擁有強大的地圖編輯、轉(zhuǎn)換、打印等功能[6]。

本文所述的方法需要參考GlobalMapper軟件中的視圖分析工具,其界面如圖8所示。

圖8中:“發(fā)射點高度(地面之上)”指的是雷達天線的高度ha;“接收點高度(地面之上)”指的是目標相對雷達的高度hr;“視角”參數(shù)中的“起始角度”指的是起始方位角θ0,“掃描的角度”指的是掃描的扇區(qū)寬度θl,掃描的步長θs默認為1°;選定點相對雷達的俯仰角φ默認為0;“可視半徑”含義相同指的是Rmax;而雷達所在位置的經(jīng)緯高參數(shù)L0、B0和H0在鼠標點擊地圖時可以由軟件自動獲取。

3.2 仿真結(jié)果對比分析

本文通過兩個典型的示例,分別在雷達仰視和下視的情況下對比了Matlab的仿真結(jié)果和GlobalMapper軟件的視圖分析結(jié)果。

1)雷達仰視的情況

參數(shù)如下:

L0=78.943321°
B0=13.165818°

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

仿真結(jié)果如圖9所示。

圖9中的右側(cè)部分為GlobalMapper軟件繪制出的可視域圖,灰色陰影區(qū)域表示可視范圍,其半徑表示可視距離?？梢钥吹阶笥覂蓚€結(jié)果的可視域輪廓基本一致,在方位90°～120°左右可視距離能夠達到50km,視線良好,其余方位的可視距離不夠50km。

2)雷達下視的情況

參數(shù)如下:

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

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

仿真結(jié)果如圖10所示。

圖10 雷達下視情況的仿真結(jié)果

可以看到此參數(shù)下可視域不連續(xù)。這是因為假設(shè)目標能夠低空貼地飛行,其障礙物遮擋情況將變得復雜。

綜上,本文所述的可視域圖繪制方法能夠很好地貼合于專業(yè)繪圖軟件GlobalMapper的繪制結(jié)果,證明了該算法的可實現(xiàn)性和有效性。

4 結(jié)束語

在雷達選址時可能遇到地物遮擋,因此需要繪制雷達站周圍一定范圍內(nèi)的可視域圖,針對于此,本文首先通過解析輸入?yún)?shù)、獲取SRTM高程數(shù)據(jù)和坐標轉(zhuǎn)換建立各個方位上的高程線方程;然后計算出高程線上各個采樣點處的最大遮蔽角;再計算出目標高度下高程線上各個采樣點處的可視性;最后,使用Matlab軟件繪制出了可視域圖。本文還將Matlab仿真結(jié)果與地圖繪制軟件GlobalMapper進行了對比,驗證了該方法的可行性和有效性。