何 強(qiáng)，宋 錚，馬斯敏

(1.電子工程學(xué)院,安徽 合肥 230037；2.西安通信學(xué)院，陜西 西安 710106)

基于MATLAB提取數(shù)字地圖的電波傳播問(wèn)題研究

何 強(qiáng)1，宋 錚1，馬斯敏2

(1.電子工程學(xué)院,安徽 合肥 230037；2.西安通信學(xué)院，陜西 西安 710106)

在拋物方程的基礎(chǔ)上引入了地形高程剖面，預(yù)測(cè)了復(fù)雜地形條件下的電波傳播特性。首先，基于MATLAB提取調(diào)用了數(shù)字地圖，給出了其主要程序流程圖以及提取調(diào)用的方法，之后采用插值方法增加了地形模型與實(shí)際地形剖面的契合度；其次，使用MATLAB的GUI功能設(shè)置了一個(gè)高程剖面系統(tǒng)，該系統(tǒng)能夠提取地球上任意兩點(diǎn)之間的地形高程剖面；最后，通過(guò)將高程剖面系統(tǒng)提取的地形高程剖面加入到拋物方程模型中進(jìn)行計(jì)算，分析了在復(fù)雜地形環(huán)境下電波傳播的特性。

數(shù)字地圖；拋物方程；地形剖面；傳輸損耗

0 引 言

數(shù)字地圖是高科技迅猛發(fā)展的一個(gè)產(chǎn)物，隨著它的不斷完善，現(xiàn)今的數(shù)字地圖幾乎能夠涵蓋所有地球上地形的特征信息。數(shù)字地圖實(shí)際上就是將地球上的地理信息數(shù)字化，可以稱之為數(shù)字地球。數(shù)字地圖是以數(shù)字的形式記錄和存儲(chǔ)的地圖，通過(guò)它能夠十分方便地提取地形的地理位置信息。由于數(shù)字地圖的時(shí)效性，使得其無(wú)論是在民用領(lǐng)域還是在軍用領(lǐng)域都有十分重要的使用價(jià)值。

因?yàn)閿?shù)字地圖所能夠帶來(lái)的巨大效益，所以其提取應(yīng)用的方法也廣受人們關(guān)注[1]。本文提出了基于MATLAB提取數(shù)字地圖的方法，通過(guò)MATLAB編程調(diào)用數(shù)字地圖具有良好的實(shí)用性，工程價(jià)值較高。

預(yù)測(cè)電波在復(fù)雜地形環(huán)境下的傳播特性是研究拋物方程法(簡(jiǎn)稱PE)的一個(gè)重要方向。由于PE的計(jì)算十分復(fù)雜，而MATLAB編程相對(duì)于其它編程軟件更為方便，并且數(shù)字地圖也可以通過(guò)MATLAB[2]來(lái)提取，所以本文選用MATLAB進(jìn)行編程計(jì)算。如果能夠?qū)佄锓匠膛c數(shù)字地圖相結(jié)合，把提取的地形高程剖面加入到拋物方程的計(jì)算當(dāng)中，那么使用拋物方程法計(jì)算的效果更具有實(shí)用性。

本文使用MATLAB提取數(shù)字地圖的地理信息并畫(huà)出相應(yīng)實(shí)際地形的高程剖面圖，然后利用插值方法提高了計(jì)算模型的精度。通過(guò)計(jì)算自由空間中電波傳輸損耗和傳播因子，分析了電波在自由空間中的傳播特性，并將地形高程剖面加入到拋物方程法的計(jì)算當(dāng)中，最后得到了電波傳輸損耗、傳播因子以及它們的偽彩圖。

1 拋物方程法的介紹

拋物方程法[3]是一種全波分析方法，它充分考慮了初始場(chǎng)、繞射效應(yīng)、折射效應(yīng)、上邊界條件、阻抗邊界條件以及不規(guī)則地形等因素對(duì)電波傳播的影響，因而能夠被廣泛應(yīng)用于分析電波傳播的問(wèn)題當(dāng)中。其中，初始場(chǎng)可以通過(guò)對(duì)方向圖函數(shù)進(jìn)行傅里葉變換得到，上邊界條件為窗函數(shù)，阻抗邊界條件為L(zhǎng)enontovich邊界條件。在胡繪斌的論文[4]中詳細(xì)推導(dǎo)了直角坐標(biāo)系、柱坐標(biāo)系以及球坐標(biāo)系下的拋物方程，本文用到的是球坐標(biāo)系下的拋物方程。在球坐標(biāo)系下討論電波的傳播過(guò)程中，設(shè)時(shí)諧場(chǎng)的時(shí)變因子為e-jωt，忽略波的后向傳播，則可以得到寬角拋物方程(WAPE)(在胡繪斌的論文中給出了推導(dǎo)過(guò)程，由于篇幅有限在這里不作介紹)：

(1)

其中坐標(biāo)轉(zhuǎn)換的具體表達(dá)式：

(2)

式中：a為地球半徑；h為距地表高度。

根據(jù)參考文獻(xiàn)[5]，通過(guò)分布傅里葉變換(SSFT)算法推導(dǎo)上式得到：

(3)

從解的形式中可以知道SSFT是一種步進(jìn)迭代算法，該算法受波長(zhǎng)限制小，穩(wěn)定度高，求解速度快。而分布傅里葉變換的核心是將折射項(xiàng)和繞射項(xiàng)相分開(kāi)，可以單獨(dú)對(duì)繞射項(xiàng)進(jìn)行傅里葉變換運(yùn)算。

(4)

LF=20lg|u(x,z)|+10lgx+10lg10λ

(5)

2 提取高程數(shù)據(jù)

2.1 地形數(shù)據(jù)

數(shù)字高程模型(DEM)是使用一組有序數(shù)列形式表示地面高程的一種實(shí)體地面模型，由于數(shù)字地圖能夠?qū)⑸酱ā⒑恿?、海面及陸地的屬性和地理位置信息表現(xiàn)出來(lái)，因此其能夠在計(jì)算機(jī)中真實(shí)地反映地形情況。

本文中，數(shù)字地圖的地形信息數(shù)據(jù)存儲(chǔ)的格式為Grid型DEM-DTED，其高程數(shù)據(jù)是以網(wǎng)格的數(shù)據(jù)進(jìn)行存儲(chǔ)的。DTED存儲(chǔ)格式是以左下角為起點(diǎn)，每一列為一個(gè)記錄單元，并且經(jīng)度、緯度以遞增的方式存儲(chǔ)高程值。在數(shù)字地圖準(zhǔn)備好了以后，MATLAB中DTED函數(shù)可以載入始末2點(diǎn)內(nèi)的經(jīng)度和緯度高程信息，并以二維數(shù)組形式儲(chǔ)存。圖1就是DTED的一個(gè)小文件，表達(dá)了經(jīng)度從東經(jīng)116°到117°，緯度從北緯30°到31°區(qū)域內(nèi)的高程數(shù)據(jù)。

2.2 插值算法

通過(guò)數(shù)字地圖，只需要知道任意2點(diǎn)的經(jīng)緯度信息，就可以提取它們之間的地形高程剖面。圖2是發(fā)射點(diǎn)(116.477 9°E，30.569 3°N)和接收點(diǎn)(116.437 1°E，30.610 2°N)之間的高程剖面，其高程剖面是由一個(gè)個(gè)離散點(diǎn)組成的，而特定點(diǎn)位的坐標(biāo)信息是唯一的。在實(shí)際環(huán)境條件下，由于所考慮的地形是不固定的，所以在數(shù)字地圖中提取的離散采樣點(diǎn)的位置也是隨機(jī)的。為了使PE中的每個(gè)步進(jìn)都有相對(duì)應(yīng)的離散點(diǎn)，也為了提高計(jì)算的精度，本文中引入內(nèi)插法進(jìn)行計(jì)算。

數(shù)字地圖的數(shù)據(jù)信息是由很多個(gè)離散點(diǎn)組成的，相鄰2點(diǎn)之間數(shù)據(jù)信息相近，關(guān)聯(lián)性較大，而內(nèi)插方法正是通過(guò)在相鄰2點(diǎn)之間做連線，再根據(jù)這2點(diǎn)的高程數(shù)據(jù)信息可以求出連線上任意一點(diǎn)的高程數(shù)據(jù)。圖3是它的原理分析圖。

A(x1,y1),B(x2,y2),C(x3,y3)是數(shù)字地圖中的存儲(chǔ)數(shù)據(jù)，A和B兩點(diǎn)之間連線上任意一點(diǎn)的坐標(biāo)為D(x4,y4)，有：

(6)

變換上式可以得到：

(7)

所以只需要知道所求高程點(diǎn)的橫坐標(biāo)x4，就可以得到其縱坐標(biāo)及高程信息。當(dāng)實(shí)際計(jì)算時(shí)，x4表示的是某個(gè)步進(jìn)的距離，它可以從MATLAB中獲得。

2.3 地形剖面的提取

2.3.1 程序流程圖

程序流程圖如圖4所示。

2.3.2 高程剖面系統(tǒng)

使用MATLAB中GUI[6]功能按照程序流程圖4設(shè)計(jì)一個(gè)系統(tǒng)，如圖5所示。該系統(tǒng)中有高程剖面[5]和傳輸路徑高程2個(gè)模塊。在高程剖面模塊中有發(fā)射站和接收站位置按鈕、收發(fā)站的經(jīng)緯度信息以及傳播距離(不輸入，會(huì)自動(dòng)算出)；傳輸路徑高程模塊在高程剖面模塊輸入地理位置數(shù)據(jù)后進(jìn)行計(jì)算就可以自動(dòng)輸出。圖6是輸入數(shù)據(jù)后的界面。

該系統(tǒng)可以提取任意2點(diǎn)之間的地形剖面信息，輸出圖形的橫坐標(biāo)為距離(km)，縱坐標(biāo)為高度(m)。該系統(tǒng)操作十分簡(jiǎn)單，只需要輸入發(fā)射站和接收站的經(jīng)緯度信息，然后點(diǎn)擊發(fā)射站和接收站位置按鈕則傳播距離以及傳輸路徑高程都會(huì)自動(dòng)算出。

3 PE與數(shù)字地圖相結(jié)合分析電波傳播問(wèn)題

3.1 地形模型的選取

根據(jù)實(shí)際地形的不同，許多學(xué)者提出了不同的地形模型，其中有地形遮蔽模型[7]、連續(xù)移位變換模型[8]、離散移位變換模型、寬角移位變換模型[9]等。而地形遮蔽模型最適合應(yīng)用于數(shù)字地圖的地形剖面中，故本文選用地形遮蔽模型。

該模型是由S.Ayasli和Newkirk等人提出的，其基本原理就是將復(fù)雜地形以階梯化的形式來(lái)表示，如圖7所示。

在分析電波傳播時(shí)，只考慮垂直地形S1面上的場(chǎng)，并將S2面上的場(chǎng)設(shè)置為零，模型的地表下場(chǎng)都設(shè)置為零。采用地形遮蔽模型時(shí)，可以通過(guò)減小步進(jìn)長(zhǎng)度以提高其計(jì)算的精度。

圖8是MATLAB提取數(shù)字地圖中發(fā)射點(diǎn)(116.257 5°E，30.332 0°N)和接收點(diǎn)(116.555 9°E，30.185 8°N)之間的地形高程剖面。圖9是通過(guò)地形遮蔽模型方式處理后的剖面圖(設(shè)定的計(jì)算距離為31 km)，其中每個(gè)臺(tái)階的長(zhǎng)度為1個(gè)步長(zhǎng)。通過(guò)圖8可以看出，該地形起伏變化較大，如果單使用地形遮蔽模型對(duì)該地形剖面進(jìn)行計(jì)算的話則精度較低，因此可以減小步進(jìn)長(zhǎng)度和使用內(nèi)插法提高其階梯化的精度，圖9是修正后的地形模型。

由圖9可以看出，通過(guò)減小步進(jìn)長(zhǎng)度和使用內(nèi)插法修正后的地形模型與原地形十分相近，所以計(jì)算的精度會(huì)得到有效提高。

3.2 計(jì)算實(shí)例

在標(biāo)準(zhǔn)大氣條件下，發(fā)射天線架高為hr=120 m，電波的工作頻率f=1 000 MHz，極化方式為垂直極化，地質(zhì)為非常干燥地，εr=3，電導(dǎo)率為σ=0.000 1 S/m。選擇圖8中的發(fā)射點(diǎn)和接收點(diǎn)位置，但是將收發(fā)2點(diǎn)之間的地形設(shè)為光滑球面地，計(jì)算距離為31 km，設(shè)定的水平步長(zhǎng)為100 m，接收天線架高以h1表示。使用離散混合傅里葉變換(DMFT)算法對(duì)拋物方程進(jìn)行計(jì)算，得到的傳播因子隨

高度以及距離的變化曲線如圖10(a)、(b)所示。

圖10(a)為傳播因子隨高度變化圖，從圖中可以看出，隨著接收天線高度的增加，傳播因子的值在零點(diǎn)左右震蕩且震蕩不斷變小，這是因?yàn)殡S著高度的增加地面反射波能量越來(lái)越小，所以其干涉效應(yīng)逐漸減弱；圖10(b)中傳播因子也是圍繞零點(diǎn)震蕩，而近距離的震蕩較遠(yuǎn)距離的震蕩更加密集，這是因?yàn)榫嚯x越近波程差變化越大，距離越遠(yuǎn)波程差變化越小。

3.3 把高程剖面數(shù)據(jù)加入到拋物方程中

仍然選用圖8的經(jīng)緯度數(shù)據(jù)，但是引入實(shí)際的復(fù)雜地形,重新進(jìn)行相應(yīng)的PE分析計(jì)算，得到了反映電波傳播特性的圖11～圖15。

從圖11可知，在山峰后的山腳處傳輸損耗最大，并且山峰后的傳輸損耗明顯大于山峰前的傳輸損耗。從偽彩圖12中也可以看出山峰后的傳輸損耗大于傳播前的傳輸損耗，對(duì)于電波在山峰后的傳播而言，此山峰類似于尖劈。

圖13計(jì)算的是傳播因子隨接收天線高度的變化情況，隨著高度的增加傳播因子逐漸變大。當(dāng)接收天線的高度高于120 m時(shí)，傳播因子在零值處震蕩，說(shuō)明此山峰的尖劈效應(yīng)開(kāi)始顯現(xiàn)。圖14反映了傳播因子隨距離的變化，從圖中可以看出，山峰前的傳播因子仍然體現(xiàn)了直接波和地面反射波的干涉效應(yīng)，但是山峰的存在導(dǎo)致山后的傳播因子明顯變小。圖15是傳播因子的偽彩圖。

4 結(jié)束語(yǔ)

首先，本文提出了使用MATLAB提取調(diào)用數(shù)字地圖地形高程剖面的方法。該方法簡(jiǎn)單實(shí)用，易于操作，具有良好的實(shí)用價(jià)值；然后，通過(guò)調(diào)整步進(jìn)長(zhǎng)度和使用內(nèi)插法提高了地形遮蔽模型與實(shí)際地形的契合度，從而提高了計(jì)算模型的精度；最后，將MATLAB提取的地形高程剖面加入到拋物方程的地形邊界條件當(dāng)中進(jìn)行計(jì)算，分析了復(fù)雜地形環(huán)境下電波傳播的特性。

[1] 白瑞杰，廖成，盛楠，張青洪.基于二維數(shù)字地圖的電波傳播問(wèn)題仿真[J].強(qiáng)激光與粒子束,2014,26(7)：191-194.

[2] 陳天偉.應(yīng)用MATLAB構(gòu)造數(shù)字高程模型工具箱[J].桂林工學(xué)院學(xué)報(bào)，2003,23(4)：453-456.

[3] 馬斯敏，宋錚，馬恒民.基于FEKO/PE混合算法的近地線天線地波場(chǎng)求解[J].安徽.微波學(xué)報(bào),2014,30(6):46-50.

[4] 胡繪斌,毛鈞杰.預(yù)測(cè)復(fù)雜環(huán)境下電波傳播特性的算法研究[D].西安.西安電子科技大學(xué)出版社,1995.

[5] 馬斯敏，宋錚.基于拋物方程的地面電波傳播特性分析[D].合肥：電子工程學(xué)院，2014.

[6] 羅華飛.MATLAB GUI設(shè)計(jì)學(xué)習(xí)手記[M].北京：航空航天大學(xué)出版社，2014.

[7] NEWKIRK M H.Recent advances in the tropospheric electromagnetic parabolic equation routine (TMPER) propagation model [C]//Battlespace Atmosph.Conf,1997:78-85.

[8] BARRIOSS A E.A terrain parabolic equation model for propagation in the troposphere[J].IEEE Transactions on Antennas and Propagation,1994,42(1):90- 98.

[9] DONOHUE D J.KUTTLER J R.Propagation modeling over terrain using the parabolic wave equation[J].IEEE Transactions on Antennas and Propagation,2000,48(2):260-277.

Study of Electric Wave Propagation Problems in Digital Map Extration Based on MATLAB

HE Qiang1，SONG Zheng1，MA Si-min2

(1．Electronic Engineering Institute of PLA，Hefei 230037，China; 2.Xi'an Communication Institute，Xi'an 710106，China)

This paper introduces the terrain elevation profile based on parabolic equation,forecasts the electric wave propagation characteristics under the complex terrain condition.This paper firstly extracts and calls the digital maps based on MATLAB,gives the main program flow chart and the extraction and call method,then uses interpolation to increase the consilience between terrain model and practical terrain profile;secondly uses GUI function of MATLAB to set an elevation profile system which can extract the terrain elevation profile between random two points on the earth;finally adds the terrain elevation profile extracted by elevation profile system into parabolic equation model and performs calculation,analyzes the characteristics of electric wave propagation in complex terrain environment.

digital map;parabolic equation;terrain profile;propagation loss

2016-03-03

TN011

A

CN32-1413(2017)01-0062-06

10.16426/j.cnki.jcdzdk.2017.01.014