瞬變電磁場模擬中的CPM L吸收邊界條件

姜彥南, 劉 文, 王 嬌,3, 張文翠（1.廣西無線寬帶通信與信號處理重點實驗室,桂林 541004；.桂林電子科技大學信息與通信學院,桂林 541004；3.桂林電子科技大學廣西信息科學實驗中心,桂林 541004）

瞬變電磁場模擬中的CPM L吸收邊界條件

姜彥南1,2,3, 劉 文2, 王 嬌2,3, 張文翠2
（1.廣西無線寬帶通信與信號處理重點實驗室,桂林 541004；2.桂林電子科技大學信息與通信學院,桂林 541004；3.桂林電子科技大學廣西信息科學實驗中心,桂林 541004）

提出用于瞬變電磁法（TEM）模擬的時域有限差分（FDTD）算法中的卷積完全匹配層（CPML）吸收邊界條件.首先,在磁場散度方程顯式處理條件下,計算磁場z分量時,構造出一種計算時域卷積項的方法,并詳細推導計算磁場z分量的表達式.而后,對均勻半空間模型進行數值計算.結果表明,提出的CPML吸收邊界條件在保證計算精度的前提下,實現了瞬變電磁法時域有限差分高效正演計算.

瞬變電磁法；時域有限差分；完全匹配層；CPML吸收邊界條件

0 引言

瞬變電磁法（Transient Electromagnetic Methods,TEM）是一種建立在人工可控源基礎上的時間域電磁勘探方法[1－2],因其具有對低阻體反應靈敏、體積效應小、縱橫向分辨率高、施工方便、效率高等優點,被廣泛應用于金屬礦勘探、煤礦水文地質調查和工程勘查等領域,是目前地質勘探中一種重要的方法.

隨著計算機技術的發展,基于有限元、有限差分等數值計算的瞬變電磁法正演研究也得到了快速推進.在有限差分正演計算方面,國外具有代表性的研究主要有：Oristaglio和Hohmann[3]開創性地采用DuFort-Frankel有限差分法研究了二維地電斷面下線源產生的瞬變電磁場；Wang和Hohmann[4－5]以地面均勻半空間瞬變電磁場的解析解為初始條件,采用改進的DuFort-Frankel有限差分形式研究了三維時間域正演算法；Commer和Newman[6－7]基于并行有限差分法研究了電偶源瞬變電磁三維正演算法.國內相關的研究開展稍晚,比如閆述[8]、劉云[9]等人利用DuFort-Frankel有限差分法分別采用總場和二次場求解法計算了二維瞬變電磁場,并對均勻半空間中異常體的瞬變電磁場分布特征進行了分析；孫懷鳳等[10]將瞬變電磁激發源直接加入安培環路定理,并考慮一次場計算和關斷時間的影響,采用時域有限差分（Finite Difference Time Domain,FDTD）法對三維模型進行了計算和分析.

由于計算機的局限性,基于有限差分的瞬變電磁數值模擬只能在有限區域內進行,這樣必須對計算區域的邊界進行截斷處理.分析現有文獻[3－10]發現,前人模擬瞬變電磁場過程中通常采用以下方式：假設計算區域足夠大,使得地下三個截斷邊界離激勵源和目標的距離足夠遠,此時地下目標擴散場在邊界處滿足Dirichlet條件,從而直接將邊界值設為零.該方式的優點是邊界條件設置簡單、方便、可靠性高,且滿足實際的物理規律；缺點是由于計算域范圍足夠大,數值正演計算時占用大量的計算內存,耗費大量的計算時間,導致基于有限差分的TEM計算效率極低.

為了克服以上缺點,在FDTD方法模擬瞬變電磁場時,我們采用縮小計算域范圍,并在地下截斷邊界處設置合適吸收邊界的形式加以實現.目前應用較多的吸收邊界條件有Mur吸收邊界條件、廖氏吸收邊界條件和完全匹配層（Perfectly Matched Layer,PML）吸收邊界條件.Mur吸收邊界條件[11]雖然簡單方便、占用內存小,但計算精度不高,在角區域存在較大誤差；二階近似的Mur吸收邊界條件雖然精度有所改善,但對低頻感應場的吸收效果差.廖氏吸收邊界條件[12]能夠很好地吸收來自計算區域的電磁波,但它的高階吸收邊界條件不穩定.PML邊界條件[13－15]給吸收邊界條件研究帶來了本質性的突破,但它采用分裂電磁場的方式來描述各向異性材料,對計算機的硬件要求很高,并且在模擬瞬變電磁波時對地下高有耗介質中的波和倏逝波的吸收效果不佳.在PML基礎上發展的UPML方法（uniaxial PML）[16,19]直接使用各向異性材料而不改變電磁場形式,能作為瞬變電磁波的吸收邊界條件,但它的遞推公式與材料特性有關,并且對大角度入射波及瞬變電磁波里低頻成份的吸收效果差.目前最具潛力的當屬卷積PML（Convolutional PML,CPML）吸收邊界條件[20],它除了具有UPML吸收邊界條件的性能外,還能對倏逝波、晚時反射波以及低頻感應場進行有效的吸收,因而CPML吸收邊界條件特別適合用做瞬變電磁法FDTD正演數值計算中截斷邊界的處理.

本文中,我們對CPML吸收邊界條件進行研究,以用于瞬變電磁場時域有限差分算法計算區域的截斷處理.

1 FDTD模擬TEM基本原理

準靜態條件下,均勻、有耗、非磁性媒質中,無源Maxwell方程組為[4]

其中E為電場強度,B為磁感應強度,σ為電導率,μ0為真空中的磁導率.由于后期瞬變電磁場變化較緩慢, （1）式中忽略了位移電流項,缺少對電場的時間導數,因此方程組（1）～（4）式不能構成顯式的時域有限差分方程,通常采用引入虛擬位移電流項的方式進行處理,則（1）式修改為

式中,γ具有介電系數的量綱,稱為虛擬介電系數,其取值必須滿足穩定性條件并且使引入的虛擬位移電流項不影響計算結果.

研究發現,在進行瞬變電磁場數值計算的過程中,公式（3）必須顯式包含在迭代過程中,否則對低頻電磁場的計算結果是錯誤的[4].因此,在下面章節中,我們將結合式（2）、（3）和（5）進行計算.對于磁場的計算,可采用先通過式（2）求解磁場x、y兩個分量,然后再通過 Δ·B＝0求解磁場z分量的方式獲得,即

在二維TMy模式下,只有Ey、Bx和Bz三個場分量,因此標量形式的Maxwell方程組為

為了便于分析,采用均勻網格將以上三式進行FDTD離散.這里僅給出與式（9）對應的離散式

值得注意的是,（10）式中Bz的求解采用從模型的底邊界向上步進到地－空邊界的遞推方式.

2 FDTD模擬TEM響應中的CPML吸收邊界條件

由式（2）、（3）和（5）可知,頻域中CPML區域內的約束方程組為

其中,算子 Δs定義為[21]

sx,sy,sz稱為坐標拉伸因子,具有以下形式

其中,ε0為自由空間的介電常數,ω為角頻率,σw是CPML區域內沿w方向的電導率,κw為網格延拓因子, αw是為改善對消逝波和低頻波的吸收而引入的參數,它們分別為[21]

式中

其中εr為背景媒質的相對介電常數,δ為網格尺寸,ρ為各場分量到計算域與CPML交界面之間的距離,d 為CPML厚度,m和mα為多項式分布的階數,σfactor為誤差控制系數,σw,max、κw,max、αw,max分別為相應參數在w方向上的最大值.

定義傅里葉逆變換關系

由傅里葉變換理論及式（15）、（19）得[21]

其中f－1表示傅里葉逆變換算子,u（t）和δ（t）分別是階躍函數和脈沖函數,t為時間.因此,在二維TMy模式下,應用于CPML區域的（11）式和（12）式轉換成時域形式為[20]

式（21）、（22）的FDTD離散表達式可采用文獻[22]的方法求得,而（13）式的CPML吸收邊界條件公式在已掌握的文獻中未見報道.下面重點推導對應于（13）式的CPML吸收邊界條件.

將（14）式代入（13）式并展開得

根據（20）式時－頻域轉換關系可得（23）式的時域形式為

考察（24）式不難發現,該式不僅左右兩邊都包含卷積項,而且卷積項形式為ζw（t）?Bw/?w,而CPML層中常規的方程表達式只有一邊包含卷積項,且形式為ζv（t）?Bw/?v（如（21）,（22）式）,因此常規的FDTD離散方法無法求解（24）式.為此,針對（24）式中卷積項的計算構造如下：

將（25）式中x和z分量的時域卷積項分別代入（24）式兩邊,并采用中心差分離散得

整理后,即得

其中,

從（27）式可以看出,在處理CPML吸收邊界時,磁場Bz的推導也是從底邊界向上步進到地－空邊界,這與計算域內磁場Bz計算過程完全一致.對比式（10）可以發現,（27）式中相對于常規FDTD迭代計算式多出了（1＋κxcx）κz/（1＋κzcz）κx和等號右邊后兩項.在非CPML區域的FDTD計算中所有對應的κw＝1以及σw＝0,從而cw＝0,（1＋κxcx）κz/（1＋κzcz）κx＝1,后兩項為0.因此,可以把非CPML區域和CPML區域的場分量迭代公式合并為一式.

由以上分析可以看出,式（21）,（22）和（24）就構成了二維TMy模式下瞬變電磁法FDTD計算的CPML吸收邊界條件.

3 算例驗證及分析

圖1 CPML吸收邊界條件截斷的均勻半空間模型Fig.1 Sketch of homogeneous half-space model truncated by CPML ABC

為了驗證分析本文提出的CPML吸收邊界條件的有效性及可靠性,考慮如圖1所示的均勻、有耗、非磁性二維半空間模型.坐標原點位于地表面上,水平方向為x軸,垂直向下為z軸正方向.空間步長Δx＝Δz＝10 m,FDTD計算區域為60×30個網格.地下介質電導率σg＝1/300（s·m－1）.在x＝－100 m和x＝100 m的地面上（z＝0）分別通以I1＝1 A和I2＝－1 A的無限長線電流源s1和s2.t＝0時刻斷開電流,初始時刻t0＝1.047 2μs；時間步長Δtn＝a,虛擬介電常數γ ＝D/μ0（Δtn/δ）2,其中a為控制系數,D為模擬空間的維度,δ為空間步長,n為時間步數,取值分別為a＝0.2,D＝2,δ＝10m,n＝0,1,2,…, Δtn和γ隨時間延遲而逐漸增大；源作為初始條件加入[3],在t＝t0時刻賦給電場初始值,t＝t0＋Δt0/2賦給磁場初始值.地－空邊界采用地面上的垂直磁場Bz向空氣層延拓半個網格的切向磁場Bx得到（如圖1中Γ0所示）,地下三個邊界則采用CPML吸收邊界條件處理.

為了獲得較好的CPML吸收性能,基于大量數值試驗并考慮計算量問題,最終選擇（16）－（18）式中的參數如下：κmax＝6,αmax＝0.05,σfactor＝300,m ＝1,mα＝1,CPML的厚度d為12個網格.設置觀察點A、B、C距離CPML區域內邊界僅一個網格,分別位于網格點（－29,1）、（－29,15）和（－29,29）處,用于定量分析CPML的處理效果.本節均基于此例進行分析.

3.1 計算結果正確性分析

線電流源在地下均勻半空間的輻射場為向下、向外擴散的“煙圈”.圖2給出了兩個不同時刻的電場Ey等值線圖.從圖中可以看出,在截斷邊界處,“煙圈”呈規則的形狀,表明本文提出的CPML吸收邊界具有良好的外向行波吸收性能.

圖2 不同時刻電場Ey等值線圖Fig.2 Contours of electric field Eyat different times

圖3 坐標原點電動勢衰減曲線Fig.3 Decay curves of EMF at origin of coordinate

圖3給出了坐標原點處垂直感應電動勢（Electromotive force,EMF.垂直感應電動勢定義為垂直磁感應強度對時間的一階導數,即 EMF ＝?Bz/?t＝－?Ey/?x）衰減曲線,其包括以下三種情況的計算結果：①解析解；②計算域足夠大,地下截斷邊界采用Dirichlet條件的數值解,即傳統的FDTD數值解；③計算域縮小至60×30個網格,地下截斷邊界采用CPML處理的FDTD數值解.由圖可看出,三種情況下的計算結果具有很好的一致性,表明本文提出的算法和編寫的計算程序是正確的.

3.2 計算結果誤差分析

為定量分析本文提出的CPML吸收邊界性能,定義相對反射誤差值（dB）為[20]

其中,ξ代表觀察點處的測試值,ξref代表觀察點處ξ的參考值,ξref,max代表觀察點處ξ的參考值在整個數值模擬時間內的最大值.

1）不同位置處的電場相對反射誤差分析

圖4是在CPML吸收邊界條件下,程序運行1 000個時間步迭代至1.1 ms,在A、B、C觀察點處計算得到的電場相對反射誤差曲線.其中,參考值采用傳統的FDTD方法獲得.由于（30）式中參考值最大值隨著深度增大而減小,因此隨著深度的增加相對反射誤差逐漸增大,但最大誤差僅為－35 dB.結果表明本文提出的CPML吸收邊界處理后的數值計算實現了較高的計算精度.

2）CPML層數不同時的電場相對反射誤差分析

圖5給出了CPML層分別為10、15和20個元胞時,觀察點B處的電場相對反射誤差曲線.其中,參考值同樣采用傳統的FDTD方法獲得.從圖中可以看出,隨著CPML層數的增加,吸收效果逐漸改善.考慮內存資源和計算時間等因素,CPML層數可在10～15之間選擇.

3）不同數值解下的地表感應電動勢相對反射誤差分析

為了進一步分析本文提出的CPML吸收邊界的效果,我們采用“偽FDTD方法”計算垂直感應電動勢.（令圖1所示的模型不含CPML吸收層,直接將邊界處的場值設為零,而后采用傳統FDTD方法進行計算.但此時不滿足計算區域足夠大的條件,故為了描述需要稱其為“偽FDTD方法”）.

圖4 觀察點A、B、C處的相對誤差Fig.4 Relative errors at points A,B and C

圖5 B處不同CPML層數的相對誤差Fig.5 Relative errors of CPML with different layer number at point B

圖6 不同數值解下的地表感應電動勢相對誤差Fig.6 Relative errors of EMF at the surface under different numerical solution

地表x＝290 m處垂直感應電動勢的三種數值解相對反射誤差如圖6所示（其中,參考值采用解析解）.可以看出,偽FDTD方法數值解的相對反射誤差大于0dB,在截斷邊界處產生了很大的反射；而采用CPML吸收邊界條件處理后,相對反射誤差大幅度減少,在截斷邊界處幾乎不產生反射,達到了與傳統FDTD數值計算接近的計算精度.結果表明本文提出的CPML吸收邊界具有很好的吸收效果.

3.3 計算資源占用及計算效率分析

對于二維TMy模式,電磁場只有三個分量,設各場分量采用雙精度型變量（每個變量占內存8個Bytes）,FDTD區域內元胞總數為N.如果在目標建模時不直接給出每個元胞電磁參數,而是給出每種介質的介質編號,則識別各個場分量離散位置處的介質參數需要1個字節的整數.如果每個網格各場分量有相同的介質參數,則有[22]：

根據（31）式,在圖1模型中未引入CPML吸收邊界條件且計算域足夠大的情況下,設其網格數為800× 400,則計算三個場量需要內存7.63 MB,運行1 000步耗時55.6秒.而引入CPML吸收邊界條件后,總網格數降為84×42,即使在完全匹配層內因處理CPML吸收邊界條件而引入了ψ變量,計算三個場量所需內存也僅為0.12 MB,運行1 000步耗時0.79秒.可以看出,本文提出的CPML吸收邊界條件可以有效地節省計算內存,大幅縮短計算時間,明顯提高數值計算效率,這使得基于FDTD正演計算的TEM快速反演成為可能.

4 結論

在瞬變電磁法的傳統FDTD正演計算中,因截斷邊界要求足夠遠,導致模擬區域范圍廣、數值計算量大、迭代時間長、計算效率低等問題.本文提出的應用于瞬變電磁法FDTD正演計算的CPML吸收邊界條件不僅保證了FDTD數值計算的精度,而且克服了傳統FDTD正演計算的缺點.

本文方法可進一步推廣到三維情形,以用于地下三維目標TEM響應的FDTD正演計算.

[1] Niu Z L.The theory of time-domain electromagnetic methods[M].Changsha：Central South University of Technology Press, 1992：1.

[2] Nabighian M N.Quasi-static transient response of a conducting half-space—An approximate representation[J].Geophysics, 1979,44（10）：1700－1705.

[3] Oristaglio M L,Hohmann GW.Diffusion of electromagnetic fields in two-dimensional earth：a finite difference approach[J].Geophysics,1984,47（7）：870－894.

[4] Wang T,Hohmann W G.A finite-difference,time-domain solution for three-dimensional electromagnetic modeling[J].Geophysics,1993,58（6）：797－809.

[5] Wang T,Tripp A C,Hohmann GW.Studying the TEM response of a 3-D conductor at a geological contact using the FDTD method[J].Geophysics,1995,60（4）：1265－1269.

[6] Newman G A,Commer M.New advances in three-dimensional transient electromagnetic inversion[J].Geophys J Int,2005, 160（1）：5－32.

[7] Commer M,Newman G.A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources [J].Geophysics,2004,69（5）：1192－1202.

[8] Yan S,Chen M S,Fu JM.Direct time-domain numericalanalysis of transientelectromagnetic fields[J].Chinese JGeophys, 2002,45（2）：275－284.

[9] Liu Y,Wang X B,Wang Y.Numericalmodeling of the 2D time-domain transient electromagnetic secondary field of the line source of the current excitation[J].Applied Geophysics,2013,10（2）：134－144.

[10] Sun H F,LiX,LiSC,etal.Three-dimensional FDTDmodeling of TEM excited by a loop source considering ramp time[J].Chinese JGeophys,2013,56（3）：1049－1064.

[11] Mur G.Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations[J].IEEE Trans Electromagn Compat,1981,23（4）：377－382.

[12] Liao Z P,Huang K L,Yang B P,etal.A transmitting boundary for transientwave analyses[J].Scientia Sinica（Series A）, 1984,27（10）：1062－1076.

[13] Berenger JP.A perfectlymatched layermedium for the absorption of electromagnetic waves[J].JComput Phys,1994,114：185－200.

[14] Berenger JP.Perfectlymatched layer for the FDTD solution ofwave-structure interaction problem[J].IEEE Trans Antennas Propagat,1996,44（1）：110－117.

[15] Zheng Z,Yang L.3D modified novel PML absorbing boundary condition for truncating magnetized plasma[J].Chinese J Comput Phys,2013,30（6）：895－901.

[16] Sacks Z S,Kingsland D M,Lee D M.A perfectly matched anisotropic absorber for use as an absorbing boundary condition [J].IEEE Trans Antennas Propagat,1995,43（12）：1406－1463.

[17] Gedney SD.An anisotropic perfectly matched layer absorbing media for the truncation of FDTD lattices[J].IEEE Trans Antennas Propagat,1996,44（12）：1630－1639.

[18] Chai Y,Sun J,Li L,etal.Unifiedmodeling and optimization of UPML absorbing boundary conditions[J].Chinese JComput Phys,2011,28（3）：413－419.

[19] Yang T C,Li H,Zhang H.Forward modeling of GPML absorbing boundary conditions in GPR simulation[J].Chinese J Comput Phys,2014,31（2）：200－206.

[20] Roden JA,Gedney SD.Convolution PML（CPML）：an efficient FDTD implementation of CFS-PML for arbitrarymedia[J].Microwave Opt Tech Lett,2000,27（12）：334－339.

[21] Taflove A,Hagness SC.Computational electro-dynamics：The finite-difference time-domain method[M].Norwood,Artech House,2005：289－310.

[22] Ge D B,Yan Y B.Finite difference time domain method of electromagnetic wave（3rd edition）[M].Xi’an：Xidian University Press,2011：206－207.

CPM L Absorbing Boundary Condition in M odeling of Transient Electromagnetic Fields

JIANG Yannan1,2,3, LIUWen2, WANG Jiao2,3, ZHANGWencui2
（1.Guangxi Key Laboratory ofWirelessWideband Communication＆Signal Processing,Guilin 541004,China；
2.School of Information and Communication Engineering,Guilin University ofElectronic Technology,Guilin 541004,China；
3.Guangxi Experiment Center of Information Science,Guilin University of Electronic Technology,Guilin 541004,China）

A scheme of convolutional perfectly matched layer（CPML）absorbing boundary condition（ABC）is proposed,which is used to truncate finite-difference time-domain（FDTD）method modeling transient electromagnetic（TEM）response.It derives specially CPML formulation dealing explicitly with divergence of magnetic induction and conceives an algorithm calculating convolutional term of z component.Then,the scheme is validated by numerical modeling of a homogeneous half-space model.The results indicate that efficient calculation of TEM is achieved.

transient electromagneticmethod；finite-difference time-domain；perfectlymatched layer；convolutional PML absorbing boundary condition

1001-246X（2015）06-0701-08

P631

A

2014－11－21；

2015－02－09

國家自然科學基金（61001020）、廣西自然科學基金（2014GXNSFAA118283）、廣西研究生教育創新計劃（GDYCS201415）和廣西信息科學實驗中心資助項目

姜彥南（1980－）,男,博士,副教授,主要從事電磁場數值計算、電磁輻射與散射研究,E-mail：ynjiang＠guet.edu.cn