彎折譜估計在拋物型方程模型電波傳播中的應(yīng)用研究

肖金光 周新力 劉曉娣

（海軍航空工程學(xué)院電子信息工程系，山東 煙臺264001）

引 言

拋物型方程（Parabolic Equation，PE）模型的分步傅里葉（Split Step Fourier Transform，SSFT）解法能夠步進求解，且步長可靈活選取，適于計算遠(yuǎn)距離電波傳播，受到了國內(nèi)外的廣泛關(guān)注［1-4］.美國“高級折射效應(yīng)預(yù)測系統(tǒng)”（Advanced Refractive Effects Prediction System，AREPS）的計算核心高級傳播模型（Advanced Propagation Model，APM）就是基于PE研制的，類似軟件還有對流層拋物方程程序（Tropospheric Electromagnetic Parabolic Equation Routine，TEMPER）［3，5］.海上發(fā)生概率較高的蒸發(fā)波導(dǎo)中，電波頻率和輻射角滿足一定條件時可實現(xiàn)超視距傳播，這種傳播機制對低高度低仰角對海輻射源影響較大，具有重要的研究價值［6］.在海風(fēng)作用下，海面表面相對于可陷獲傳播波段的電波波長呈光滑或粗糙表面特性.PE通過邊界阻抗引入下邊界粗糙海面或者地形的影響，而邊界阻抗取值依賴于電波在表面上的入射余角.常用的入射余角解法有射線追蹤（Ray Tracing，RT）和平面波譜估計（Plane Wave Spectral Estimation，PWS）兩種.APM和TEMPER在蒸發(fā)波導(dǎo)條件下采用RT［7］，國內(nèi)研究PE模型時，對入射余角延續(xù)了類似的思路［8］.RT法所得入射余角在部分區(qū)域可能多值或無值造成入射余角取值震蕩或不得不直接忽略表面粗糙的問題，而大氣的超折射作用使得PWS的平面波前假設(shè)不成立，估計誤差非常大［7，9］.為了計算海雜波，2012年Ali Karimian等提出了彎折波譜估計（Curved Wave Spectral Estimation，CWS），通過將接收陣列的陣列流型與大氣折射聯(lián)系起來，克服了PWS對平面波前的要求，利用PE計算結(jié)果估計蒸發(fā)波導(dǎo)條件下的表面入射余角［10］.但是，進行阻抗邊界PE模型計算需要用到入射余角，對當(dāng)前步進的PE結(jié)果進行入射余角估計，就會成了一個“雞生蛋、蛋生雞”問題.

為解決上述問題，本文提出單獨運行一次在光滑平面條件下的PE，用CWS估計所需的入射余角，然后用于粗糙海面的PE計算，并通過數(shù)值試驗驗證了方法的可行性.

1 阻抗邊界條件下的拋物方程模型

1.1 拋物方程模型

笛卡爾坐標(biāo)系下，時諧因子為exp（-i wt），二維標(biāo)量波方程忽略后向傳播并作近軸近似，可得標(biāo)準(zhǔn)拋物方程（Standard Parabolic Equation，SPE）為［1］

式中：k為電磁波波數(shù)；n（x，z）為大氣折射率；u（x，z）為場分量.考慮地球曲率的影響（re為地球半徑），SPE的SSFT解為［11］

顯然式（2）可以步進迭代解算，兩指數(shù)項分別表征了傳播媒質(zhì)的折射效應(yīng)和障礙物繞射效應(yīng).F、F-1為正、反傅里葉變換，變換點數(shù)N由奈奎斯特準(zhǔn)則zmaxpmax＝Nπ確定，zmax和pmax分別為需輸出的高度和變換域最大值.若β是以輻射源點為頂點的最大計算域錐角，相鄰兩點高度差Δz需滿足［1，11］

初始場由天線方向圖的傅里葉逆變換求得，上邊界采用Tukey窗函數(shù)濾波，在需要輸出計算值區(qū)域上方添加吸收層，層內(nèi)用窗函數(shù)濾波使場緩慢衰減為0，避免邊界上的強反射“污染”計算域.下邊界為光滑表面時，可由鏡像原理將FFT進一步簡化為正弦變換或余弦變換，但對于粗糙表面或地形等阻抗邊界，需要引入滿足SPE解的輔助函數(shù)

式中α為邊界阻抗.根據(jù)Dirichlet邊界條件，z＝0時w＝0，場分量為u的電磁波在阻抗邊界上的傳播就等效于為w在光滑阻抗邊界上的傳播，可得混合傅里葉變換（Mixed Fourier Transform，MFT）解.若引入式（4）的差分形式，同樣可以步進求解PE，這種方法稱為離散混合傅里葉變換（Discrete Mixed Fourier Transform，DMFT）法［8］.

傳播因子Fa定義為接收點場強E和自由空間接收點場強E0之比，由PE計算結(jié)果可求得［3］

則傳輸損耗L（x，z）為

1.2 邊界阻抗與入射余角

阻抗邊界條件下，式（2）的求解與邊界阻抗α密切相關(guān).阻抗邊界條件可描述為［8］

粗糙海面條件下有［8，12］

式中：Γ為電波在海面上的反射系數(shù)；θ為其入射余角；ρ為粗糙度衰減因子.若復(fù)介電常數(shù)為ξ′，則

式中：HPol和VPol分別表示水平極化和垂直極化.根據(jù)米勒-布朗模型，粗糙度衰減因子ρ為［13］

式中I0是第一類零階修正貝塞爾函數(shù).瑞利粗糙參數(shù)γ＝2khrmssinθ，表示電波以入射余角θ入射到均方根高度為hrms的海面上時，入射波與反射波的相位差，由Philips海譜，風(fēng)速為w時hrms＝0.005 1 w2.γ趨于0時ρ趨于l，則光滑阻抗邊界條件下α為

由于能夠在蒸發(fā)波導(dǎo)中陷獲傳播的電波入射余角一般都小于1°，此時α又可以進一步近似為

由式（7）～（11）可知，α是θ的函數(shù)，因此入射余角θ對粗糙海面條件下的PE模型求解至關(guān)重要.

2 入射余角求解

2.1 RT和PWS譜估計法

RT本質(zhì)上是幾何光學(xué)近似，其泰勒級數(shù)近似解法較適于小仰角電波傳播情況.根據(jù)折射率梯度將大氣分層，計算射線軌跡［14］，與電波的頻率無關(guān).

對大氣折射率n引入平坦地球近似，以將電波沿地球表面的傳播等效為水平傳播，則平坦地球修正大氣折射率為

大氣修正折射指數(shù)M為（nmod-1）×106.在蒸發(fā)波導(dǎo)中有［3］

式中：表面粗糙參數(shù)z0為1.5×10-4m；d為波導(dǎo)高度；M（z0）取為330M單位.取蒸發(fā)波導(dǎo)高度24m，天線高度20m，基于RT的射線軌跡如圖1所示.

由圖1可以看出：RT法用射線近似電波傳播路徑，在部分距離上沒有射線到達(dá)表面，這種部分區(qū)域的無值現(xiàn)象在天線高度較高時比較明顯，天線高度較低時會有所改善；而部分地方可能存在兩個甚至多個值，這種多值性在天線高度較低時會更加突出.取24m高蒸發(fā)波導(dǎo)中天線高度為10m，RT估計的入射余角如圖2所示.

圖1 20m天線的射線軌跡

圖2 10m天線RT估計的入射余角

由圖2可知入射余角在部分距離上存在可能多值，這在天線安裝位置較低時更加明顯.這種在相鄰區(qū)域甚至是同一距離上，因多條射線到達(dá)而形成的入射余角多值，恰恰表明了RT不能反映多路徑傳播電波的疊加效應(yīng)，不能表征衍射、繞射機制，無射線到達(dá)區(qū)域的衍射、繞射和反射直接被忽略了.對于多值的取舍及造成取值震蕩，文獻［15］取RT估值最大值包絡(luò)作為最終擬合值.

時域信號的頻譜表示信號在頻率上的能量分布，而空間譜表示的是信號在空間上的能量分布，從空間譜可得信號波達(dá)方向（Dirrect of Arrive，DOA）.PWS假設(shè)電波具有平面波前，取空域上的部分場值（可以通過補零增加序列長度），經(jīng)傅里葉變換后就是入射波的空間譜，搜索譜峰值即得到入射波主能量的DOA［16］.APM中對于非蒸發(fā)波導(dǎo)情況下粗糙表面的DOA估計，通過單獨運行一次光滑表面條件下的PE，取下表面附近的部分?jǐn)?shù)據(jù)，估計到達(dá)下表面的電波的入射余角，另外在混合模式下，PE計算區(qū)域電波溢出方向的估計也采用PWS.但是蒸發(fā)波導(dǎo)條件下，大氣超折射作用使得PWS的平面波前假設(shè)不成立，致使估計誤差非常大.

2.2 CWS譜估計法［10］

高度z處的電波傳播速度為

式中c0為真空電波傳播速度.以θz入射到高度z處電波的波數(shù)為

由費涅耳定理，水平波數(shù)為

式中θ是海表面z1處的入射余角.則垂直波數(shù)為

垂直方向上場的相位變化就是垂直波數(shù)的積分.取貼近表面的Na點PE計算復(fù)場值進行譜估計，則第l點相對海面的相位差為

則空間譜為

對于線性陣列，在-90°～90°范圍內(nèi)搜索式（18）的峰值，空間譜取得峰值時，對應(yīng)角度θ即為入射波在參考高度上的入射余角，CWS通過式（12）～（18）將接收陣元相對參考高度的相位變化與大氣折射率剖面和電波波長聯(lián)系起來，從而可以應(yīng)用于非平面波前的情況［10］.而PWS實質(zhì)上是CWS中垂直波數(shù)為常數(shù)的特例，第l點相對海面的相位差為

可得PWS空間譜為

由式（18）可知，CWS中，相同大氣折射率條件下，電波頻率越高，垂直波數(shù)越大，垂直接收陣列上場的相位變化越大，這一與電波波長相關(guān)聯(lián)的特性是RT所不具備的.

2.3 三種方法估計特性比較

對RT、PWS、CWS三種入射余角方法的估計特性進行列表比較，如表1所示.

由表1可知CWS不僅表征了大氣非均勻特性，而且與電波波長有關(guān)，估值具有唯一性，規(guī)避了RT和PWS的缺陷，因此用于海面電波傳播的入射余角估計更具優(yōu)勢.

表1 RT、PWS、CWS估計特性比較

3 CWS用于PE計算

能否在一次計算中直接應(yīng)用CWS呢？如果采用文獻［10］的思路，直接對PE結(jié)果進行入射余角估計，就會在計算當(dāng)前距離場值時需要基于該場值計算結(jié)果估計的參數(shù)，這顯然是一個“雞生蛋、蛋生雞”的悖論.

針對這一矛盾，本文提出單獨運行一次在光滑平面條件下的PE，用CWS估計所需距離上的入射余角，然后用于粗糙海面的PE計算，步驟如下：

1）PE模型相關(guān)參數(shù)設(shè)置；

2）根據(jù)波導(dǎo)剖面，用CWS估計各陣元相對于參考高度的相位差；

3）運行光滑表面條件下的PE模型；

4）根據(jù)3）的結(jié)果，用2）所得陣列流型，搜索譜峰值，對應(yīng)的角度就是電波在該距離處海表面上的入射余角；

5）基于所估計的入射余角進行粗糙海面的PE模型計算.

算法可行性分析：

1）TEMPER對存在地形的情況也是單獨運行一次在光滑平面條件下的PE，以進行PWS，在地形處取RT和SE的最大值，因為大氣超折射作用使得PWS的平面波前假設(shè)不成立，導(dǎo)致誤差太大而無法將PWS應(yīng)用于蒸發(fā)波導(dǎo)情況［7］.而CWS實質(zhì)上是通過陣列流型的相位補償，將入射波等效為平面波了，因此從原理上具有可行性.

2）由于CWS對光滑表面的PE計算場值進行估計，因此可以估計任意距離上入射余角，且估計值具有唯一性，彌補了RT部分區(qū)域存在多值或無值的兩個重要缺陷.

3）避免了計算當(dāng)前距離場值時需要基于該場值計算結(jié)果估計的參數(shù)這一悖論.

4）光滑平面條件下的PE，不需要考慮有限邊界阻抗的問題，且初始化參數(shù)設(shè)置與第3步相同.運行光滑海面條件下的PE，計算量相對于RT要大，如果存在地形情況，相對于TEMPER所采取的在RT和PWS兩種方法所得結(jié)果做選擇，因為省略了RT，計算量反而是少的.因而從計算量的角度上具有應(yīng)用的可行性.

5）蒸發(fā)波導(dǎo)陷獲傳播電波的入射余角一般小于1°，因此在超視距傳播區(qū)域的入射余角估計中可以將譜峰搜索的范圍限定在一個較小的半空間角度內(nèi)，如0～1°，可以進一步降低計算量.

4 數(shù)值試驗分析

4.1 參數(shù)設(shè)置

設(shè)輻射源頻率為10GHz，高斯天線高20m，仰角為0°，蒸發(fā)波導(dǎo)高24m，粗糙海面時取風(fēng)速為5 m／s，需輸出高度為0～100m.根據(jù)式（3）陣元間距為0.16m，需1 024點FFT，步長取為Δr＝4π（Δz）2／λ即296.5m.陣元數(shù)選為［10］

式中θBW是第一零點波瓣寬度（rad），對應(yīng)角度值取為0.17°，則陣元數(shù)為126.

4.2 基于光滑海面PE的CWS性能分析

由式（17）、（19）、（21）及設(shè)定的參數(shù)，取角度θ搜索范圍中一點0.47°為例，PWS和CWS各陣元相對參考海面的相位差φl（0.47）如圖3所示.

圖3 PWS和CWS相對參考面的相位差

由圖3可知，由于蒸發(fā)波導(dǎo)中的大氣陷獲折射作用，各陣元間的相位差不再均勻，垂直接收陣列的底部相位梯度大，而頂部梯度小，這與式（14）所示大氣折射指數(shù)剖面曲線特性一致.CWS估計所用的陣元相位陣列表征了這一作用機制，是其相對于PWS的本質(zhì)區(qū)別.

對光滑海面條件下的PE計算場值作PWS和CWS估計，并與RT所估計的入射余角作比較，如圖4所示.

圖4 光滑海面入射余角估計

由圖4可知，CWS與RT的入射余角估計值走勢吻合較好，CWS估值略小于RT，而PWS誤差較大.而且非常明顯的，RT法在開始階段的入射余角估計值密集，但在40～80km范圍內(nèi)和末端估計值較為稀疏，這與圖1射線到達(dá)的位置是一致的.與RT不同，PWS、CWS是對光滑表面PE計算場值進行估計，可以根據(jù)需要獲得所需步進點上的估計值.

圖4中角度估值的最大值與輻射源高度和PE步長有關(guān)，步長小、輻射源布設(shè)位置高時，就會較大，反之亦然.相應(yīng)的，在譜峰搜索時，可以根據(jù)幾何關(guān)系或者RT法進行粗略估算，將搜索范圍限定在一個適當(dāng)小的角度范圍內(nèi).

光滑海面CWS（風(fēng)速為0）、粗糙海面CWS（對風(fēng)速為10m／s的計算結(jié)果進行CWS估計）與RT三種方法超視距部分的估計值均值隨輻射源頻率的變化如圖5所示.

圖5 估計均值與頻率的關(guān)系

由圖5可知，RT估值不隨頻率和風(fēng)速變化，CWS估值比RT的估值小，電波頻率越低，兩者之間的差值越小，頻率越高差值越大.由式（18）可知，頻率越高，垂直波數(shù)受到大氣非均勻特性影響越大，而RT與電波的頻率沒有關(guān)系.CWS通過使用場值進行角度估計，表征了多徑信號在空間中的相互調(diào)制機制.比較光滑海面CWS（風(fēng)速為0）、粗糙海面CWS（風(fēng)速為10m／s）可知，風(fēng)速越大，CWS估值越大，但總體上風(fēng)速對CWS結(jié)果影響較小，這也是本文中用光滑海面PE結(jié)果代替粗糙海面PE結(jié)果（暫不可得）進行CWS估計的原因.

此外，CWS和RT估值均值隨波導(dǎo)高度和天線高度的變化呈相同變化規(guī)律，在此不再贅述.

4.3 基于CWS的粗糙海面PE性能分析

基于光滑海面條件下PE模型計算結(jié)果，進行CWS估計，估計的入射余角用于粗糙海面條件下的PE模型計算.10GHz輻射源的電波傳播損耗如圖6所示（對于高于180dB的值和低于110dB的值均歸入相近的區(qū)間表示）.對不同頻率的輻射源，分別基于RT法和CWS估計入射余角進行PE計算，在不同距離和高度上進行比較，分別如圖7～9所示.為了清晰顯示差異，圖8、9中繪出了兩者傳播損耗的差值（CWS-RT）.

圖6 基于CWS入射余角的電波傳播損耗

由圖7～9可知，10GHz信號源基于CWS的傳播損耗比RT大0.1dB，走勢一致（圖7），而20 GHz時兩者差值最大可達(dá)數(shù)dB（圖8、9），且隨距離增大而增大，在較低高度上比較明顯.這與式（18）的內(nèi)涵一致，較低頻率電波的垂直波數(shù)變化較小，電波受大氣非均勻特性的影響也較小，當(dāng)頻率變大，垂直波數(shù)變大，電波受大氣非均勻特性的影響也變大.同樣是表征大氣非均勻特性，CWS與電波波長有關(guān)，而RT與電波波長無關(guān).20GHz頻率下，CWS的損耗總體上比RT的大，接收高度越高，差值曲線越趨于平緩；天線較低時低高度上的差值較大.這是因為天線高度較低時，RT估值存在多值性和跳躍性，不能表征多徑信號疊加的缺陷更加明顯.

圖7 20m高10GHz信號源的傳播損耗

圖8 20m高20GHz信號源的傳播損耗

5 結(jié) 論

圖9 5m高20GHz信號源的傳播損耗

入射余角是拋物方程法求解阻抗邊界電波傳播中的重要參數(shù).海面蒸發(fā)波導(dǎo)情況下，目前采用的RT法估計入射余角存在部分區(qū)域存在多值或無值以及與電波頻率無關(guān)的缺陷，而PWS因為平面波前假設(shè)不成立而誤差較大.CWS表征了大氣非均勻特性對電波波前的影響，同時避免了計算當(dāng)前值需要用到基于當(dāng)前值的估計參數(shù)這一個悖論，本文提出了一種基于CWS入射余角估計進行PE模型海面電波傳播計算的方法.通過對光滑海面的PE計算結(jié)果進行CWS譜估計，將估計的入射余角用于粗糙海面的PE模型計算，仿真和分析表明基于CWS的PE計算結(jié)果與基于RT的走勢契合較好，較低頻率電波的垂直波數(shù)變化較小，電波受到大氣非均勻特性的影響也較小，兩者差值較??；在較高頻率時，CWS計算的低空傳播損耗明顯較大，更好地反映了海面低空大氣超折射對較高頻率電波的顯著作用.新方法可以估計任意距離上入射余角，且估計值具有唯一性，規(guī)避了RT在部分區(qū)域多值或無值而相應(yīng)造成的取值震蕩，及其與電波頻率無關(guān)的缺陷.新方法的計算量相比于純RT法要大，相比地形條件下的RT和PWS選擇法要小.新方法為蒸發(fā)波導(dǎo)和地形并存條件下，一致地求解粗糙面入射余角提供了解決方法，后續(xù)將予以研究.

［1］LEVY M F.Parabolic Equation Methods for Electromagnetic Wave Eropagation［M］.London：IEE Press，2000：20-35.

［2］王紅光，張 蕊，康士峰，等.大氣波導(dǎo)傳播的拋物方程模型研究綜述［J］.裝備環(huán)境工程，2008，5（1）：11-15.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）

［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］BARRIOS A E.Considerations in the development of the advanced propagation model（APM）for U.S.Navy applications［C］／／Proceedings of the International Conference on Radar RADAR’2003Radar Conference，2003：77-82.

［5］BROOKNER E，CORNELY P R，LOK Y F.AREPS and TEMPER-getting familiar with these powerful propagation software tools［J］.2007.

［6］黃小毛，張永剛，王 華，等.蒸發(fā)波導(dǎo)中電磁波異常傳播特征研究及其應(yīng)用［J］.電子與信息學(xué)報，2006，28（8）：1508-1512.HUANG Xiaomao，ZHANG Yonggang，WANG Hua，et al.A ray tracing algorithm for duct environment［J］.Journal of Electronics ＆Information Technology，2006，28（8）：1508-1512.（in Chinese）

［7］DOCKERY G D.An overview of recent advances for the TEMPER radar propagation model［J］.IEEE Transactions on Antennas and Propagation，2007，44（1）：896-905.

［8］郭立新，李宏強，楊 超，等.改進DMFT算法研究粗糙海上蒸發(fā)波導(dǎo)中的電波傳輸特性［J］.電波科學(xué)學(xué)報，2009，24（3）：414-421.GUO Lixin，LI Hongqiang，YANG Chao，et al.Characteristics of radio wave propagation in the evaporation duct environment over the rough surface by the improved DMFT algorithm［J］.Chinese Journal of Radio Science，2009，24（3）：414-421.（in Chinese）

［9］SPRAGUE R A，PATTERSON W L，BARRIOS A E.Advanced Propagation Model（APM）Version 2.1.04 Computer Software Configuration Item（CSCI）Cocuments［R］.San Diego：Space and Naval Warfare Systems Command，2007.

［10］KARIMIAN A，YARDIM C，GERSTOFT P，et al.Multiple grazing angle sea clutter modeling［J］.IEEE Transactions on Antennas and Propagation，2012，10（9）：4408-4417.

［11］姚景順，楊世興.拋物方程模型在海上電波傳播中的應(yīng)用［J］.電波科學(xué)學(xué)報，2009，24（3）：493-497.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）

［12］劉愛國，察 豪.海上蒸發(fā)波導(dǎo)條件下電磁波傳播損耗實驗研究［J］.電波科學(xué)學(xué)報，2008，23（6）1119-1203.LIU Aigou，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）

［13］FREUND D E，WOODS N E.Forward radar propagation over a rough sea surface：an numerical assessment of the Miller-Brown approximation using a Horizonally polarized 3GHz line source［J］.IEEE Trans on Antennas Propagat，2006，54（4）：1292-1304.

［14］孫 方，王紅光，康士峰，等.大氣波導(dǎo)環(huán)境下的射線追蹤算法［J］.電波科學(xué)學(xué)報，2008，23（1）：179-183.SUN Fang，WANG Hongguang，KANG Shifeng，et al.A ray tracing algorithm for duct environment［J］.Chinese Journal of Radio Science，2008，23（1）：179-183.（in Chinese）

［15］DOCKERY G 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.

［16］湯 俊.統(tǒng)計信號處理［M］.北京：清華大學(xué)出版社，2006.