肖金光，周新力，劉曉娣

（1.海軍航空工程學院 研究生管理大隊，山東 煙臺264001；2.海軍航空工程學院 電子信息工程系，山東 煙臺264001）

0 引 言

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

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

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

設電磁場時諧因子為e－iωt，忽略后向傳播并作近軸條件下的窄角近似可得標準拋物方程 （standard parabolic equation，SPE）為

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

其適于求解遠距電波傳播問題的SSFT 解為［3］

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

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

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

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

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

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

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

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

2 虛部增量法的理論推導

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

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

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

不妨令

則

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

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

3 虛部增量法的算法流程設計及計算量分析

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

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

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

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

步驟4 拋物方程步進求解和最高點場值置零；

步驟5 重復步驟3、步驟4，直至達到計算域距離終點。

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

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

圖1 基于虛部增量法的PE計算流程

4 仿真分析

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

圖2 虛部增量修正折射率

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

圖3 復折射因子的實部

圖4 復折射因子的虛部

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

為了驗證應用虛部增量法進行PE模型的步進求解的效果，取10km 處場的垂直剖面如圖5所示。

圖5 傳播距離10km 處的場強幅度

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

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

5 結(jié)束語

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

［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）.［王紅光，張蕊，康士峰，等.大氣波導傳播的拋物方程模型研究綜述 ［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）. ［姚景順，楊世興.拋物方程模型在海上電波傳播中的應用 ［J］.電波科學學報，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）.［劉愛國，察豪.海上蒸發(fā)波導條件下電磁波傳播損耗實驗研究 ［J］.電波科學學報，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］.西安：西安電子科技大學，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］.宇航學報，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.