PE模型SSFT解電波傳播問題的虛部增量法

肖金光，周新力，劉曉娣

（1.海軍航空工程學(xué)院 研究生管理大隊(duì)，山東 煙臺(tái)264001；2.海軍航空工程學(xué)院 電子信息工程系，山東 煙臺(tái)264001）

0 引 言

近年來，電磁環(huán)境問題越來越為軍民用戶所重視。采用數(shù)值法求解電波傳播問題時(shí)，由Maxwell方程簡(jiǎn)化而來的橢圓型波動(dòng)方程，一點(diǎn)的求解意味著必須同時(shí)求出計(jì)算域內(nèi)全部點(diǎn)的解，實(shí)時(shí)求解較大計(jì)算域時(shí)計(jì)算量將極其龐大，拋物型方程 （parabolic equation，PE）忽略電波后向傳播而只計(jì)算沿某坐標(biāo)軸以較小角度形成的錐形區(qū)域內(nèi)的前向傳播。其解法有有限差分 （finite difference，F(xiàn)D）法、有限元 （finite element，F(xiàn)E）法和分步傅里葉 （split step Fourier transform，SSFT）法等。其中SSFT 能夠方便的實(shí)現(xiàn)步進(jìn)求解，步長(zhǎng)選取靈活，適用于電波遠(yuǎn)距離傳播計(jì)算［1］。基于PE 的SSFT 解法，美國(guó)研制了高級(jí)傳播模型（advanced propagation model，APM），并以APM 為計(jì)算核心研制了 “高級(jí)折射效應(yīng)預(yù)測(cè)系統(tǒng)” （advanced refractive effects prediction system，AREPS），用于美軍的戰(zhàn)場(chǎng)電磁態(tài)勢(shì)評(píng)估［2，3］。國(guó)內(nèi)眾多研究所和院校等研究機(jī)構(gòu)的專家學(xué)者對(duì)大氣波導(dǎo)中PE計(jì)算模型關(guān)鍵技術(shù)、應(yīng)用及實(shí)驗(yàn)驗(yàn)證進(jìn)行了持續(xù)的跟蹤研究［4－6］。

對(duì)于對(duì)流層電波傳播的拋物方程求解，需要知道初始場(chǎng)、上下邊界條件和大氣折射情況。其中上邊界需要滿足第三類邊界條件——Sommerfeld輻射條件，PE 的SSFT 解法需采取措施避免邊界上的強(qiáng)反射影響到計(jì)算域。傳統(tǒng)辦法是在需要求解區(qū)域上加吸收層，層內(nèi)采用窗函數(shù)進(jìn)行濾波［2，6，7］，這種方法的優(yōu)點(diǎn)是概念清楚，步驟清晰易理解，但在每一SSFT 步進(jìn)上都要進(jìn)行一次濾波。不難看出，窗函數(shù)濾波的功能和計(jì)算方法都具有單一性，因此本文提出一種等效方法——虛部增量法，將窗函數(shù)濾波的過程綜合到環(huán)境折射效應(yīng)計(jì)算中去，不再需要單獨(dú)進(jìn)行濾波計(jì)算。在理論推導(dǎo)的基礎(chǔ)上，設(shè)計(jì)了實(shí)現(xiàn)步驟，分析了該方法的計(jì)算性能優(yōu)勢(shì)，并對(duì)2種方法的等效性做了仿真分析。

1 PE模型的SSFT解及吸收層窗函數(shù)法濾波

設(shè)電磁場(chǎng)時(shí)諧因子為e－iωt，忽略后向傳播并作近軸條件下的窄角近似可得標(biāo)準(zhǔn)拋物方程 （standard parabolic equation，SPE）為

式中：u（x，z）——水平距離x、垂直高度z 處的場(chǎng)解，k——自由空間波數(shù)，n——大氣折射率。

其適于求解遠(yuǎn)距電波傳播問題的SSFT 解為［3］

式中：Δx——一次步進(jìn)的水平距離變化，顯然式 （2）是一種步進(jìn)迭代算法。F、F－1分別表示傅里葉變換和傅里葉反變換，變換的點(diǎn)數(shù) N 由 Nyquist 準(zhǔn)則確定［4］，如式（3）所示

式中：zmax和pmax——輸出高度和變換域最大值。

對(duì)于考慮地球曲率的修正折射指數(shù)M（x，z）有

將式 （4）代入式 （2）可得

式（5）中前后兩指數(shù)項(xiàng)分別表征了傳播媒質(zhì)的折射效應(yīng)和障礙物繞射效應(yīng)，不妨分別稱之為折射因子和繞射因子。這樣只要知道初始場(chǎng)、上下邊界條件就可以按照式（5）步進(jìn)求解。初始場(chǎng)可以通過天線方向圖的傅里葉逆變換求得。如果傳播環(huán)境下邊界是光滑平面，還可將FFT 進(jìn)一步簡(jiǎn)化為正弦變換，若下邊界為粗糙水面或陸地，可將下表面視為阻抗邊界，采用DMFT 技術(shù)，但這2 種情況下，PE模型SSFT 解法的基本求解思路不變，在此不展開論述。

應(yīng)用PE進(jìn)行電波傳播計(jì)算時(shí)，上邊界必須滿足Sommerfeld輻射條件，即場(chǎng)在無窮遠(yuǎn)處變?yōu)榱?，然而?jì)算域不可能是無限大的，必須在需要求解區(qū)域的上方添加一個(gè)吸收層，層內(nèi)電波完全被吸收，避免截?cái)噙吔缟系陌l(fā)生強(qiáng)反射造成需求解區(qū)域場(chǎng)值的 “污染”，造成計(jì)算誤差［8］。常用吸收邊界條件有基于媒質(zhì)的和基于差分方程的兩類，前者通過添加一個(gè)損耗媒質(zhì)將向外傳播的電波吸收［9］。PE 的SSFT 算法中只能使用損耗媒質(zhì)法，在需求解區(qū)域上加吸收層，用窗函數(shù)對(duì)吸收層內(nèi)的場(chǎng)進(jìn)行濾波，保持場(chǎng)在需求解區(qū)域內(nèi)保持不變，而在吸收層平滑衰減至0，即在最大高度處完全吸收。APM 中使用窗函數(shù)為 Cosine－tapered（Tukey）窗［2］

窗函數(shù)法相當(dāng)于給折射指數(shù)添加了復(fù)數(shù)部分，使得吸收層中的傳播媒介成為損耗媒介［8］。如果能夠?qū)⒋昂瘮?shù)的濾波作用等效成折射指數(shù)的一個(gè)增量，從而，將濾波的步驟綜合到折射效應(yīng)計(jì)算中去，可以減少計(jì)算步驟和計(jì)算量。下一節(jié)對(duì)這一思路進(jìn)行理論推導(dǎo)。

立冬始降溫，天冷好養(yǎng)腎。按照中醫(yī)季節(jié)養(yǎng)生的理論，冬季對(duì)應(yīng)五臟中的腎。立冬時(shí)心肺氣弱，腎氣強(qiáng)盛，飲食宜減辛苦，以養(yǎng)腎氣。在飲食上依然要遵循“秋冬養(yǎng)陰”的原則，也就是說，少食生冷之物，但也不宜進(jìn)食燥熱之物，有的放矢地食用一些滋陰潛陽、熱量較高的膳食為宜，同時(shí)也要多吃新鮮蔬菜以避免維生素的缺乏。

2 虛部增量法的理論推導(dǎo)

不妨設(shè)初始場(chǎng)和第n步的解分別為u0（x，z）和un（x，z）。濾波過程可以并入到式 （6）中

相當(dāng)于在每次步進(jìn)之前，先進(jìn)行濾波。因?yàn)榇昂瘮?shù)的濾波區(qū)域?qū)嵸|(zhì)上只作用于吸收層，離散計(jì)算時(shí)，只對(duì)吸收層計(jì)算，而該部分并不作為計(jì)算結(jié)果輸出，因此濾波運(yùn)算完全可以前移到上一次步進(jìn)結(jié)束而不會(huì)對(duì)輸出計(jì)算結(jié)果造成任何影響，即

如果折射率剖面在水平上是不變的，即折射率n不隨水平距離x變化，那么只需要將窗函數(shù)w（z）和折射因子整合為一個(gè)因子即可，每當(dāng)折射率剖面變化 （多折射率剖面）時(shí)，需要重新計(jì)算該因子。對(duì)于多折射率剖面的情況，APM 使用線性插值法求解各個(gè)步驟上的折射率剖面，則在每一步驟上該因子都需要重新計(jì)算，即濾波過程每次都需要執(zhí)行N／4 次乘法，對(duì)應(yīng)于M 步PE，就要執(zhí)行NM／4乘法。

不妨令

則

這樣，由w（z）可求得ΔM ，必須要指出的是，ΔM 是虛數(shù)，而M（x，z）是實(shí)數(shù)，但由于式 （8）環(huán)境因子指數(shù)部分最終結(jié)果是復(fù)數(shù)，因此在計(jì)算上從實(shí)數(shù)到復(fù)數(shù)的轉(zhuǎn)變沒有影響，令

這一變化相當(dāng)于給修正折射指數(shù)一個(gè)小虛部增量，相當(dāng)于添加了損耗媒質(zhì)，使得電磁波在該吸收層中完全被吸收，故將這種處理方法稱為修正折射率虛部增量法，簡(jiǎn)稱虛部增量法。

3 虛部增量法的算法流程設(shè)計(jì)及計(jì)算量分析

總結(jié)虛部增量法的實(shí)現(xiàn)步驟為：

步驟1 拋物方程模型及參數(shù)初始化，設(shè)定初始場(chǎng)、傅里葉變換的點(diǎn)數(shù)和下邊界條件；

步驟2 根據(jù)窗函數(shù)求得虛部增量；

步驟3 根據(jù)當(dāng)前步進(jìn)上的大氣修正折射指數(shù)及步驟2所得的虛部增量，求得等效大氣修正折射指數(shù)；

步驟4 拋物方程步進(jìn)求解和最高點(diǎn)場(chǎng)值置零；

步驟5 重復(fù)步驟3、步驟4，直至達(dá)到計(jì)算域距離終點(diǎn)。

基于虛部增量法的PE計(jì)算流程如圖1所示。

通過將窗函數(shù)等價(jià)為折射率指數(shù)的虛部增量，一次濾波過程只需要N／4次加法即可實(shí)現(xiàn)，對(duì)應(yīng)于M 步PE，需執(zhí)行NM／4次加法運(yùn)算。將NM／4復(fù)數(shù)乘法用加法運(yùn)算即可實(shí)現(xiàn)，這一優(yōu)化在進(jìn)行電波傳播實(shí)時(shí)求解，追求運(yùn)算速度時(shí)計(jì)算量的節(jié)約達(dá)到的效果非常明顯。對(duì)于插值獲得大量折射率剖面，如每一步進(jìn)上的折射率剖面，則有M 個(gè)修正折射指數(shù)剖面，使用窗函數(shù)濾波法和虛部增量法的運(yùn)算量 分 別 為NM2／4 次復(fù)數(shù)乘法和NM2／4 復(fù)數(shù)加法，2 種方法所需運(yùn)算量的差別是可觀的，特別是長(zhǎng)距離電波傳播時(shí)PE步數(shù)較多情況，對(duì)降低PE 模型進(jìn)行電波傳播求解的計(jì)算量很有意義。

圖1 基于虛部增量法的PE計(jì)算流程

4 仿真分析

假設(shè)輻射源工作頻率為10GHz，天線類型為高斯型，歸一化輻射強(qiáng)度，天線仰角為0，天線高度10 m，輸出高度為100m，則由PE模型的參數(shù)設(shè)置原則，在需要求解區(qū)域上方添加高度為100／3 m 的吸收層，根據(jù)式 （3），垂直剖面上需要1024點(diǎn)FFT 運(yùn)算。濾波窗函數(shù)為式 （4）所示Tukey窗，不妨設(shè)海面蒸發(fā)波導(dǎo)高度為14m，風(fēng)速為5m／s，采用米勒布朗模型求解粗糙海面上的邊界阻抗，用DMFT 法求解PE模型［10］。由所述虛部增量法，可得修正折射率的虛部增量。修正折射率剖面M（x，z）是實(shí)數(shù)，添加虛部增量后為M′（x，z）成為復(fù)數(shù)。實(shí)質(zhì)上，蒸發(fā)波導(dǎo)折射率剖面M（x，z）就是M′（x，z）的實(shí)部，因而，只繪出了M′（x，z），如圖2所示。

圖2 虛部增量修正折射率

采用窗函數(shù)法和虛部增量法兩種方法分別對(duì)折射因子進(jìn)行處理，并與濾波前的復(fù)折射因子做比較如圖3、圖4所示。

圖3 復(fù)折射因子的實(shí)部

圖4 復(fù)折射因子的虛部

圖3、圖4表明，濾波前、窗函數(shù)濾波和虛部增量法3種情況下的復(fù)折射因子，在需要求解區(qū)域內(nèi)復(fù)折射因子是完全重合的，在需要求解區(qū)域上部添加的吸收層內(nèi)，窗函數(shù)法和虛部增量法2種方法所得得到的復(fù)折射因子完全一致，均平滑衰減為0。

為了驗(yàn)證應(yīng)用虛部增量法進(jìn)行PE模型的步進(jìn)求解的效果，取10km 處場(chǎng)的垂直剖面如圖5所示。

圖5 傳播距離10km 處的場(chǎng)強(qiáng)幅度

如圖5 所示，在垂直方向上，垂直高度3／4 以上（對(duì)應(yīng)于1024點(diǎn)FFT，就是768～1024 點(diǎn)區(qū)域）的吸收層內(nèi)，電波波振面上場(chǎng)強(qiáng)幅度平滑衰減的趨勢(shì)明顯，并最終衰減為0，這是采用虛部增量法后，負(fù)折射因子在吸收層內(nèi)平滑衰減對(duì)電波的作用結(jié)果。必須要說明的是此處所謂 “平滑衰減”是指電波在原有信號(hào)幅度基礎(chǔ)上的平滑衰減，是信號(hào)包絡(luò)的平滑衰減。這表明虛部增量法可使場(chǎng)在吸收層內(nèi)平滑的衰減為0，從而滿足了Sommerfeld輻射條件。

因此，仿真表明虛部增量法能達(dá)到與窗函數(shù)法相同的吸收效果，并降低了計(jì)算量需求。

5 結(jié)束語

針對(duì)拋物方程模型的分步傅里葉變換法求解電波傳播時(shí)，上邊界吸收層窗函數(shù)濾波法在每一步進(jìn)上都需要進(jìn)行一次濾波的問題，提出了與窗函數(shù)法等效的修正折射率虛部增量法，理論推導(dǎo)了這兩種方法的等價(jià)關(guān)系，設(shè)計(jì)了計(jì)算流程，并對(duì)2種方法的等效性進(jìn)行了仿真分析。通過將窗函數(shù)濾波等價(jià)為修正折射率的虛部增量，修正折射率虛部增量法節(jié)約了計(jì)算量，在多折射率剖面插值和長(zhǎng)距離電波傳播計(jì)算的情況下，較有意義。雖然本文的方法是基于Cosine－tapered窗函數(shù)推導(dǎo)的，但同樣適用于其它類型窗函數(shù)情況。后續(xù)研究可以尋找更為高效的滿足邊界條件的窗函數(shù)，并利用本文的等價(jià)關(guān)系，得到等價(jià)的折射指數(shù)虛部增量，應(yīng)用于PE的步進(jìn)求解。

［1］WANG Hongguang，ZHANG Rui，KANG Shifeng，et al.Overview on parabolic equation model research for atmospheric ducting propagation ［J］.Equipment Environmental Engineering，2008，5 （1）：11－15 （in Chinese）.［王紅光，張蕊，康士峰，等.大氣波導(dǎo)傳播的拋物方程模型研究綜述 ［J］.裝備環(huán)境工程，2008，5 （1）：11－15.］

［2］Sprague R A，Patterson W L，Barrios A E.Advanced propagation model（APM）version 2.1.04computer software configuration item（CSCI）documents［R］.Space and Naval Warfare Systems Command San Diego CA，2007.

［3］Sirkova I.Brief review on PE method application to propagation channel modeling in sea environment ［J］.Central European Journal of Engineering，2012，2 （1）：19－38.

［4］YAO Jingshun，YANG Shixing.A terrain parabolic equation model for propagation over the ocean ［J］.Chinese Journal of Radio Science，2009，24 （3）：493－497 （in Chinese）. ［姚景順，楊世興.拋物方程模型在海上電波傳播中的應(yīng)用 ［J］.電波科學(xué)學(xué)報(bào)，2009，24 （3）：493－497.］

［5］LIU Aiguo，CHA Hao.Experiment study of electromagnetic wave propagation loss in oceanic evaporation duct［J］.Chinese Journal of Radio Science，2008，23 （6）：1119－1203 （in Chinese）.［劉愛國(guó)，察豪.海上蒸發(fā)波導(dǎo)條件下電磁波傳播損耗實(shí)驗(yàn)研究 ［J］.電波科學(xué)學(xué)報(bào)，2008，23 （6）：1119－1203.］

［6］YANG Chao.Critical technologies of wave propagation in atmospheric duct and its inversion ［D］.Xi’an：Xidian University，2010 （in Chinese）. ［楊超.大氣波導(dǎo)中電磁波傳播及反演關(guān)鍵技術(shù) ［D］.西安：西安電子科技大學(xué)，2010.］

［7］LI Dexin，YANG Rijie，GUAN Wei，et al.Research on two－way parabolic equation modeling under irregular terrain environment［J］.Journal of Astronautics，2012，33（2）：235－240（in Chinese）.［李德鑫，楊日杰，管巍，等.不規(guī)則地形條件下的雙向拋物方程模型研究［J］.宇航學(xué)報(bào)，2012，33（2）：235－240.］

［8］Apaydin G，Sevgi L.Numerical investigations of and path loss predictions for surface wave propagation over sea paths including hilly island transitions［J］.IEEE Trans on Antenna and Propagation，2010，58 （4）：1302－1314.

［9］Anwar Z.Efficient absorbing boundary conditions for modeling wave propagation ［D］.North Carolina State University，2005：1－3.

［10］Dockery GD，Awadallah RS，F(xiàn)reund DE，et al.An overview of recent advances for the TEMPER radar propagation model［C］／／Proceedings of IEEE Radar Conference，2007：896－905.