任 明，錢志華

(杭州電子科技大學天線與微波技術(shù)研究所，杭州310018)

電波傳播一直是工程電磁場理論和環(huán)境電磁特性研究領(lǐng)域中最為人們廣泛關(guān)注和研究的方向之一。拋物型方程法是從Maxwell 方程組出發(fā)做前向傳播近似，從而把橢圓型的波動方程簡化為拋物型方程［1］。拋物型方程方法在精確預(yù)測電波傳播路徑損耗中得到了廣泛的應(yīng)用，使得電波傳播數(shù)值計算有了新的突破。本文采用有限差分算法，可以直接計算出任意邊界上的場，因此在處理不規(guī)則地形邊界時比較方便，同時解決了大規(guī)模矩陣求解難的問題，為計算復(fù)雜環(huán)境下大尺度寬范圍的電波傳播損耗提供了一種有效方法。

1 拋物型方程的Crank－Nicolson 離散

本文設(shè)時諧場因子為e－jωt，電磁場分量在無源空間滿足波動方程［1］:

式中，k0，n 分別為自由空間波數(shù)和大氣折射率。設(shè)想電磁波沿主軸(x 軸)方向傳播，提取沿x 軸方向相位發(fā)生快速變化的這一項，即令u=e－jk0xψ，

考慮到u 關(guān)于x 的二階導(dǎo)數(shù)很小，忽略此項，式(2)就轉(zhuǎn)換成標準的三維拋物型方程:

為了數(shù)值求解式(3)，我們采用Crank－Nicolson(CN)差分方案［2］。因此:

上標k 和下標m，n 分別是x，z 和y 方向上的位置示意，Δx，Δy，Δz 分別是x，y 和z 方向上的空間步長。代入式(3)，將關(guān)于k+1 和k 的項分別置于方程的兩邊，得到

不難發(fā)現(xiàn)，k+1 平面內(nèi)的場可以通過k 平面的場遞推得到。如果已知初始場，就可以借助相鄰平面間場值的遞推求出任意感興趣平面上的場分布。定義UM×N為表征二維場點的矩陣，則式(4)可以表示成一般矩陣形式:

其中A，B，D，E 分別是如下所示的三對角方陣:

需要特別強調(diào)的是，式(6)只適用于完全位于離散區(qū)域內(nèi)的場點。如圖1 所示，意味著m 不能?。? 或M，n 不能取0 或N。因此，總的待求場點數(shù)實為M×(N－1)。為建立完整的求解方程組，還需結(jié)合初始場分布，以及地表邊界條件、吸收邊界條件經(jīng)離散得到輔助方程。幸運的是，用有限差分法處理上述條件非常方便，詳見后續(xù)內(nèi)容。

圖1 二維場點的位置示意圖

2 激勵源

遞推所需初始場由激勵源產(chǎn)生，本文采用了高斯電流源。假設(shè)在初始距離處，電磁場分量由一個電流源J 輻射產(chǎn)生［3］:

式中I0l 為電流矩(取1)，ht是發(fā)射天線的架設(shè)高度，σz決定著高斯方向圖在z 方向的3 dB 波束寬度。δ(y－y0)是狄拉克函數(shù)，意味著天線的位置在y軸上的投影為y=y0。

3 吸收邊界條件

3.1 列昂托維奇(Leontovich)地面邊界條件

三維場的數(shù)值計算需要耗費巨大的計算機資源。因此，有必要對計算區(qū)域進行一定的限制。為此，在計算區(qū)域的截斷邊界處必須提供吸收邊界條件。

圖1 中處于m=0 的場點，應(yīng)滿足一定的空/地邊界條件。電波傳播區(qū)域的下邊界通常并非理想導(dǎo)電平面，而是滿足阻抗邊界條件，又稱列昂托維奇邊界條件［4］。即

其中η 是歸一化表面阻抗，它由地表面的相對介電常數(shù)和導(dǎo)電率決定，還取決電磁場分量的極化方式，利用地面下方的虛擬場點(見圖1)以及中心差分思想，可將式(8)離散為:

代入式(6)，消去虛擬場點項，得到波動方程融合地表作用后的差分方程:

3.2 非局部吸收邊界條件(NLBC)

這里介紹3DPE 方程離散算法常用的非局部吸收邊界條件［2］。這種邊界條件最大的好處是其占用輔助場點數(shù)較少，節(jié)約計算機存儲單元。

不妨以計算區(qū)域頂部吸收邊界條件為例，左右吸收邊界條件的推導(dǎo)只需簡單地交換坐標變量。借助拉普拉斯變換，嚴格按照文獻［2］中的思路，可以得到:

在大氣折射率n=1 的條件下，取φ1=22°，φ2=45°，NLBC 對電磁波的吸收效果最好。不難發(fā)現(xiàn)，式(12)涉及到卷積計算，其物理意義就是:在當前距離處，邊界條件的滿足依賴于所有以前距離處的場。

式(12)的離散借助輔助場點(圖1 中標有“○”的點)，在相鄰遞推平面區(qū)間作卷積運算時將二階導(dǎo)數(shù)項簡單看作常數(shù)，即不隨距離發(fā)生變化，并取作:

將式(14)、式(15)代入式(6)，最終得到頂部NLBC的差分方程:

遵循同樣的思路，可以分別推導(dǎo)出左右邊界處NLBC 滿足的方程，注意左邊界與右邊界的方程之間差一個負號。需特別指出的是，在左右邊界NLBC 方程中，常數(shù)將變?yōu)镃3、C4(只需要將正弦改為余弦即可)。到此，圖1 中標有“* ”的點滿足的差分方程就全部推導(dǎo)完畢。由于場點分布呈二維矩形，圖1 中“Δ”的4 個角點，它們對應(yīng)離散方程聯(lián)合地面條件和NLBC 單獨推導(dǎo)即可。

4 矩陣預(yù)處理與快速求解技術(shù)

考察式(5)，其中U 為表征二維場點的M×N 矩陣，A，B，D，E 均為三對角系數(shù)方陣。若已知初始場，k+1 平面內(nèi)的場將通過k 平面的場遞推得到。求解這種類型的矩陣方程，傳統(tǒng)的做法是先將方程轉(zhuǎn)換成新的矩陣方程Tx=b(這里T 是M×N 階方陣)，然后運用迭代法。顯然，傳統(tǒng)做法對計算機的存儲量要求很高，求解速度慢，不能適應(yīng)大范圍、遠距離電波傳播特性預(yù)測的要求。為克服這一缺點，本文考慮先將矩陣進行Schur 分解，再實現(xiàn)快速求解。

對于式(5)，令DU+UE=C，則

這種方程也稱作Sylvester 方程。矩陣預(yù)處理技術(shù)如下，先將矩陣A，B 進行Schur 分解:

式中W，Y 為酉矩陣，即W*W=WW*=I，Y*Y=YY*=I。RA和RB均為上三角矩陣，其主對角元素分別是矩陣A 和B 的特征值，特征值對應(yīng)的特征向量分別構(gòu)成W，Y 的列。

將式(18)代入式(17)，得:WRAW*U+UYRBY*=C，然后在方程兩邊分別左乘矩陣W*，右乘矩陣Y，同時利用酉矩陣的性質(zhì)，有:RAW*UY+W*UYRB=W*CY

令W*UY=，W*CY=，于是方程轉(zhuǎn)化為

此時問題轉(zhuǎn)化為已知矩陣^C，求解矩陣^U。顯然，式(19)的求解非常簡單且快速，只需直接回代，因為矩陣RA和RB均是上三角矩陣。一旦求得^U，則最初的待求矩陣為:

5 算例與討論

為檢驗我們的算法，以電波在平靜海平面上的傳播為例，計算電波傳播損耗值。并將結(jié)果與雙射線法［5］比較。假設(shè)發(fā)射天線高度為ht=5 m，信號工作頻率f=1 GHz，接收天線與發(fā)射天線之間水平距離為d=1 000 m，海水電參數(shù)為:εr=80，σ=4 S/m。取空間步長為Δx=4λ，Δy=3λ，Δz=2λ，場矩陣U 在z，y 方向上的點數(shù)分別為M=120，N=30，大約迭代835 步。

圖2、圖3 分別展示了CNFD 算法結(jié)合NLBC 邊界計算垂直極化和水平極化下的電波傳播損耗隨接收天線高度變化的結(jié)果，CNFD 算法與雙射線法吻合得很好，但隨著接收天線高度的增加，CNFD 結(jié)果的出現(xiàn)不穩(wěn)定但總體上還是比較吻合，驗證了算法的正確性。

圖2 垂直極化下二種方法的比較

圖3 水平極化下二種方法的比較

最后，我們將矩陣預(yù)處理技術(shù)與傳統(tǒng)的迭代法求解技術(shù)進行了比較研究，發(fā)現(xiàn):以CPU 主頻為2 GHz 的單機為例，同樣求解130×30 個場點(其他參數(shù)與圖2、圖3 中的算例完全相同)，采用迭代法需要8 h(收斂精度取10－8)，而借助矩陣預(yù)處理技術(shù)，只需2 min。顯然，矩陣經(jīng)過預(yù)處理后求解速度大大提高。另一方面，先將矩陣進行Schur 分解，還可以大大改善計算機存儲量的限制情況。以單機2 G 內(nèi)存為例，采用迭代法大約可以求解4 000個場點，而采用矩陣預(yù)處理技術(shù)，可以求解約32 000個場點，提高了整整7 倍，這為預(yù)測大尺度寬范圍的電波傳播特性提供了可能。

6 結(jié)論

本文研究了用3DPE 方法預(yù)測電波傳播的衰減問題，結(jié)合CNFD 技術(shù)，結(jié)果與雙射線法結(jié)果吻合很好。而且矩陣預(yù)處理技術(shù)的運用是非常成功的，不但大大改善了計算機存儲量的限制狀況，同時實現(xiàn)了場值的快速遞推，為后續(xù)的程序調(diào)試和研究工作節(jié)省了大量時間。

［1］ Levy M F. Parabolic Equationmethods for Electromagnetic Wave Propagation［M］.London:IEE Press，2000:6－45.

［2］ Zelley C A，Constantinou C C. A Three－Dimensional Parabolic Equation Applied to VHF/UHF Propagation over Irregular Terrain［J］. IEEE Transactions on Antennas ans Propagation，1999，47(10):1586－1596.

［3］ 胡繪斌. 預(yù)測復(fù)雜環(huán)境下電波傳播特性的算法研究［D］. 長沙:國防科學技術(shù)大學，2006.

［4］ 胡繪斌，柴舜連，毛鈞杰. 寬角拋物方程在阻抗邊界條件下的應(yīng)用［J］.電波科學學報，2005，20(4):500－504.

［5］ 胡繪斌，柴舜連，毛鈞杰.基于三維PE 方法的雷達波傳播損耗估計［J］.微波學報，2005，21(2):4－7.

［6］ Hoffman J D. Numerical Methods for Engineers and Scientists［M］.New York:McGraw－Hill，1992.

［7］ Taflove A，Hagness S C. Computational Electrodynamics:The Finite－Difference Time－Domain Method［M］. 2nd ed. Artech House，Inc.，2000.

［8］ 李永順.基于拋物線方程求解電大尺寸目標雙(多)站RCS 方法的研究［D］.合肥:安徽大學，2005.

［9］ 李方，察豪.拋物線方程法研究不規(guī)則地形對電波傳播的影響［J］.火力與指揮控制，2010，35(8):114－116.

［10］ 孔勐，吳先良.Pade(1，0)拋物線方程方法求解三維電磁散射問題［J］.安徽大學學報:自然科學學報，2008，32(32):44－47.

［11］ 孔勐，胡玉娟，吳先良.標準拋物線方程SSFT 算法在電波傳播問題中的應(yīng)用［J］.安徽大學學報:自然科學學報，2011，35(3):2－4.

［12］ 唐海川，黃小毛，吳繩正，等.海上大氣負折射現(xiàn)象及其對電波傳播影響分析［J］.微波學報，2009，25(6):43－48.

［13］ 王文祥.微波工程技術(shù)［M］.北京:國防工業(yè)出版社，2009.

［14］ 謝益溪.無線電波傳播原理與應(yīng)用［M］. 北京:人民郵電出版社，2008:61－120.

［15］ 熊皓.無線電波傳播［M］.北京:電子工業(yè)出版社，2000:403－413.