障礙物環(huán)境下電波傳播的分區(qū)計(jì)算方法

李安琪 尹成友 甘泳機(jī)

（國(guó)防科技大學(xué)電子對(duì)抗學(xué)院，合肥 230037）

引 言

當(dāng)今社會(huì)，在各個(gè)領(lǐng)域，不論是信息通信還是國(guó)防建設(shè)，研究不同地形環(huán)境下的電波傳播問(wèn)題變得越來(lái)越重要，因此考慮到實(shí)際地形情況，研究存在障礙物下的電波傳播分析與計(jì)算非常必要.近年來(lái)，用于計(jì)算預(yù)測(cè)電波傳播的方法主要有：模式理論法[1-3]、光學(xué)近似方法[4-5]、半經(jīng)驗(yàn)公式法（ITU-R標(biāo)準(zhǔn)）[6]、時(shí)域有限差分法[7]、拋物線方程(parabolic equation, PE)近似法[8-10]. 其中PE是由波動(dòng)方程近似化簡(jiǎn)之后得到的，屬于一種迭代步進(jìn)方法，是研究電波傳播問(wèn)題時(shí)最常用的方法，具有穩(wěn)定性好和計(jì)算精度高的優(yōu)點(diǎn). 其他常用的方法，如模式理論法只適用于某些特定的情況，射線追蹤法的精度不高，矩量法(method of moment, MoM)的計(jì)算量大，只適用于計(jì)算短距離傳播的情況.

PE求解電波傳播問(wèn)題如圖1所示. 在使用經(jīng)典PE方法求解時(shí)將地面與障礙物視為整體，等效為起伏地形來(lái)處理，這種計(jì)算方法存在一定的缺陷：一是沒(méi)有考慮障礙物與普通地形不同，兩者介電常數(shù)存在差異；二是障礙物的存在導(dǎo)致電波傳播的仰角過(guò)大，普通的PE近似不再成立[11]；三是沒(méi)有考慮障礙物內(nèi)部場(chǎng)值變化對(duì)空間中場(chǎng)值分布的影響. 因此實(shí)際上計(jì)算整個(gè)空間內(nèi)場(chǎng)強(qiáng)的分布情況時(shí)，我們需要額外考慮存在障礙物的情況.

圖1 PE法求解復(fù)雜環(huán)境電波傳播示意圖Fig. 1 Schematic diagram of wave propagation in complex environment

目前對(duì)于這個(gè)問(wèn)題的主要處理方法來(lái)自于文獻(xiàn)[12]，將障礙物上方空間計(jì)算的下邊界往障礙物之下延伸一小段，內(nèi)部空間的計(jì)算也往障礙物之上延伸至無(wú)窮遠(yuǎn)處，這樣障礙物的上方和內(nèi)部區(qū)域都變成了無(wú)窮大的區(qū)域，均可以使用離散混合傅里葉變換(discrete mixed Fourier transform, DMFT)進(jìn)行求解，同時(shí)將連接處區(qū)域的場(chǎng)值使用兩種算法的平均值代替. 顯然這種方法存在一定的欠缺：向下延伸一小段的長(zhǎng)度并沒(méi)有給出計(jì)算方法，也沒(méi)有較好的對(duì)這一小段長(zhǎng)度進(jìn)行求解的辦法；該方法選取了邊界上的一小段求平均值來(lái)保持場(chǎng)值連續(xù)，但這段的距離也沒(méi)有較好的辦法直接得出，因此計(jì)算精度不高；沒(méi)有考慮到兩個(gè)區(qū)域之間的互相影響，所以這種方法存在一定缺陷，我們需要尋求其他的方法來(lái)解決這個(gè)問(wèn)題.

而當(dāng)前已知的計(jì)算電波傳播PE的求解方法主要有：有限差分(finite difference, FD)方法[13-14]、有限元方法[15-17]、分步傅里葉變換(split step Fourier transform, SSFT)方法[18-21]. 其中，F(xiàn)D法和有限元方法都屬于數(shù)值計(jì)算方法. SSFT方法計(jì)算速度快，計(jì)算性能穩(wěn)定，遠(yuǎn)距離電波傳播問(wèn)題上優(yōu)勢(shì)明顯.有限元方法適用于低頻信號(hào)，不適用于大氣波導(dǎo)問(wèn)題的求解，且計(jì)算量隨著距離的增大而增大，也不適用于遠(yuǎn)距離問(wèn)題的求解. 雖然FD方法和SSFT方法相比，計(jì)算時(shí)間和計(jì)算效率都不高，但可以更為靈活地處理邊界條件.

本文將區(qū)域分解思想應(yīng)用于該問(wèn)題的求解，各子域分別采用不同算法進(jìn)行計(jì)算，在目標(biāo)區(qū)域采用差分法求解，并對(duì)邊界進(jìn)行處理，保證邊界場(chǎng)值連續(xù)，在原先PE方法上提高了計(jì)算效率，擴(kuò)展了適用范圍.

1 問(wèn)題描述

區(qū)域分解算法是指將整個(gè)區(qū)域上的大問(wèn)題分解為若干子區(qū)域上的小問(wèn)題分別進(jìn)行求解，獲得每個(gè)小問(wèn)題的結(jié)果之后，通過(guò)級(jí)聯(lián)即可獲得整個(gè)區(qū)域上的解[22]. 這樣每個(gè)區(qū)域上的問(wèn)題不但可以獨(dú)立求解，而且可以選擇適合于該子域的最有效的方法，不同子域選擇不同解法，并且求解規(guī)模較小. 假設(shè)地面為理想導(dǎo)體的情況下，對(duì)于電波傳播來(lái)說(shuō)，在遇到障礙物時(shí)，可以將分析的區(qū)域分解為無(wú)障礙物區(qū)和有障礙物區(qū)，具體如圖2所示的 Ω1、Ω2區(qū).

圖2 區(qū)域分解示意圖Fig. 2 Schematic diagram of domain decomposition

在 Ω1區(qū)電波傳播滿足PE，其上邊界采用吸收邊界，可以認(rèn)為其到無(wú)窮遠(yuǎn)，下邊界連接在障礙物頂部，可以采用阻抗邊界條件，因此該區(qū)域PE的求解可以采用DMFT步進(jìn)求解.

在Ω2區(qū)電波傳播一樣滿足PE，其下邊界是理想導(dǎo)體，對(duì)于水平極化波是第一類邊界條件，而垂直極化屬于第二類邊界條件；其上邊界與Ω1區(qū)的下邊界相連，同樣滿足阻抗邊界條件. 由于Ω2區(qū)上下封頂，PE不便于使用DMFT進(jìn)行求解，因此該區(qū)域可以采用FD求解.

2 雙向PE DMFT的計(jì)算

2.1 雙向PE的推導(dǎo)

假設(shè)求解時(shí)諧場(chǎng)為 exp(jωt)，其中ω為角頻率.在二維直角坐標(biāo)系中我們對(duì)場(chǎng)量 Ψ 在不同極化情況下分別定義為方向. 則場(chǎng)量Ψ 滿足二維波動(dòng)方程

電磁波沿軸向傳播，這里選擇x軸為正向傳播

式中：k是媒質(zhì)均勻環(huán)境中的波數(shù)，Ω1區(qū)中k為自由空間波數(shù)k0，Ω2區(qū) 中k為障礙物介質(zhì)波數(shù)代表修正后的大氣折射率. 在上述理論分析的基礎(chǔ)上，假設(shè)折射率n不 隨距離x變化而變化，將式(2)因式分解后得

分別引入兩種不同的相位解調(diào)因子[3]

式中，ΨfΨb分別代表前向場(chǎng)和后向場(chǎng). 因此本文采用的都是雙向PE算法，將式(4)分別代入因式分解后所得式，求得前向傳播和后向傳播PE如下：

要得到方程(2)的解，方程(5)的前向和后向就必須在同一系統(tǒng)中進(jìn)行求解，給定初始值和底部、頂部的邊界值，就可以用步進(jìn)法求解. 在求解完成后我們還會(huì)采用對(duì)前后場(chǎng)進(jìn)行相位修正的辦法來(lái)進(jìn)一步改善計(jì)算精度，解決整體場(chǎng)值相位不匹配的問(wèn)題[23].

2.2 邊界條件的討論

前面推導(dǎo)中都假設(shè)地面是理想導(dǎo)體情況下處理的，所以在平坦地面區(qū)域可以直接采用SSFT方法求解. 但在介質(zhì)障礙物區(qū)域，對(duì)于 Ω1區(qū)的上邊界，該假設(shè)不再成立，因此需要另外處理. 一種簡(jiǎn)便的方法是采用阻抗邊界條件來(lái)獲取邊界上的場(chǎng)值關(guān)系. 媒質(zhì)1和媒質(zhì)2的分布如圖3所示.

圖3 阻抗邊界示意圖Fig. 3 Schematic diagram of impedance boundary

2.2.1 廣義阻抗邊界條件

對(duì)于水平極化波情況：Ψ=Ey，由Maxwell旋度方程可得

因此有邊界條件為

假設(shè)媒質(zhì)2是均勻的，其復(fù)相對(duì)介電常數(shù)為

同樣對(duì)媒質(zhì)2采用PE求解，假設(shè)在分界面下初始場(chǎng)為0，近似認(rèn)為在分段內(nèi)不是x的函數(shù)，對(duì)式(4)中前向傳播方程中的x分量進(jìn)行Laplace變換，則

式中：k是自由空間波數(shù)；U(s,z) 為u(x,z)關(guān)于x的Laplace變換結(jié)果.

因?yàn)槊劫|(zhì)2是均勻、有耗的，因此其電磁波在z→?∞是有界的，甚至是指數(shù)衰減的，有

α(s)應(yīng)該具有負(fù)虛部，將式(10)代入(9)，并利用譜域的思想，得到

這樣媒質(zhì)2在邊界處滿足

利用媒質(zhì)1與媒質(zhì)2的邊界切向場(chǎng)的連續(xù)性，可以得到媒質(zhì)1邊界滿足

根據(jù)上面的分析，對(duì)于水平極化情況，比較式(12)和(13)可得到

同理可得垂直極化情況. 上述推導(dǎo)的關(guān)系式稱為廣義阻抗邊界條件.

2.2.2 Leontovich阻抗邊界條件

上述阻抗邊界條件是在譜域下得到的，假設(shè)一個(gè)平面波入射到表面上：

對(duì)其進(jìn)行Laplace變換得到

當(dāng)對(duì)應(yīng)的電磁波為消逝波，我們可以忽略. 當(dāng)掠入射角較小時(shí)，對(duì)應(yīng)窄角PE情況，對(duì)于水平極化波，有

同理可得垂直極化波情況. 顯然這種邊界條件適合 Ω1區(qū)的下邊界.

2.2.3 Ω2區(qū)域上邊界條件的確定

對(duì)于Ω2區(qū)的上邊界來(lái)說(shuō)，其邊界由圖3可以看出，此時(shí)它與 Ω1區(qū) 相比，兩種媒質(zhì)正好顛倒相當(dāng)于空氣，不具有大的模值時(shí)，上述Leontovich阻抗邊界條件不再適用. 但由前面分析知，在介質(zhì)的分界面上，對(duì)于水平極化波有如式（7）所示的情況. 同理可得垂直極化波情況.

同時(shí)由上述分析可得，Ω1區(qū)的下邊界滿足Leontovich 阻抗邊界條件. 由于式（17）中δ(s)與s無(wú)關(guān)，對(duì)式（13）進(jìn)行Laplace反變換可以得到水平極化波情況下 Ω1區(qū)的下邊界條件，同理可得垂直極化波情況下的Ω1區(qū)下邊界條件，因此可以表示為

假設(shè)Ω1區(qū) 貼近邊界處折射率由邊界切向場(chǎng)的連續(xù)性和式(18)可得 Ω2區(qū)的上邊界條件為

由于上下區(qū)域在邊界處交匯，邊界處的值在步進(jìn)過(guò)程中會(huì)同時(shí)代入兩種算法進(jìn)行迭代計(jì)算，我們選擇將兩種算法計(jì)算出的邊界值取平均之后，再作為初值代入下一次的迭代過(guò)程. 這樣既確保了邊界切向場(chǎng)連續(xù)，又考慮了上下區(qū)域電波傳播過(guò)程中場(chǎng)值之間發(fā)生的耦合作用，避免突然截?cái)嗨鶎?dǎo)致的計(jì)算誤差.

2.3 雙向PE的DMFT方法求解

式(4)中對(duì)Q因子的處理方法有很多種，這里主要從寬角近似出發(fā)，將Q按Feit-Fleck模型展開(kāi)，對(duì)應(yīng)的Q展開(kāi)式為[24-25]

將Q代入到式(4)中的Feit-Fleck寬角近似下的雙向波動(dòng)方程，通過(guò)DMFT方法求得步進(jìn)公式uf(x+?x,z),ub(x??x,z)：

3 PE的FD計(jì)算方法

3.1 寬角PE的差分近似

使用FD法來(lái)求解PE，實(shí)際就是將PE中的一階偏微分和二階偏微分都用差分形式改寫(xiě)，從而通過(guò)步進(jìn)方法求解，獲得整個(gè)計(jì)算區(qū)域內(nèi)的場(chǎng)值. 記：

式 中，zp為 垂 直 網(wǎng) 格 點(diǎn). 令x0、x1、···、xm、···、xN表示水平網(wǎng)格點(diǎn)，x方向即為步進(jìn)方向，網(wǎng)格劃分如圖4所示. 為了進(jìn)行中心差商處理，設(shè)的中點(diǎn)為則可以得到點(diǎn)處的對(duì)x一階差分和對(duì)z的二階差分:

在PE的求解中，首先需要得到關(guān)于u的標(biāo)量方程. 標(biāo)準(zhǔn)PE的限制主要在于難以準(zhǔn)確計(jì)算大傳播仰角的波束，由此可見(jiàn)標(biāo)準(zhǔn)PE是一個(gè)窄角處理方法. 因此，為實(shí)現(xiàn)寬角情況下的電波傳播計(jì)算的PE方法，采用pade近似得到Q的近似表達(dá)式，即

取a=0.75,b=0.25，稱之為Claerbout逼近. 下面開(kāi)始對(duì)場(chǎng)強(qiáng)進(jìn)行求解.

圖4 FD網(wǎng)格劃分示意圖Fig. 4 Schematic diagram of finite difference mesh generation

3.2 電磁波前向傳播的情況

首先討論前向差分時(shí)前向傳播方程，利用Claerbout展開(kāi)：

將式（26）代入（4）中的正向傳播方程后進(jìn)一步化簡(jiǎn)得到計(jì)算方程

整理化簡(jiǎn)可得

式(28)為第p=1,···,N?1個(gè)方程. 為了完整起見(jiàn)，還需要包括上下邊界的差分方程.

3.2.1 理想導(dǎo)體地面下邊界的處理

對(duì)于水平極化，由于其值恒為0，相對(duì)簡(jiǎn)單，其邊界方程為

對(duì)于垂直極化，則需要導(dǎo)出寬角PE的邊界條件.

3.2.2 介質(zhì)內(nèi)部上邊界的處理

此時(shí)仍采用寬角差分方程，由式(27)和(4)，取z=N,x=(xm?1+xm)/2點(diǎn)進(jìn)行差分處理，依照式（19）有

將式（23）和（32）代入式（27）得到計(jì)算方程

這樣式(28)、(31) 和(33)共同構(gòu)成了水平極化時(shí)0～N的前向差分PE.

PE差分解遞推結(jié)構(gòu)用矩陣可表達(dá)為

式中：Um是在xm處的場(chǎng)向量，

對(duì)于水平極化波，Am是一個(gè)三對(duì)角矩陣，

這種Crank-Nicolson FD結(jié)構(gòu)，由于具有三對(duì)角特征，方程計(jì)算方便，而且具有非常好的數(shù)值特性，可以采用追趕算法快速計(jì)算，則方程(34)可以改寫(xiě)為

根據(jù)式(36)和(37)求解該矩陣，得到前向傳播的Um.

3.3 電磁波后向傳播的情況

后向傳播的迭代情況同前向傳播一致，所以后向傳播的迭代矩陣和前向的完全一致.

通過(guò)前面的推導(dǎo)我們得到了電波傳播矩陣，要求解差分方程，傳統(tǒng)的方法是對(duì)Am求逆，但計(jì)算量大且速度慢，目前比較方便的方法是 LU算法，就是將正定的三對(duì)角系數(shù)矩陣A分解為上三角矩陣U和下三角矩陣L. 本文在 LU算法的基礎(chǔ)上采用追趕算法，這種算法非常簡(jiǎn)單，且具有數(shù)值穩(wěn)定的特點(diǎn)，在求解三對(duì)角矩陣方程中應(yīng)用也很廣泛[26].

將其具體應(yīng)用到FD矩陣方程中，實(shí)現(xiàn)步驟如下：

1)計(jì)算Vm=BmUm?1；

2)對(duì)Am進(jìn)行Crout分解，參數(shù)均已知；

3)求解方程組Ly=Vm；

4)求解方程組Ux=y，最終得到x.

通過(guò)以上4個(gè)步驟可由第m?1步的場(chǎng)值導(dǎo)出第m步的場(chǎng)值，因此只要給定每一步的初始場(chǎng)即可得到任意位置的場(chǎng)值.

4 仿真結(jié)果

4.1 單矩形介質(zhì)障礙物情況下的電波傳播分析

圖5為單矩形障礙物仿真示意圖，假設(shè)地面為理想導(dǎo)體，發(fā)射源工作頻率30 MHz、高80 m、波束寬度 20?、水平極化，輻射源等效輻射功率（equivalent radiation power, EIRP）為W，Z0為自由空間波阻抗. 矩形模型位于距離發(fā)射源752 m處，介質(zhì)障礙物模型寬16 m、高100 m，媒質(zhì)的電導(dǎo)率為10?5、相對(duì)介電常數(shù)為4. 計(jì)算得到的整體場(chǎng)強(qiáng)分布如圖6所示. 從圖6可以看到，地面反射與直達(dá)波之間形成干涉，導(dǎo)致波瓣分裂現(xiàn)象. 此外，介質(zhì)障礙物產(chǎn)生的反射波與直達(dá)波形成干涉，導(dǎo)致各波瓣上形成干涉條紋. 障礙物內(nèi)部也有相應(yīng)的干涉存在.

圖5 單矩形障礙物仿真示意圖Fig. 5 Simulation diagram of a single rectangular medium obstacle

圖6 單矩形介質(zhì)障礙物情況下的Ey分布圖Fig. 6 Distribution diagram of Ey in the case of a single rectangular medium obstacle

圖7為三種算法在相同給定觀察高度80 m下的整體場(chǎng)值比較，其中 PE算法是文獻(xiàn)[12]計(jì)算結(jié)果，新PE算法是本文所提新算法結(jié)果. 可以看出，PE算法和新PE算法結(jié)果都與MoM的吻合較好，但新PE算法似乎吻合更好.

圖7 單矩形介質(zhì)障礙物情況下80 m處的Ey比較圖Fig. 7 Comparison of Ey at 80 m in the case of a single rectangular medium obstacle

為進(jìn)一步觀察這兩種方法在不同區(qū)域的場(chǎng)值計(jì)算精度，對(duì)圖7提取局部場(chǎng)值如圖8所示. 由圖8（a）障礙物之前的局部場(chǎng)值可以看出，在等效源模型的基礎(chǔ)上，障礙物之前新PE算法相位修正后的計(jì)算結(jié)果與MoM計(jì)算結(jié)果幾乎一致，驗(yàn)證了新PE算法計(jì)算結(jié)果的精確性. 圖8（b）中，新PE算法與MoM在障礙物內(nèi)部區(qū)域計(jì)算的場(chǎng)值略有差異，而和文獻(xiàn)[12]中的PE算法近乎一致. 這種誤差來(lái)源于PE算法本身的缺陷，因?yàn)檩椛湓淳哂幸欢ǖ牟ㄊ鴮挾龋ㄔ凑丈涞秸系K物表面時(shí)不是垂直入射，而且角度偏差越大，障礙物內(nèi)的計(jì)算誤差越大. 圖8（c）中，靠近障礙物部分時(shí)新PE算法與MoM的場(chǎng)值略有差異，但越是遠(yuǎn)離障礙物時(shí)，二者的結(jié)果就越一致. 在計(jì)算過(guò)程中，MoM平均用時(shí)159.6 s，文獻(xiàn)[12]中的PE算法用時(shí)16.9 s，而新PE算法則用時(shí)7.6 s，計(jì)算速度優(yōu)勢(shì)明顯.

圖8 單矩形介質(zhì)障礙物情況下Ey的局部場(chǎng)值比較Fig. 8 Comparison of local field values of Ey in the case of a single rectangular medium obstacle

4.2 雙矩形介質(zhì)障礙物情況下的電波傳播分析

圖9為雙矩形障礙物仿真示意圖，假設(shè)地面為理想導(dǎo)體，發(fā)射源工作頻率30 MHz、高80 m、波束寬度 20?、 水平極化，輻射源EIRP為W，Z0為自由空間波阻抗. 雙矩形模型位于距離發(fā)射源750 m和900 m處，介質(zhì)障礙物模型寬22 m、高100 m，媒質(zhì)的電導(dǎo)率為1 0?5、相對(duì)介電常數(shù)為3. 雙矩形介質(zhì)障礙物下的整體場(chǎng)值分布如圖10所示. 可以看出，直達(dá)波、地面反射波以及介質(zhì)障礙物的反射波形成的波束與圖6類似.

圖9 雙矩形障礙物仿真示意圖Fig. 9 Simulation diagram of double rectangular medium obstacles

圖10 雙矩形介質(zhì)障礙物下的Ey分布圖Fig. 10 Distribution diagram of Ey in the case of double rectangular medium obstacles

圖11所示為接收高度為80 m處雙矩形介質(zhì)障礙物下新PE算法和MoM兩種算法得到的整體場(chǎng)值比較. 可以看出，在雙矩形介質(zhì)模型下，兩種算法所得整體場(chǎng)值依然吻合得較好.

圖11 雙矩形介質(zhì)障礙物情況下80 m處的Ey比較圖Fig. 11 Comparison of Ey at 80 m in the case of double rectangular obstacles

圖12所示為雙矩形介質(zhì)障礙物下新PE算法和MoM兩種算法得到的局部場(chǎng)值比較. 可以看出，和單矩形模型類似，在障礙物之前新PE算法相位修正后的計(jì)算結(jié)果與MoM計(jì)算結(jié)果幾乎一致，障礙物內(nèi)部區(qū)域計(jì)算的場(chǎng)值略有差異，障礙物后的靠近障礙物部分兩種算法的場(chǎng)值略有差異，隨著離障礙物距離越遠(yuǎn)，二者的結(jié)果越一致.

圖12 雙矩形介質(zhì)障礙物情況下Ey的局部場(chǎng)值比較Fig. 12 Comparison of local field values of Ey in the case of double rectangular medium obstacles

5 結(jié) 論

本文根據(jù)區(qū)域分解原理，提出了一種針對(duì)障礙物環(huán)境下新的雙向PE算法. 該算法在無(wú)障礙物區(qū)域采用普通的DMFT算法，在障礙物區(qū)域進(jìn)行區(qū)域分解，在障礙物上方的無(wú)限大區(qū)域采用DMFT算法，在障礙物內(nèi)部的封閉區(qū)域則采用FD算法，有效解決了障礙物中電波傳播問(wèn)題. 在相同波源、天線等參數(shù)設(shè)置的基礎(chǔ)上，用較為精確的MoM驗(yàn)證了近距離下新雙向PE算法在存在單矩形和雙矩形障礙物情況下的電波傳播計(jì)算的精確性和有效性. 該算法與文獻(xiàn)[12]已有的雙向PE算法相比，新算法理論分析更為清晰，計(jì)算速度和精度都有一定提高. 本文研究的是規(guī)則形狀障礙物地形下的電波傳播問(wèn)題，實(shí)際障礙物不一定是規(guī)則的形狀，計(jì)算更加復(fù)雜，因此對(duì)不規(guī)則障礙物地形下的電波傳播問(wèn)題還需要作進(jìn)一步研究.