任 明，錢志華

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

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

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

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

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

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

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

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

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

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

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

圖1 二維場點(diǎn)的位置示意圖

2 激勵(lì)源

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

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

3 吸收邊界條件

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

三維場的數(shù)值計(jì)算需要耗費(fèi)巨大的計(jì)算機(jī)資源。因此，有必要對(duì)計(jì)算區(qū)域進(jìn)行一定的限制。為此，在計(jì)算區(qū)域的截?cái)噙吔缣幈仨毺峁┪者吔鐥l件。

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

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

代入式(6)，消去虛擬場點(diǎn)項(xiàng)，得到波動(dòng)方程融合地表作用后的差分方程:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 算例與討論

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

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

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

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

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

6 結(jié)論

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

［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］. 長沙:國防科學(xué)技術(shù)大學(xué)，2006.

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

［5］ 胡繪斌，柴舜連，毛鈞杰.基于三維PE 方法的雷達(dá)波傳播損耗估計(jì)［J］.微波學(xué)報(bào)，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］ 李永順.基于拋物線方程求解電大尺寸目標(biāo)雙(多)站RCS 方法的研究［D］.合肥:安徽大學(xué)，2005.

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

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

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

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

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

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

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