用三維拋物方程的CNFD 技術計算傳播損耗*

任 明，錢志華

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

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

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

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

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

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

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

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

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

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

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

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

2 激勵源

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

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

3 吸收邊界條件

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

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

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

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

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

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

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

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

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

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

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

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

4 矩陣預處理與快速求解技術

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

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

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

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

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

令W*UY=，W*CY=，于是方程轉化為

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

5 算例與討論

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

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

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

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

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

6 結論

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

［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］ 胡繪斌. 預測復雜環境下電波傳播特性的算法研究［D］. 長沙:國防科學技術大學，2006.

［4］ 胡繪斌，柴舜連，毛鈞杰. 寬角拋物方程在阻抗邊界條件下的應用［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］ 李方，察豪.拋物線方程法研究不規則地形對電波傳播的影響［J］.火力與指揮控制，2010，35(8):114－116.

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

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

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

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

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

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