一種基于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進行了對比,驗證了該方法的可行性和有效性。