孟書生

（鄭州科技學(xué)院基礎(chǔ)部，河南鄭州450064）

海上大氣波導(dǎo)傳播拋物方程模型的算法研究

孟書生

（鄭州科技學(xué)院基礎(chǔ)部，河南鄭州450064）

根據(jù)電波傳播的角譜理論原理，考慮兩種極化波在導(dǎo)體、阻抗和粗糙邊界層面上具有不同的反射系數(shù)，在初始場(chǎng)和分步傅里葉變換求解的過(guò)程中加入反射系數(shù)影響，從而得到可一致計(jì)算水平或垂直極化波傳播特性的新算法.通過(guò)仿真計(jì)算電波在蒸發(fā)波導(dǎo)、粗糙海面情況下的傳播損耗，并與APM（Advanced Propagation Model）的預(yù)測(cè)結(jié)果進(jìn)行比較，驗(yàn)證了該算法的有效性.

大氣波導(dǎo)；拋物方程；角譜；邊界條件

海面存在大氣波導(dǎo)時(shí)，電波傳播損耗預(yù)測(cè)的精度對(duì)岸基或艦載雷達(dá)性能評(píng)估結(jié)果影響很大.隨著計(jì)算電磁學(xué)的發(fā)展，人們能夠通過(guò)數(shù)值求解拋物方程方法實(shí)現(xiàn)近海面復(fù)雜對(duì)流層環(huán)境中的電波傳播特性預(yù)測(cè).拋物方程是由Leontovich和Fock于1946年首先提出［1］，但是直到1973年，Hardin在海洋聲學(xué)中提出拋物方程的SSFT（Split－step Fourier Transform）算法［2］后，拋物方程才在海洋聲學(xué)、光傳播和無(wú)線電波傳播等領(lǐng)域得到了長(zhǎng)足發(fā)展.20世紀(jì)80年代末至今，Sirkova I，Mikhalev M等人廣泛地開(kāi)展了基于拋物方程的對(duì)流層電波傳播模型與算法研究［3］，使拋物方程方法成為目前解決復(fù)雜大氣和地（海）面條件下遠(yuǎn)距離電波傳播問(wèn)題最有效的手段.研究表明，平靜海面上低入射余角水平極化波傳播可按導(dǎo)體邊界條件處理，但是垂直極化波傳播和兩種極化波在海浪較大時(shí)，傳播損耗計(jì)算需考慮阻抗或海面粗糙度的影響［4］.Hitney根據(jù)波導(dǎo)模理論在粗糙海面條件下的計(jì)算結(jié)果，得到修正拋物方程模型SSFT算法的經(jīng)驗(yàn)公式［5］.Dockery等人采用DMFT（Discrete Mixed Fourier Transform）計(jì)算阻抗或粗糙海面條件下的傳播損耗［6－7］.Hannah在研究自由空間GPS信號(hào)傳播時(shí)，基于拋物方程模型仿真多徑波導(dǎo)環(huán)境下，邊界為任意導(dǎo)體、阻抗或粗糙面時(shí)可采用一致的形式和過(guò)程計(jì)算水平或垂直極化波傳播特性的新算法［8］.筆者對(duì)新算法進(jìn)行了仿真計(jì)算，并與APM（Advanced Propagation Model）的預(yù)測(cè)結(jié)果進(jìn)行了比較.

1 拋物方程模型

拋物方程是波動(dòng)方程的近似形式.取電磁場(chǎng)時(shí)諧因子為e－iwt，以標(biāo)量ψ表示電磁場(chǎng)的某一場(chǎng)分量.在空間直角坐標(biāo)系（x，y，z）下，假設(shè)ψ與y無(wú)關(guān)，即其只隨水平距離（x軸）和高度（z軸）變化，則水平極化波中只有Ey為非零電場(chǎng)分量，可設(shè)ψ＝Ey，而垂直極化波中只有Hy為非零磁場(chǎng)分量，則可設(shè)ψ＝Hy.電波傳播過(guò)程中，ψ滿足二維標(biāo)量波動(dòng)方程［9］：

式（1）中：k0為自由空間波數(shù)，n為大氣折射指數(shù).

定義沿x軸正向傳播的波函數(shù)為：

把式（2）代入式（1）并進(jìn)行不同的近似，就可以得到所謂的標(biāo)準(zhǔn)拋物方程和一種常用的寬角拋物方程形式，分別如式（3）和式（4）所示：

標(biāo)準(zhǔn)拋物方程是窄角拋物方程中的一種，在傳播仰角小于15°時(shí)具有很好的計(jì)算精度.上述形式的寬角拋物方程在計(jì)算傳播仰角小于30°時(shí)都有較好的精度，兩者都可以用SSFT算法求解.

2 SSFT算法

拋物方程的數(shù)值算法主要有三種：SSFT、有限差分法和有限元法.SSFT相對(duì)于其他算法，距離向步進(jìn)Δx可以取很大，并且利用基于FFT（Fast Fourier Transform），因此計(jì)算速度快.有限差分法和有限元法在處理地形邊界條件時(shí)具有一定優(yōu)勢(shì)，采用何種算法與具體應(yīng)用有關(guān).在大氣波導(dǎo)條件下研究近海面遠(yuǎn)距離的電波傳播問(wèn)題，出于計(jì)算效率方面的考慮，宜采用SSFT算法.

2.1 SSFT求解形式

標(biāo)準(zhǔn)拋物方程的SSFT求解形式如式（5），

式（5）中：p為變換域變量，F(xiàn)和F－1分別表示傅里葉變換和逆變換，下同.

式（4）中所表示的寬角拋物方程的SSFT求解形式如式（6），

拋物方程中傅里葉變換定義為：

與傅里葉變換把時(shí)間域函數(shù)變換到頻域類似，拋物方程中傅里葉變換是把空間函數(shù)變換到角度域，其物理意義可以用平面波角譜來(lái)解釋.與惠更斯原理不同，角譜法提供了解釋傳播過(guò)程的另一種觀點(diǎn)，它認(rèn)為在某一平面上的波場(chǎng)，可以用沿各不同方向傳播的平面波疊加而成，所有這些平面波組成該波長(zhǎng)的角譜，角譜中的每一平面波各自獨(dú)立地傳播.空間某一點(diǎn)的波場(chǎng)就是角譜中的平面波傳到該處后疊加的結(jié)果.

式（5）～（8）中p＝k0sinθ，θ為平面波分量與水平面夾角，式（8）表示某一平面上的波場(chǎng)用沿各不同方向傳播的平面波疊加而成.

2.2 初始場(chǎng)

拋物方程必須要在一定的初始條件和邊界條件下才能求解，其初始場(chǎng)一般通過(guò)對(duì)發(fā)射天線方向圖函數(shù)進(jìn)行傅里葉逆變換得到.

初始位置處場(chǎng)的垂直分布與天線方向圖f（p）為傅里葉變換對(duì)：

考慮地面效應(yīng)，設(shè)場(chǎng)的鏡像為A0（x0，z），則其傅里葉變換為f0（p），即：

若與x軸夾角為θ的平面波分量反射系數(shù)為R（p），則實(shí)際上變換域鏡像場(chǎng)則為R（p）f0（p），有：

初始場(chǎng)則為：

如果發(fā)射天線高度為z0米，仰角為θelev，那么用f（p－p0）e－ipz0代替f（p）即可，其中：p0＝k0sin（θelev）.

2.3 邊界條件

2.3.1 地表邊界

平靜海面可視為光滑阻抗表面，此時(shí)水平和垂直極化的反射系數(shù)分別為：

其中：ε＝εr＋i60λσ是復(fù)介電常數(shù)，εr為相對(duì)介電常數(shù)，σ為電導(dǎo)率.海面εr和σ值可根據(jù)CCIR建議（1986）計(jì)算獲得.

海上風(fēng)速較大時(shí)，需考慮海面粗糙度的影響，粗糙海面反射系數(shù)Rr為：

式中：RS為光滑表面的反射系數(shù)，J0為零階貝塞耳函數(shù)；σh為海面均方根高度：

u表示風(fēng)速，單位m／s.理想導(dǎo)體邊界條件下，水平極化反射系數(shù)RH＝－1，垂直極化反射系數(shù)RV＝1.

2.3.2 上邊界

由于采用電磁場(chǎng)數(shù)值計(jì)算技術(shù)，必須考慮截?cái)噙吔绲挠绊?通過(guò)加窗函數(shù)的方法處理上邊界，Cosine－taper窗函數(shù)形式如式（18）［10］：

式中：zmax為計(jì)算域的最大高度.

2.4 算法流程

與以往算法不同，筆者提出的算法除了在計(jì)算初始場(chǎng)時(shí)考慮了海面反射系數(shù)的影響，更重要的是在SSFT算法步進(jìn)求解的過(guò)程中也考慮海面反射系數(shù)因素，如式（19）和圖1所示.不同的邊界條件影響均體現(xiàn)在反射系數(shù)參數(shù)中，從而可一致地計(jì)算導(dǎo)體、阻抗和粗糙面邊界條件下的傳播損耗.

計(jì)算流程如圖1所示：

圖1 海上拋物方程算法流程

3 仿真計(jì)算與比較

仿真計(jì)算海面蒸發(fā)波導(dǎo)環(huán)境中，在不同風(fēng)速條件下，比較本算法與APM（Advanced Propagation Model）的預(yù)測(cè)結(jié)果.APM是美國(guó)的先進(jìn)對(duì)流層折射效應(yīng)預(yù)測(cè)系統(tǒng)（Advanced Refractive Effects Prediction System，AREPS）的核心計(jì)算模式，在包含大氣波導(dǎo)的低高度區(qū)域內(nèi)采用拋物方程模型，阻抗和粗糙邊界條件下APM應(yīng)用DMFT算法計(jì)算傳播損耗［11］.

仿真參數(shù)分別為：頻率8 000GHz，極化方式水平極化，蒸發(fā)波導(dǎo)高度13.5m，天線高度10m，天線仰角0°，垂直半功率波束寬度3°，高斯性天線方向圖.

圖2 10m高度處兩種算法的傳播損耗隨距離變化情況

圖3 80km距離處兩種算法的傳播損耗隨高度變化情況

圖4 160km距離處兩種算法的傳播損耗隨高度變化情況

圖2 ～圖4分別給出了不同風(fēng)速條件下傳播損耗計(jì)算結(jié)果.從圖中可見(jiàn)，海面風(fēng)速為0（光滑阻抗邊界）和10m／s（粗糙邊界）時(shí)，算法計(jì)算結(jié)果與APM模式計(jì)算結(jié)果都很一致.由圖2可見(jiàn)，固定高度10m處，不同風(fēng)速下，遠(yuǎn)距離區(qū)域內(nèi)兩種方法計(jì)算傳播損耗差值在2dB以內(nèi)；由圖3和圖4可知，以上仿真參數(shù)條件下，距離80km和160km處兩種方法得到的傳播損耗隨高度變化也很一致，相差一般在3dB以內(nèi)；蒸發(fā)波導(dǎo)高度內(nèi)，兩者結(jié)果相差小于2dB.以上結(jié)果表明了新算法的有效性.

4 結(jié)論

根據(jù)電波傳播的角譜理論原理，基于SSFT算法研究在對(duì)流層大氣波導(dǎo)環(huán)境下，在初始場(chǎng)和步進(jìn)求解的過(guò)程中考慮不同邊界面反射系數(shù)影響，得到了在大氣波導(dǎo)環(huán)境下可一致計(jì)算水平或垂直極化波傳播特性的新算法.通過(guò)仿真計(jì)算電波在蒸發(fā)波導(dǎo)、粗糙海面情況下的傳播損耗，并把本文算法與APM的預(yù)測(cè)結(jié)果進(jìn)行比較，兩者結(jié)果基本一致.另外，該算法在各種條件下都能獲得穩(wěn)定解，表明了筆者的新算法具有普適性和有效性的特點(diǎn).

［1］Leontovich M A，F(xiàn)ock V A.Solution of propagation of electromagnetic waves along the Earths’surface by the method of parabolic equations［J］.J.Phys.USSR，1946，10（1）：13－23.

［2］Hardin R H，Tappert F D.Application of the splitstep Fourier method to the numerical solution of nonlinear and variable coefficient wave equations［J］.Society for Industry and Applied Mathematics Review，1973，15（2）：423.

［3］Sirkova I，Mikhalev M.Parabolic Wave Equation Method Applied to the Tropospheric Ducting Propagation Problem：A Survey［J］.Electromagnetics，2006，26（2）：155－173.

［4］Kuttler J R，Dockery G D.Theoretical description of the parabolic approximation／Fourier split-step method of representing electromagnetic propagation in the troposphere［J］.Radio Science，1991，26（2）：381－393.

［5］Hitney H V.Modeling surface effects with the parabolic equation method［A］.Hitney H V.Surface and Atmospheric Remote Sensing［C］.San Diego：Geoscience and Remote Sensing Symposium，1994，4（1）：2322－2325.

［6］Dockery D，Kuttler J R.An improved impedance boundary algorithm for Fourier split-step solutions of the parabolic wave equation［J］.IEEE Transactions on Antennas and Propagation，1996，44（12）：1592－1599.

［7］Guillet N，F(xiàn)abbro V，Bourlier C，et al.Low grazing angle propagation above rough surface by the parabolic wave equation［A］.Guillet N.Surface and Atmospheric Remote Sensing［C］.San Diego：Geoscience and Remote Sensing Symposium，2003，7（7）：4186－4188.

［8］Hannah B M.Modeling and Simulation of GPS Multipath Propagation［D］.Brisbane：Queensland university of technology，2001.

［9］Kuttler J R，Janaswamy J.Improved Fourier transform methods for solving the parabolic wave equation［J］.Radio Science，2002，37（2）：1021－1031.

［10］Barrios A E.A terrain parabolic equation model for propagation in the troposphere［J］.IEEE Transactions on Antennas and Propagation，1994，42（1）：90－98.

［11］Barrios A E.Considerations in the development of the advanced propagation model（APM）for U.S.Navy applications［A］.Barrios A E.Proceedings of the International Radar Conference［C］.San Diego：Atmospheric Propagation Branch，2003：77－82.

Study of algorithm for atmospheric duct propagation with parabolic equation model over the sea

MENG Shusheng
（Department of Basic，Zhengzhou Institute of Science ＆Technology，Zhengzhou 450064，China）

According to the principle of the theory of angular spectrum of electromagnetic wave propagation，considering two kinds of polarized wave have different reflection coefficients on the conductor，impedance and rough surface boundary layer，adding the effect of reflection coefficients to the solving process of the initial field and the Split-step Fourier Transform solution，so the new algorithm which can be used to consistently calculate the propagation characteristics of horizontal or vertical polarization wave is obtained.The propagation loss is computed by this new algorithm in evaporation duct and rough sea surface condition.The performance of the new algorithm is then demonstrated with the result of other model.

atmospheric duct；parabolic equation；angular spectrum；boundary condition

TN011

A

1671-9476（2017）02-0066-04

10.13450／j.cnkij.zknu.2017.02.015

2016-08-24；

2016-10-18

國(guó)家高技術(shù)發(fā)展研究計(jì)劃（863計(jì)劃）“大氣波導(dǎo)實(shí)施探測(cè)技術(shù)研究”（No.2008AA093001）

孟書生（1980－），男，河南正陽(yáng)人，講師，碩士，研究方向：海洋信息探測(cè)與處理.E－mail：mss0420＠163.com，Tele：13592698724