存在障礙物時電波傳播拋物線方程分析及其驗證?

魏喬菲 尹成友 范啟蒙

(電子工程學院,脈沖功率激光國家重點實驗室,合肥 230037)

存在障礙物時電波傳播拋物線方程分析及其驗證?

魏喬菲 尹成友?范啟蒙

(電子工程學院,脈沖功率激光國家重點實驗室,合肥 230037)

(2017年1月14日收到;2017年3月20日收到修改稿)

雙向拋物線方法主要用于起伏地形下電波傳播問題的計算,該算法本身無法處理地面存在障礙物,尤其是真實環境下障礙物與地面為不同媒質的情況.因此本文提出一種用于存在障礙物時電波傳播計算的拋物線方程新算法.該方法采用區域分解,對不同障礙物區域的場值進行分區計算,并對計算結果進行相位修正,從而實現該情況下空間中場值的計算.在此基礎上,使用矩量法來精確驗證拋物線方法的計算精度.通過實例分析,證明了存在障礙物時新算法的精確性,為之后求解真實環境下的電波傳播問題提供了參考.

拋物線方程,電波傳播,區域分解,矩量法

1 引 言

近年來,雙向拋物線方法[1?3]已經逐步成為對流層電波傳播問題研究中最常用的方法.該算法可以計算不規則地形特征和不同電磁參數的地表結構對電波傳播的影響,在不規則地形的處理方面,拋物線方法[4?7]的計算精度受限于起伏地形的傾角[8?10],傾角超過一定范圍后,雙向拋物線方法的計算存在較大誤差.這種下邊界條件的處理方式只是將障礙物當起伏地形處理,沒有考慮障礙物內存在的場值對空間中場值的影響.因此,當地面存在障礙物時,拋物線方法計算存在困難.

在拋物線方法的計算精度驗證方面,通常在不規則地形下使用AREPS軟件[11]或高頻近似法[12?14]來驗證拋物線方法計算的正確性,高頻近似法主要有菲涅耳積分法及幾何光學結合一致性幾何繞射理論法兩種方式.這兩種方法都是基于射線追蹤法來計算空間中場強的分布,因為這兩種高頻近似法在計算過程中都忽略被障礙物反射的場值,所以主要用于前向拋物線方法計算的驗證.與此同時,障礙物情況越復雜,空間中繞射場的場值計算越困難.除此之外,當波源傳播角度較大時,為提高近場區域的驗證精度,Omaki等[15]提出用時域有限差分方法[16]驗證拋物線方法計算結果的精確性.

結合上述情況,本文在原有起伏地形下雙向拋物線方法計算方式的基礎上,提出了一種新的拋物線方程算法,該算法可以較為精確地計算地面存在障礙物時空間中的場值大小.當地面存在障礙物時,根據區域分解原理[17]將障礙物處的場值計算分為兩部分,并分別增加一段重疊區域,以滿足障礙物上邊界處場值計算的連續性.針對不同區域的場值計算,使用不同的拋物線計算方法.除此之外,在場值的迭代計算中對總場計算時的相位差異和障礙物傳播中造成的相位差異兩方面都進行了相位修正.在新算法的驗證方面,本文用矩量法[18]來驗證近距離下雙向拋物線方法求解電波傳播問題的計算精度.最后通過一系列實例分析驗證了新算法計算的精確性.

2 障礙物下拋物線方程的處理方法

設電磁場的時諧因子為exp(jωt),本文主要計算討論的是地面存在障礙物時雙向拋物線方法的具體求解,假設地面為理想導體,障礙物近似如圖1所示.

圖1 障礙物(矩形)模型Fig.1.The model of obstacle(rectangle).

圖1中ht,hr,h分別代表源的發射高度、觀察點的接收高度和障礙物的高度.在此基礎上,本節主要就雙向拋物線方法計算中的波源設置、基本原理以及對計算結果的相位修正三個方面進行分析說明.

2.1 初始場設置的等效源模型

常用的拋物線方法的波源設置是利用遠場與近場之間的關系反推初始場表達式,這種初始場u(0,z)的設置方法存在波束較窄或天線架高不能滿足u(0,z)=0(z<0)時,初始場會延伸到地面以下的問題,并且這種方法沒有建立起初始場與天線電流分布之間的直接聯系.結合上述情況,本文使用等效源模型來逼近初始場設置[19,20].

已知拋物線方程由初始場推導出的自由空間遠區場u(x,z)的解析形式如下[21]:

其中H(2)1表示一階第二類漢克爾函數,ρ為源點到場點的距離,k為自由空間的波數.現假設有一個在y方向無限長的有具體分布數值橫向磁流帶,電流方向為z方向,

則該電流產生的電場為

由(3)和(1)式,可以看出兩者具有相似的形式.根據對偶定理,同理可推導出等效電流產生的磁場表達式,由此推出初始場與等效源的關系如下:

在(4)式的基礎上,可以通過等效源直接逼近拋物線方法的初始場設置,這種初始場設置建立了其與天線電流分布之間的直接聯系,不僅解決了之前方法存在的問題,而且為后面利用矩量法分析驗證拋物線方程的計算精度提供了統一的入射源環境.

2.2 雙向拋物線方法的基本原理

在球面周向均勻起伏的地形下,通常使用雙向拋物線方法來求電波傳播問題.已知根據平坦地球模型,可得電磁場滿足二維標量波動方程:

其中U表示電磁場量Ey,Hy;m代表修正后的大氣折射率.為解調(5)式中的相位快變部分,在波動方程的基礎上引入軸向衰減因子如下[1?3]:

其中UF,UB代表前向場和后向場.將(6)式分別代入(5)式,并使用Feit-Fleck寬角近似處理,得到雙向波動方程為

當地面為理想導體時,通過分步傅里葉方法求得雙向拋物線方法的步進公式如下:

其中F和F?1分別代表傅里葉變換和逆變換.可以看出在忽略傳播方向的前提下,推導出的步進公式是完全一致的.為進一步分析雙向拋物方程在不規則地形下的具體實現方法,以圖1中的障礙物模型為例,文獻[1]中雙向拋物線方法的計算方式如圖2所示.

圖2 單矩形障礙物下空間場值的處理 (a)前向傳播場;(b)后向傳播場Fig.2.The processing of spatial fi eld under single rectangular obstacle:(a)The forward fi eld;(b)the backward fi eld.

以圖2中對單矩形障礙物的處理方式,可得空間中場值設定如下:

由(9)式可以看出,當迭代計算到障礙物時,障礙物之下的場值被強行置零,只計算上空間的場值.文獻[1]中的這種計算方式將障礙物視為起伏地形,沒有考慮障礙物內透射場的存在,也沒有考慮不同極化下障礙物邊界處場值的處理.

2.3 雙向拋物線方法的相位補償

在(9)式中場值分類的基礎上,雙向拋物線方法求解空間中總場如(6)式所示,是將前向和后向波統一到一個坐標原點得出的,在疊加計算時忽略了前向場與后向場相位不一致的問題.從圖2(b)中對后向傳播場的處理中可以看出后向波源是在x1處產生的,根據障礙物反射面x1處的邊界條件可得

其中Γ為自由空間向障礙物垂直入射下的反射系數.在迭代求解空間中任意一點x處的總場時,可以看出前向波是從原點出發,后向波是從x1出發,此時x點的場強應該為

對比(11)式與(6)式得,在求解總場時,需要先進行相位補償將前向波與后向波的相位統一后疊加.將(11)式代入(6)式得

通過(12)式可以修正雙向拋物線方法總場計算時的相位差異,但是在此之前的計算處理都是將障礙物處理為起伏地形,沒有從根本上處理存在障礙物時空間中的場值該如何計算的問題.因此本文提出一種新的拋物線方程計算方法,將障礙物處場值的計算進行單獨處理.

3 存在障礙物時拋物線方程的處理方法

當地面存在障礙物時,經典的拋物線方法是將障礙物與地面視為整體,將起伏地形之下的場值全部置零,并根據地面媒質參數的不同,使用不同的分步傅里葉變換法.這種計算方法存在以下兩方面的問題:第一方面是因為該方法的計算精度受限于障礙物的傾角大小;第二方面是沒有考慮障礙物內部場值對空間中場值分布的影響.本文提出了一種區域分解算法.

3.1 區域分解算法

以圖2中對單矩形障礙物的處理為例,經典拋物線方法的計算區域是從地面一直到吸收層.本文在此基礎上,對障礙物區域場值的求解進行了分區處理,當步進求解到障礙物處,場值的計算分為兩部分,第一部分是障礙物之上的場值計算,第二部分是障礙物之下的場值計算.根據區域的不同,采用不同的離散傅里葉變換方法,并對障礙物上邊界處的場值計算加以修正.根據區域分解原理,對場值計算處理如下.

由圖3可以看出,障礙物區域的場值計算主要分為三部分.

1)障礙物上半空間的場值計算.由圖3可以看出,計算障礙物上半空間的場值時,該區域中的下邊界從障礙物的上邊界B點往障礙物之下延伸了一段,這種延伸方式是為了避免對區域突然截斷時造成的計算誤差.對重疊區延伸長度的選擇與計算精度有關,通常為了計算方便,直接取障礙物高度的一半.在此基礎上,此時上空間的下邊界不再是理想導體而是均勻介質,所以上半空間的迭代計算相當于阻抗邊界下的混合傅里葉求解.對下邊界的處理等效為Leontovich阻抗邊界條件[21]

根據(13)式,通過離散混合傅里葉變換法就可以求解障礙物之上包括重疊區域的場值.

圖3 障礙物處場值計算的區域分解Fig.3.Domain decomposition of fi eld calculation in obstacle zone.

2)障礙物內部空間的場值計算.此時介質障礙物內場值計算的上邊界也發生了改變,為了提高計算精度,直接將障礙物內部的場值計算擴充到整個空間中,因為地面仍是導體,所以下邊界的處理使用導體邊界條件,傅里葉變換采用正弦變換計算,最后出障礙物恢復空間場值時只提取障礙物高度h以下區域的場值.

3)求解完障礙物上下空間的場值后,如圖3所示,在求解障礙物u(x±?x,z)時,需要對介質障礙物上邊界C處的場值進行修正.該點處場值計算如下:

(14)式是將連接的邊界點處的場值使用兩種區域算法的平均值代替,這樣確保邊界切向場連續.

綜上,在整個空間的迭代計算中遇到障礙物處就將計算區域按圖3所示分成兩部分計算,其余平坦區域的場值計算方法不變.

3.2 介質中拋物線方程的處理及相位補償

在圖3的基礎上,障礙物之下的場值計算等效于介質體內的電波傳播問題的求解.因為地面仍為理想導體,所以在障礙物內的迭代計算與自由空間的迭代計算相比區別在于空間中的折射率發生了改變,等效為(6)式的解調因子的改變:

在(16)式步進求解的基礎上,還需要考慮進出障礙物時電磁波的透射作用,場值修正為

當障礙物為理想導體時,因為障礙物內不存在場值,所以只需要考慮上半空間的場值求解.又因為此時上半空間的下邊界為導體,所以根據極化模型的不同,混合傅里葉變換法也簡化為正弦或余弦傅里葉變換法,并在恢復場值時,直接將下半空間的場值賦零.

4 雙向拋物線方程矩量法驗證分析

為了驗證所提算法的精度,本文提出在近距離下用矩量法精確求解空間中的場強分布,從而驗證存在障礙物的情況下新算法計算結果的精確性.因為障礙物為理想導體,也等效為一種均勻介質,所以下面主要考慮矩形障礙物為均勻介質的情況,并分別從兩種極化角度出發分析討論矩量法的處理過程.

如圖1所示的障礙物,將地面影響等效到格林函數G中,只對障礙物采用矩量法分析進行處理.由于障礙物為均勻介質,所以可以用Poggio-Miller-Chang-Harrington-Wu(PMCHW)方程和體等效兩種方法來求解空間中的場強分布,本文選擇PMCHW方程組來聯立求解上述問題:

其中Hi,Ei,J,M分別代表入射磁場、入射電場以及障礙物表面的等效電流與等效磁流;z1,z2分別指自由空間和均勻介質障礙物內的波阻抗.L(X),K(X)都是簡化表達的算子,定義如下:

聯立(19)式中的算子,通過(18)式中等效電磁流與入射場之間的關系可以實現任意極化下空間中的場值求解.因為存在導體地面,考慮地面鏡像作用,對于電磁流源,(19)式中的二維格林函數為

其中ρi,ρ′i代表源點及其鏡像點的位置矢量,G0表示零階第二類漢克爾函數,鏡像格林函數的正負取決于等效電磁流的鏡像原理.因為在求解過程中存在二階導數的計算,所以對于縱向電磁流采用點匹配脈沖基求解,橫向電磁流基函數和檢驗源均采用三角基函數矩量法求解.

在TE波(垂直極化)照射下,障礙物在空間中產生的散射磁場只有y分量.因此使用二維等效電流來逼近拋物線方法的初始場,在垂直極化下,障礙物表面會產生環向的等效電流,流向同圖1,障礙物表面的等效磁流為y向,進而求得空間中總場如下:

在TM波(水平極化)照射下,空間中只有y向電場,因此使用二維的等效磁流來等效初始場,在水平極化下,障礙物表面的等效電流為y向,等效磁流變成環向,積分方程中的各項分量的求解可以通過對偶原理從(21)式直接求得,惟一的區別在于考慮鏡像作用時,鏡像電流的等效與鏡像磁流的等效相反.同理求得空間中總電場如下:

5 計算實例和結果分析

本文分析的問題基于圖1所示的矩形障礙物模型,用矩量法計算近距離下電波傳播的具體場值,從而驗證在存在障礙物情況下新算法計算結果的精確性,并在障礙物為理想導體的情況下將新算法與文獻[2]的PEtool拋物線方法進行對比分析.在下述實例分析中,矩量法的剖分間隔為0.5 m,新算法的步進間隔是1 m,并且所有實例都是在i5處理器、32位操作系統、4G內存的計算機設備環境下運行仿真的.

實例一:發射源工作頻率30 MHz,高度為200 m,其發射波束仰角0?,波束寬度20?,水平極化,輻射源等效輻射功率為W,Z0為自由空間波阻抗.障礙物為雙矩形理想導體,矩形模型形狀一致,第一個矩形模型位于距離發射源750 m處,第二個矩形模型位于距離發射源900 m處,障礙物模型寬50 m,高度100 m.計算得到在接收高度為80 m處場強隨水平距離變化的場強振幅分布如圖4(a)所示,圖4(b)為圖4(a)的局部展開結果,其中MoM代表矩量法計算結果,PE2代表本文提出的存在障礙物情況下基于雙向拋物線方程的新算法的計算結果,PEtool代表文獻[2]算法求得的空間場值.

分析圖4(a)可以看出:障礙物后新算法的計算結果與矩量法相比幅度基本一致,障礙物之前的場值振蕩較為劇烈,造成這種現象的原因是地面反射波與直達波疊加相互作用產生的;障礙物之間的場值振蕩是因為障礙物之間存在反射場的多次疊加.障礙物之后場值的抬高是因為隨距離增加,導體障礙物的繞射場從深陰影區逐步過渡到直達波能夠照射到的亮區.由圖4(b)可以看出PEtool算法的計算結果與新算法和矩量法相比起伏的相位出現了偏移,幅度也不完全一致.從而說明邊界處理和相位修正對雙向拋物線方法的計算精度有較大的改善,驗證了新算法在障礙物為理想導體情況下的精確性.

圖4 (網刊彩色)水平極化雙矩形模型場值對比 (a)整體場值;(b)局部場值Fig.4.(color online)The contrasts of fi elds for horizontal polarization under double rectangular models:(a)The whole fi elds;(b)the local fi elds.

實例二:驗證水平極化下圖1中單矩形介質障礙物模型下新方法計算結果的精確性,介質的參數為εr=4,σ=10?5s/m,將發射源的波束寬度更改為30?,單矩形介質障礙物距離波源750 m,障礙物寬度變為20 m,其余參數與實例一相同,計算的接收高度在80 m處場強隨水平距離變化的場強振幅分布如圖5所示,圖5(b)為圖5(a)的局部展開結果.

分析圖5可以看出,當障礙物為均勻介質時,新算法的計算結果與矩量法相比幅度基本一致,尤其是障礙物之后800 m以上,在波形變換凹口處也完全一致,因此可以得出當存在障礙物時,使用分區原理以及介質內的相位分布修正方法可以有效地計算此時空間中的場值大小,也驗證了新算法在障礙物為均勻介質情況下的精確性.

在該實例中,矩量法的計算時間為335 s,新算法的計算時間為5 s,在計算精度幾乎一樣的情況下,新算法的計算時間比矩量法快了將近67倍,從而驗證了新算法計算的有效性.

圖5 (網刊彩色)單矩形介質障礙物場值對比 (a)整體場值;(b)局部場值Fig.5.(color online)The contrasts of fi elds under single rectangular dielectric obstacle:(a)The whole if elds;(b)the local fi elds.

實例三:在實例一參數的基礎上將發射源的波束寬度改為30?,雙矩形寬度改為20 m,其余參數不變,雙矩形均勻介質障礙物的介電常數同實例二,計算水平極化下接收高度為80 m處場強隨水平距離的變化如圖6所示,圖6(b)為圖6(a)的局部展開結果.

從圖6可以看出,在雙矩形均勻介質障礙物模型下,新算法與矩量法的計算結果相比,場值變換趨勢基本一致.由圖6(b)中雙矩形障礙物之間的場值比較可以看出,新算法與矩量法的計算結果相比誤差非常小,從而體現了新算法計算的精確性,此時矩量法所需的時間為1097 s,而新算法的計算時間僅為65 s.當地形變換更為復雜時,矩量法的計算難度會迅速增大甚至無法計算,新算法的優勢性會更進一步體現.在此基礎上,對于雙向情況下新算法對障礙物內的計算誤差有待進一步研究.

圖6 (網刊彩色)雙矩形介質障礙物場值的對比 (a)整體場值;(b)局部場值Fig.6.(color online)The contrasts of fi elds under double rectangular dielectric obstacles:(a)The whole if elds;(b)the local fi elds.

6 結 論

經典的拋物線方法只能用于求解起伏地形下的電波傳播問題的計算,當地面存在障礙物時,尤其是兩種媒質介質參數不同的情況下,雙向拋物線方法無法計算.本文基于雙向拋物線方法的基本原理,提出了一種新的區域分解算法,對存在障礙物時空間中的場值計算進行了區域分解處理,并對均勻介質內雙向拋物線方法的求解給出了計算公式.通過上述方法,使得基于拋物線方程的新算法可以處理地面存在障礙物的情況.在此基礎上,對求解總場時的前向場與后向場的相位以及出介質障礙物邊緣時上下空間的場值疊加時的相位都分別進行了修正,修正方法有效改善了新算法計算的空間中場值的相位匹配性.最后,結合矩量法在近距離下對算法在存在障礙物環境下的電磁傳播計算結果進行了驗證分析,通過實例分析,驗證了基分區原理的雙向拋物線方程新算法計算出場值的精確性以及優越性.

[1]Ozlem O 2009 IEEE Trans.Antenn.Propag.57 2706

[2]Ozlem O,Gokhan A,Mustafa K,Levent S 2011 Comput.Phys.Commun.182 2638

[3]Wang K,Long Y L 2012 IEEE Trans.Antenn.Propag.60 4467

[4]Zhang P,Bai L,Wu Z S,Guo L X 2016 IEEE Trans.Antenn.Propag.Mag.58 31

[5]Wang D D,Xi X L,Pu Y R,Liu J F,Zhou L L 2016 IEEE Trans.Antenn.Wireless Propag.Lett.15 734

[6]Yuan X J,Lin W G 1993 Chin.Phys.Lett.10 57

[7]Omaki N,Yun Z Q,Iskander M F 2012 2012 IEEE International Conference on Wireless Information Technology and Systems(ICWITS)Maui,USA,November 11–16,2012 p1

[8]Kuttler J R 1999 IEEE Trans.Antenn.Propag.47 1131

[9]Donohue D J,Kuttler J R 2000 IEEE Trans.Antenn.Propag.48 260

[10]Beilis A,Tappert F D 1979 J.Acoust.Soc.Am.66 811

[11]Wang Y J,Guo L X,Li Q L 2016 11th International Symposium on Antennas,Propagation and EM Theory(ISAPE)Guilin,China,October 18–21,2016 p404

[12]Ozgun O,Sevgi L 2012 Aces J.27 376

[13]Gokhan A,Levent S 2013 IEEE Trans.Antenn.Propag.Mag.55 244

[14]Zhu J,Yin C Y,Wei Q F 2016 J.Microwaves 32 32(in Chinese)[祝杰,尹成友,魏喬菲 2016微波學報 32 32]

[15]Omak N,Yun Z Q,Iskander M F 2012 Antennas and Propagation Society International Symposium(APSURSI)Chicago,USA,July 8–14,2012,p1

[16]Lu J,Zhou H C 2016 Chin.Phys.B 25 90203

[17]Pvel V,Pvel P 2007 IEEE Antenn.Wireless Propag.Lett.6 152

[18]Sheng X Q 2004 Computational Electromagnetic Theory(Beijing:Science Press)pp49–53(in Chinese)[盛新慶2004計算電磁學要論(北京:科學出版社)第49—53頁]

[19]Zhu J,Yin C Y 2016 J.Microwaves 32 26(in Chinese)[祝杰,尹成友 2016微波學報32 26]

[20]Yin C Y,Zhu J,Wei Q F 2016 37th Progress in Electromagnetics Shanghai,China,August 8–11,2016 p1655

[21]Levy M 2000 Parabolic Equation Methods for Electromagnetic Wave Propagation(London:IEE Press)pp149,287–291

PACS:41.20.JbDOI:10.7498/aps.66.124102

Research and veri fi cation for parabolic equation method of radio wave propagation in obstacle environment?

Wei Qiao-FeiYin Cheng-You?Fan Qi-Meng

(National Key Laboratory of Pulsed Power Laser Technology,Electronic Engineering Institute of PLA,Hefei 230037,China)

14 January 2017;revised manuscript

20 March 2017)

In recent years,the two-way parabolic equation method(2WPE)has been widely utilized for studying the tropospheric ground-wave propagation under the irregular terrain.This algorithm can deal with the in fl uences of the irregular terrain characteristic and the di ff erent electromagnetic parameters of the surface structure on wave propagation.However,there are still some defects in 2WPE method.Firstly,the method considers the irregular terrain and obstacles as a whole,so it cannot deal with the situation where the medium parameters of obstacles and the ground are di ff erent.Secondly,its calculation precision is limited with the inclination of the undulating terrain:if there are obstacles the upper bound of the inclination is easily broken through.Therefore,in this paper,a novel two-way parabolic equation method is proposed for analyzing the radio wave propagation in obstacle environment.According to the principle of domain decomposition,the obstacle zones are divided into two domains in the new algorithm,and the two subdomains are calculated,respectively.Meanwhile,in order to avoid the calculation error caused by the abrupt truncation of the obstacle zone,the fi eld at the upper boundary of obstacles is modi fi ed to ensure the continuity of tangential fi eld.To further improve the accuracy of the new algorithm,according to the historical transmission paths,we exactly retrieve the phases of each forward and backward wave,especially when stepping in and out of the obstacles.Furthermore,the method of moment(MoM)is used to verify the calculation accuracy of the new algorithm in obstacle environment.Although the accuracy of the MoM is very high,it also requires a great deal of calculation resources:it can only be employed to compute the fi elds in short distance.To overcome the difficulty,we use the image principle in the obstacle environment and do not subdivide the ground into segments;therefore the veri fi cation accuracy can be improved.On this basis,to unify the source setting of the new algorithm and the MoM,the equivalent source model is used to set the initial fi eld.Finally,through numerical experiments,the simulation results of both methods agree very well,so the e ff ectiveness of the boundary correction and the phase correction which are presented in this paper are both veri fi ed.The accuracy and superiority of the new algorithm in obstacle environment are also demonstrated.To sum up,the novel two-way parabolic equation method can be used to accurately calculate the fi eld of the space in the obstacle environment,and lays the foundation for the fi eld calculation of radio wave propagation in real environment.

parabolic equation,radio wave propagation,domain decomposition,method of moment

10.7498/aps.66.124102

?總裝備部預研基金(批準號:51333020201)資助的課題.

?通信作者.E-mail:cyouyin@sina.com

?2017中國物理學會Chinese Physical Society

http://wulixb.iphy.ac.cn

*Project supported by the General Equipment Department Pre-Research Foundation,China(Grant No.51333020201).

?Corresponding author.E-mail:cyouyin@sina.com