頻率域電磁波優化25點有限差分正演模擬

中國電子科技集團公司第二十研究所 李宏亮 黨 杰

時間域電磁波有限差分方法(FDTD)是一種常用的電磁波數值模擬方法，通過將電磁波進行網格化處理，使得電場和磁場均勻的分布在網格上，從而求解電磁波方程。本文基于麥克斯韋方程組，開展了基于頻率域電磁波有限差分的相關研究。為了驗證頻率域有限差分的可行性以及提高數值模擬精度，主要采用了兩種方法：第一種方法是在麥克斯韋方程組的基礎上直接推導出頻率域電磁波方程，避免了時間域遞推造成的累計誤差；第二種方法是采用優化25點差分方法構建頻率域電磁波方程算子，不僅保證了模擬的精度，而且能有效的抑制數值頻散。此外，為了驗證算法的可行性，在PML邊界條件下進行了頻率域電磁波的數值模擬。

隨著電磁波理論研究與實際應用的推廣，麥克斯韋方程組的數值解法也被廣泛的關注和研究。已經提出的方法有頻率域的矩量法、有限元法、有限差分法以及時間域的傳輸線矩陣法、有限差分等方法。各種方法都有自己的優勢以及應用上的局限性，要么計算過于復雜，要么占用內存大、計算速度慢，在實際應用中往往需要將幾種方法混合使用。

有限差分方法的中心思想是用差分代替微分，通過將研究區域網格化處理，使電場和磁場有規律的分布在網格點處，從而得到網格區域的電場和磁場值。本文通過推導麥克斯韋方程組得到頻率域二維電磁波方程，然后整理成矩陣方程形式，采用優化25點有限差分方法構建頻率域差分算子，對算子進行稀疏化存儲，使得頻率域數值模擬更加簡便快捷，同時也證明了該方法的有效性。

1 頻率域電磁波優化25點有限差分

1.1 頻率域電磁波方程推導

麥克斯韋方程組的微分形式是這樣的，如式(1)：

式(1)中，E表示電場強度，B表示磁場強度；ρ表示電荷密度，J表示電流密度；ε0和μ0分別表示真空中的介電常數和磁導率；是矢量微分算子，和分別表示求取散度和旋度。

為簡化后續推導，考慮真空、均勻介質(即ρ=0，J=0)情況下，麥克斯韋方程組變為式(2)：

式(2)方程兩側同時取旋度，并作進一步簡化，得到式(3)：

考慮二維情況，可以將式(3)改寫為式(4)：

式(4)兩側對時間做傅氏變換，即：

式(5)中，ω表示角頻率；和分別是電場強度和磁場強度的頻率域表示。

1.2 基于優化25點有限差分方法的算子離散化

圖1所示為優化25點有限差分方法中各個權系數分布示意圖。在求解時，需要用到該點周圍的25個點，而這25個點對于要求解的二階差分算子“貢獻程度”是不同的。如圖1(a)，每一行需要用到5個點來構造2個差分算子，根據每一行中各列元素距離待計算點的距離遠近，設置加權系數c、d，將他們結合起來作為這一行的差分算子，因此可以得到5個差分算子。然后，再將這5個差分算子通過加權系數b1、b2、b3進行加權平均即可得到最終的離散形式。同理，也通過類似的方法可以獲得，如圖1(b)所示。圖1(c)是方程左端加速項的差分格式，其權系數以差分點為圓心，分布在一組同心圓上。

圖1 優化25點差分示意圖

為了使網格色散和數值各向異性盡可能最小，必須確定一組最優的加權系數使相速度和群速度盡可能的接近。求取這些系數可以采用梯度法或牛頓法，通過給定一組初始值，不斷地迭代搜索以使該系數能夠同時使相速度和群速度接近。

本文采用Min利用高斯牛頓法計算出來的最優加權系數：

1.3 吸收邊界條件及離散方程

電磁波在空間中傳播時相當于是在一個半無窮的介質中傳播，沒有邊界的影響。但在數值模擬過程中，由于計算區域有限，而且人為的引入了反射界面，這就會導致出現一些在實際中不會出現的反射波，影響最終的模擬結果。因此，必須構造一個吸收邊界，以期望盡可能的抑制這種邊界反射，從而更加真實的模擬電磁波的傳播過程，研究電磁波的傳播特性。

本文采用PML（如圖2所示）吸收邊界來抑制邊界反射。當地震波傳播到邊界時，隨著進入到吸收層內傳播距離的增加，電磁波不斷地衰減，這樣當波傳播到最外層邊界時，其能量被完全衰減就不會產生邊界反射。當電磁波傳播到區域1時，沿x和z方向都衰減；當波傳播到區域2時，沿x方向不衰減，沿z方向衰減；當波傳播到區域3時，沿x方向衰減，沿z方向不衰減。

圖2 PML吸收邊界示意圖(引自王守東,2003)

具體的，將這種邊界施加到電磁波方程時，需要做一個坐標拉，這樣施加PML吸收邊界的電磁波方程就可以表示為式(9)：伸的變換

用d來表示衰減因子，使得PML吸收區域內d≠0，而有效計算區域內d=0。定義，得：

則：

其中：

式中A是在0.5～1.79之間取值的常數，通常需要通過不斷地試驗確定最佳的參數，后文中數值模擬時采用A=1.3；f0是主頻；LPML是PML層的厚度；lf是PML邊界與內部有效計算區域的距離。

將衰減因子代入電磁波方程(5)式，得到：

根據文獻，為簡化后續推導，通常省略式(12)中一階偏導項，亦即：

將差分加權系數以及式(10)～式(13)帶入到式(9)中，施加PML邊界條件，合并同類項后得到最終的離散方程為：

式中參數為：

方便起見，將式(14)和式(15)聯合起來整理成矩陣方程形式：

式(16)中，A是由差分系數構成的一個具有25條對角線的系數矩陣，該矩陣是一個大型的、非對稱的稀疏矩陣；U是由和構成的一個列向量；F是施加的激發源。若網格區域為N=m×n，則A是一個2N×2N的矩陣，U是2N×1的向量，F是2N×1的向量。

2 數值模擬

2.1 數值模擬

為了驗證各種激發源加載的有效性，在均勻各向同性介質中進行數值模擬實驗。模型的網格大小為160×160，x和z方向的網格間距為dx=dz=0.005m。圖3是激發源主頻為3GHz，在模型中心位置處激發得到的單頻切片。

圖3 集中力源單頻切片

2.2 精度對比

為了測試差分點數對數值模擬精度的影響，在相同參數下，采用網格大小201×201，網格間距Δx= Δz=0.005m，激發源主頻為3GHz，放置在模型中心處。根據式(16)分別采用9點差分方法和25點差分方法構建系數矩陣，進行電磁波正演模擬。圖4是將頻率域各個頻點單頻切片匯總后進行反傅里葉變換得到的時間域某時刻波場快照。可以看出，通過優化25點差分方法構建正演算子得到的模擬結果較之常規方法而言能夠很好地抑制數值頻散。

圖4 波場快照

結論與展望：本文在麥克斯韋方程組的基礎上，推導出了真空中電磁波的頻率域表達式，并采用優化25點有限差分方法對電磁方程進行了差分，得到了基于優化25點有限差分方法的離散電磁場方程。通過添加PML邊界條件，進行了頻率域電磁波數值模擬，得到了以下認識和結論：

（1）采用25點有限差分方法進行電磁波數值模擬具有更高的精度；

（2）由于頻率域系數矩陣只與差分系數有關，故系數矩陣A只需要計算一次，這使得頻率域電磁波數值模擬計算速度快，很適合并行計算，為理論研究和工程實踐提供了一個較好的實現思路，為電磁波成像和反演奠定了一定的基礎。

當然，這種方法從目前來看仍有一定的限制，比如從二維模型轉向實際三維模型時系數矩陣需要的巨大內存或者說是系數矩陣稀疏化如何處理，以及實際工程應用中大量的頻率域數據如何轉化等問題都有待進一步考慮。