多源并發(fā)下Mur二階吸收邊界和非分裂遞歸卷積完全匹配層對比研究

崔凡，陳毅，薛晗鵬，彭蘇萍，杜云飛

(1.中國礦業(yè)大學(xué)(北京) 地球科學(xué)與測繪工程學(xué)院，北京 100083;2.中國礦業(yè)大學(xué)(北京) 煤炭資源與安全開采國家重點(diǎn)實(shí)驗(yàn)室，北京 100083)

0 引言

目前，在時(shí)域有限差分方法中運(yùn)用邊界條件來截?cái)嗑W(wǎng)格存在很多方法。Bayliss和Turkel[1]采用外行波的模擬法作為吸收邊界條件，實(shí)現(xiàn)了對電磁波的簡單吸收；Engquist和Majda[2]提出了基于單向波動方程的Engquist-Majda吸收邊界條件；Mur[3]給出了波動方程的各階近似及差分形式，表明了在高階近似時(shí)，對電磁波的吸收效果較好，然而對于各向異性介質(zhì)和色散介質(zhì)，以二階近似Mur作為吸收邊界條件會在網(wǎng)格截?cái)嗵幇l(fā)生強(qiáng)反射。之后，Berenger[4]提出了完全匹配層吸收邊界條件，他將電磁波在吸收邊界區(qū)域進(jìn)行波場分裂，并對各個(gè)分裂場賦予不同的損耗值，這相當(dāng)于在FDTD網(wǎng)格外加入一層吸收媒介，它具有不依賴于外向傳播的電磁波入射角及頻率的波阻抗，因此，進(jìn)入完全匹配層內(nèi)的電磁波快速地衰減，實(shí)現(xiàn)了對不同頻率、不同角度的入射波較好的吸收效果。有很多研究發(fā)現(xiàn)[5-7]，Berenger完全匹配層作為吸收邊界條件比旁軸近似吸收邊界條件、Higdon吸收邊界條件[8]、廖氏吸收邊界條件[9]和指數(shù)衰減吸收邊界條件[10]具有更好地吸收效果。但是，Berenger的PML理論違背了波的折射原理，在Maxwell方程中是不能實(shí)現(xiàn)的，同時(shí)，電磁場的分裂增加了數(shù)值計(jì)算的難度，計(jì)算效率低，而且只對行波具有吸收效果，對空域中衰減的隱失波、低頻波以及入射角度較小的掠射波吸收效果差[11]。為了改善邊界條件的吸收效果，Sacks等[12]，Gedney等[13]提出了單軸各向異性完全匹配層作為吸收邊界條件，該方法在數(shù)學(xué)上與Berenger提出的經(jīng)典完全匹配層等價(jià)，它是基于Maxwell方程的，不分裂波場，而是在吸收系數(shù)中引入線性吸收因子，實(shí)現(xiàn)了對隱失波、低頻波以及在PML層內(nèi)凋落波等干擾波的吸收，該方法作為良好的吸收邊界條件在差分算法中被廣泛使用。如肖明順等[14]、詹應(yīng)林等[15]將各向異性完全匹配層應(yīng)用于二維探地雷達(dá)的正演中；中南大學(xué)馮德山等[16-18]實(shí)現(xiàn)了在二維和三維空間中將各向異性完全匹配層引入交替方向隱式時(shí)域有限差分算法中(ADI-FDTD)；吉林大學(xué)的李靜等[19]人開展了三維各向異性完全匹配層作為吸收邊界條件的高階時(shí)域有限差分正演。同一時(shí)間段內(nèi)，Chew等[20]提出了拉伸坐標(biāo)完全匹配層理論；Kuzuoglu等[21]在復(fù)頻域內(nèi)的完全匹配層變量Sk中引進(jìn)低頻分量吸收參數(shù)，將復(fù)平面的極點(diǎn)從實(shí)軸移動到虛軸，加強(qiáng)了對低頻波和隱失波的吸收。然而到現(xiàn)研究階段，所有的吸收邊界條件都是在單個(gè)激勵(lì)源的條件下對電磁波進(jìn)行數(shù)值模擬。

本文基于Maxwell方程推導(dǎo)了一種在時(shí)域有限差分算法中運(yùn)用遞推卷積求解電磁場的方法，該方法在離散網(wǎng)格條件下通過遞推形式計(jì)算卷積，不分裂變量，直接計(jì)算卷積，避免了對卷積求解的復(fù)雜計(jì)算。并將該方法運(yùn)用到多源并發(fā)下對電磁波的吸收，為陣列探地雷達(dá)數(shù)值仿真做了鋪墊。

1 基本原理

1.1 坐標(biāo)伸縮變化下Maxwell卷積方程

非遞歸卷積完全匹配層主要求解是在麥克斯韋方程的頻率域進(jìn)行，對連續(xù)和時(shí)諧場，麥克斯韋旋度支配方程可寫為：

(2)

在PML中電導(dǎo)率不是由物理空間模型的基本參數(shù)決定的，而是設(shè)置電導(dǎo)率參數(shù)使網(wǎng)格截?cái)嗟姆瓷渥钚?。在此意義上，電導(dǎo)率σ具有任意性，將一個(gè)從屬于特定位置的相對介電常數(shù)來規(guī)范電導(dǎo)率，在分裂場PML條件下定義：

(3)

當(dāng)頻率趨于0時(shí)，Sw沒有意義。這將導(dǎo)致低頻引發(fā)數(shù)值異常，為修正這個(gè)問題，增加額外的參數(shù)來保證頻率趨于0時(shí)Sw值的有限性。引入更一般化的坐標(biāo)伸縮因子Si：

(4)

式中：Ki是改善PML對表面波吸收的參數(shù)；σi為PML層內(nèi)k方向電導(dǎo)率參數(shù)；αi為改善PML對低頻分量的吸收參數(shù)；σi和αi都為大于零的正實(shí)數(shù)；Ki≥1。在三維空間下，i∈{x,y,z}，每一個(gè)Si項(xiàng)總是與i∈{x,y,z}方向的道數(shù)是成對的。有耗介質(zhì)中，在TM極化模式下(Hx,Hy,Ez)，三維空間下電場頻域方程表達(dá)式為：

(5)

將方程(5)通過傅里葉逆變換轉(zhuǎn)化到時(shí)域，由于坐標(biāo)伸縮因子與頻率無關(guān)，所以在時(shí)域內(nèi)方程右邊存在卷積，即：

(7)

對坐標(biāo)伸縮因子的倒數(shù)求傅里葉逆變換有：

其中δ(t)是Dirac沖擊函數(shù)，u(t)是單位階躍函數(shù)。定義：

i=(x,y,z) 。

(9)

因?yàn)镈irac沖擊函數(shù)和其他函數(shù)的卷積仍為原函數(shù)，即δ(t)*f(t)=f(t)，將式(10)代入式(6)有：

(11)

式(11)為電場z分量的時(shí)域遞推卷積方程，直接在遞推方程(14)中計(jì)算時(shí)域卷積效率是很低的，為有效計(jì)算卷積，下面給出遞歸卷積求解原理和在時(shí)域有限差分方法中的實(shí)現(xiàn)步驟。

1.2 基于遞歸卷積的非分裂場PML的FDTD實(shí)現(xiàn)

(12)

式中：Eui表明這個(gè)函數(shù)將會在電場Eu分量上更新，并與i(i=x,y,z)方向的空間導(dǎo)數(shù)有關(guān)。假設(shè)積分變量τ在式(12)中是連續(xù)變化的，由于時(shí)域有限差分方法的離散性，所以Hv場是離散變化的。用連續(xù)變量t表示Hv時(shí)，只取離散值，即

(13)

卷積包含f(qΔt-τ)函數(shù)。在時(shí)間步為q=0時(shí)，函數(shù)值為f(-τ)，如圖1b所示。在所有采樣點(diǎn)fn都關(guān)于原點(diǎn)反轉(zhuǎn)對稱，在-Δt≤τ<0上函數(shù)值為f1，在-2Δt≤τ<-Δt上函數(shù)值為f2，在-3Δt≤τ<-2Δt上函數(shù)值為f3，其余各時(shí)刻的函數(shù)值以此類推。

圖1 函數(shù)f(t)階梯表示法Fig.1 Stepped representation of function f(t)

圖1c展示了在特定時(shí)刻q=4時(shí)f(qΔt-τ)的函數(shù)值。如圖所示，擴(kuò)展到τ=0右側(cè)右邊的第一個(gè)脈沖具有fq的值，擴(kuò)展到τ=Δt右邊的值為fq-1，擴(kuò)展到τ=2Δt右邊的值為fq-2,以此類推。此時(shí)移動函數(shù)可表示為：

(14)

對磁場求導(dǎo)數(shù)有：

(15)

(16)

式(16)中，脈沖函數(shù)pk(t)用來確定積分上下限?？紤]如下特殊積分：

(17)

把式(16)中的ζi(t)離散沖擊響應(yīng)定義為：

(18)

式中：

i=(x,y,z) ，

(19)

將式(17)、(18)、(19)、(20)代入式(16)中有：

(21)

將k=0項(xiàng)和其他項(xiàng)分離有：

(22)

令n=k-1代替指數(shù)項(xiàng)，即k=n+1：

(23)

用k來表示指數(shù)項(xiàng)n有：

(24)

對比式(24)和式(21)有：

(25)

空間離散下的電場遞推式，按照Yee氏網(wǎng)格將式(25)進(jìn)行時(shí)間和空間上的離散，可得：

(26)

上式中Z0i(k)項(xiàng)含有卷積的計(jì)算，離散卷積計(jì)算十分復(fù)雜，但是同時(shí)Z0i(k)項(xiàng)是簡單的指數(shù)形式，從而它們的和可以通過遞歸卷積來得到，引入一組新的輔助表達(dá)式φi：

(i=x,y)

(27)

將φi代入式(26)經(jīng)整理后得到電場Ez分量的遞推表達(dá)式：

(28)

式(28)中：σ為電導(dǎo)率，ε為介質(zhì)的介電常數(shù)，m=(i,j,k)，在仿真時(shí)令模擬區(qū)域的K值為1，在PML層內(nèi)使K值漸進(jìn)變化，PML層為有限的厚度單元，σ、K和α在PML層中單調(diào)變化：

(29)

(30)

(31)

1.3 Mur二階吸收邊界條件的FDTD實(shí)現(xiàn)

在二維直角坐標(biāo)系下，行波的平面波解為：

f(x,y,t)=Aexp[j(ωt-αxx-αyy)] ，

(32)

假設(shè)在x=0的位置設(shè)置截?cái)噙吔纾瑒t在x≥0的區(qū)域會同時(shí)存在入射波和反射波。因此式(32)可寫為：

(33)

可將上式分解為左行波f-(入射波)和右行波f+(反射波)之和：

(34)

令微分算子L使得L±f±=0，其中：

(36)

對左行波使得L-f-|x=0中并令f=Ez，化簡后得到：

(37)

(38)

在(i+1/2,j)點(diǎn)和t=(n+1/2)Δt時(shí)刻對式(62)進(jìn)行差分離散并進(jìn)行線性插值可得到左截?cái)噙吔缣嶮ur的二階吸收邊界條件遞推式：

(39)

式(39)中，v為電磁波傳播速度，μ為磁導(dǎo)率，Δx、Δy、Δt分別為x、y方向的空間步長和時(shí)間步長。式(39)給出的Mur二階吸收邊界條件遞推公式?jīng)]有在邊角點(diǎn)處對電場值進(jìn)行修正，因此，在4個(gè)邊角位置處發(fā)生明顯的波反射，如圖2所示，給出了激勵(lì)源位于模型區(qū)域中心位置處，Ez波場在600個(gè)時(shí)間步的變化，在第400個(gè)時(shí)間步時(shí)，由于沒有對邊角進(jìn)行修正，此時(shí)4個(gè)角點(diǎn)位置處都發(fā)生了明顯的反射，在500和600時(shí)間步時(shí)刻的研究區(qū)域內(nèi)存在4個(gè)角點(diǎn)反射的波場。

圖2 未修正角點(diǎn)不同時(shí)間步的Ez波場Fig.2 Ez wave field at different time steps of uncorrected corners

1.4 Mur二階吸收邊界條件的角點(diǎn)修正

(40)

創(chuàng)建新坐標(biāo)系ηoξ，新坐標(biāo)系于原始坐標(biāo)系逆時(shí)針旋轉(zhuǎn)45°，在新坐標(biāo)系下Mur一階近似公式可表示為：

(41)

圖3 TM波左下角點(diǎn)修正示意Fig.3 Schematic diagram of TM wave lower left corner correction

由于P點(diǎn)與原點(diǎn)在同一離散網(wǎng)格內(nèi)，忽略波振幅衰減，利用線性插值可得到：

。(42)

根據(jù)式(40)結(jié)合式(42),當(dāng)Δx≠Δy時(shí)，有：

(43)

根據(jù)式(43)，修正角點(diǎn)位置處的電場值取決于該角點(diǎn)相鄰一個(gè)時(shí)間步位置處前一時(shí)刻的電場值和該點(diǎn)的速度值。

2 模型驗(yàn)證

2.1 均勻介質(zhì)

為驗(yàn)證在均勻介質(zhì)下Mur二階近似吸收邊界條件和基于遞歸卷積完全匹配層作為吸收邊界條件在多源并發(fā)下對電磁波的吸收效果，設(shè)置模擬區(qū)域網(wǎng)格大小為200×200，四周PML層厚度為10個(gè)網(wǎng)格。模擬區(qū)域的物性參數(shù)為：εr=3.5,ur=1,σ=2.5 mS/m。兩個(gè)激勵(lì)源采用主頻為900 MHz的布萊克曼—哈里斯脈沖，兩個(gè)激勵(lì)源的位置沿x軸相距80個(gè)空間步長，空間步長為0.006 m，時(shí)間步長為 0.015 ns。不同時(shí)間步的Ez波場如圖4所示。圖4為應(yīng)用遞歸卷積完全匹配層作為吸收邊界條件的Ez波場快照，仿真時(shí)間持續(xù)400個(gè)時(shí)間步長。在第100時(shí)間步時(shí)(圖4a)，兩個(gè)并發(fā)激勵(lì)源電磁波相遇，經(jīng)相互干涉作用后繼續(xù)擴(kuò)散；在時(shí)間步200時(shí)(圖4b)，部分電磁波進(jìn)入PML層被吸收掉； 時(shí)間步300時(shí)(圖4c)兩個(gè)激勵(lì)源產(chǎn)生的擴(kuò)散電磁波到達(dá)彼此對面PML層，此時(shí)部分波到達(dá)邊界進(jìn)入PML層無任何反射發(fā)生；時(shí)間步400時(shí)(圖4d)，波擴(kuò)散離開研究區(qū)域，在整個(gè)仿真持續(xù)時(shí)間段內(nèi)，無任何明顯反射發(fā)生，電磁場在網(wǎng)格截?cái)辔恢锰幈籔ML層有效吸收掉了。

圖5給出了Mur二階吸收邊界條件不同時(shí)間步Ez波場快照，二階Mur吸收邊界條件是運(yùn)用微分二階近似來求解單向行波方程，它是將微分方程進(jìn)行泰勒級數(shù)展開后去掉高階項(xiàng)的近似算法，此時(shí)反射系數(shù)不為零，因此不能實(shí)現(xiàn)對波的良好吸收。在兩個(gè)激勵(lì)源并發(fā)的條件下，電磁波相互干涉，重構(gòu)后的波場形成小方形包絡(luò)，如圖5c箭頭所示位置，隨著時(shí)間推移，小方形包絡(luò)在水平方向上擴(kuò)展開來，以平面波束形式達(dá)到網(wǎng)格邊界位置。這時(shí)，Mur二階吸收邊界對電磁波吸收效果較差，會在邊界處發(fā)生反射。在前300個(gè)時(shí)間步內(nèi)(圖5a,5b,5c)，由單個(gè)激勵(lì)源產(chǎn)生的電磁波在均勻介質(zhì)中擴(kuò)散，以球面的形式到達(dá)網(wǎng)格邊界處， 在第400個(gè)時(shí)間步時(shí)(圖5d)，經(jīng)過干涉后的波到達(dá)邊界，此時(shí)部分波返回研究區(qū)域；在時(shí)間步500時(shí)，波已經(jīng)擴(kuò)散離開研究區(qū)，然而因吸收不完全導(dǎo)致部分反射波返回研究區(qū)域，如圖5e、5f所示，此時(shí)波在邊界位置發(fā)生多次反射，形成諧振，對下一時(shí)間段內(nèi)的有效波場造成二次干擾。

圖4 均勻介質(zhì)遞歸卷積完全匹配層不同時(shí)間步Ez波場快照Fig.4 Snapshots of Ez wave field at different time steps of the homogeneous medium recursive convolution perfectly matched layer

圖5 均勻介質(zhì)Mur二階吸收邊界不同時(shí)間步Ez波場快照Fig.5 Snapshots of Ez wave field at different time steps of Mur second-order absorbing boundary in homogeneous medium

2.2 低洼模型

在復(fù)雜地質(zhì)構(gòu)造條件下，在地下介質(zhì)空間中擴(kuò)散的電磁波會產(chǎn)生反射、繞射現(xiàn)象。繞射波會以不同的角度入射到匹配層中，為體現(xiàn)兩種不同的吸收邊界對波的吸收效果，構(gòu)造了復(fù)雜低洼地質(zhì)模型，模擬現(xiàn)實(shí)地質(zhì)環(huán)境中的正、逆斷層。設(shè)置模型區(qū)域網(wǎng)格大小為300×400，模擬區(qū)域的物性參數(shù)為：上層介質(zhì)相對介電常數(shù)ε1=8，電導(dǎo)率σ1=0.003 S/m，下層介質(zhì)相對介電常數(shù)ε2=5，電導(dǎo)率σ2=0.001 S/m。所有激勵(lì)源采用主頻為900 MHz的布萊克曼—哈里斯脈沖，空間步長為0.006 m，時(shí)間步長為 0.015 ns，26個(gè)激勵(lì)源的位置沿x軸分布，每個(gè)激勵(lì)源相距5個(gè)空間步長，其模型的空間分布如圖6所示。對于遞歸卷積完全匹配層在四周設(shè)置10個(gè)網(wǎng)格大小的PML層厚度。為方便更加清晰顯示波場信息，下面將使用灰度圖顯示不同時(shí)間步Ez波場的信息。

圖6 低洼模型結(jié)構(gòu)示意Fig.6 Schematic diagram of low-lying model structure

在多個(gè)激勵(lì)源并發(fā)的情況下，采用相同主頻的脈沖電磁波，此時(shí)在近地表區(qū)域，由單個(gè)激勵(lì)源產(chǎn)生的球面波在空間上進(jìn)行波場重構(gòu)，子波波前形成平面電磁波束，如圖7a,7b所示。平面電磁波束在地下空間擴(kuò)散，在第400個(gè)時(shí)間步時(shí)遇到第二層介質(zhì)產(chǎn)生繞射波，如圖7c箭頭指示位置；第500個(gè)時(shí)間步時(shí)，低洼模型底部產(chǎn)生的反射波如圖7d標(biāo)記所示，此時(shí)反射波的能量明顯高于繞射波能量；兩種吸收邊界條件下，電磁波在介質(zhì)區(qū)域內(nèi)傳播規(guī)律相同，有相同的波場快照；在第700個(gè)時(shí)間步時(shí)，平面電磁波束到達(dá)底部邊界，基于遞歸卷積完全匹配層作為吸收邊界條件的波場快照如圖8a所示，此時(shí)平面電磁波束進(jìn)入完全匹配層被有效地吸收掉，而以二階Mur作為吸收邊界條件在底部網(wǎng)格截?cái)嗵幇l(fā)生了反射，部分波返回介質(zhì)區(qū)域，如圖8b圓圈位置，這將對能量相對較弱的繞射波造成了二次干擾，影響有效波場。

為了對比在不同偏移距下，兩種吸收邊界條件對多源并發(fā)電磁波的吸收效果，設(shè)置了寬角法觀測低洼模型，將26個(gè)發(fā)射天線(激勵(lì)源)置于低洼位置正上方，每個(gè)發(fā)射天線間隔兩個(gè)空間步長，接收天線沿x軸移動，采集200道雷達(dá)數(shù)據(jù)。天線參數(shù)依然采用900 MHz主頻的布萊克曼—哈里斯脈沖，多個(gè)發(fā)射天線無延時(shí)同時(shí)發(fā)射相同信號脈沖電磁波，得到在多源并發(fā)下，以遞歸卷積完全匹配層作為吸收邊界條件的數(shù)據(jù)記錄如圖9a所示。此時(shí)，由于采用的是相同極化方向的激勵(lì)信號，直達(dá)波、反射波與單個(gè)激勵(lì)源雷達(dá)數(shù)據(jù)記錄(如圖10a所示)相比加大了子波時(shí)寬，增強(qiáng)了來自地下界面的反射信號，此時(shí)在直達(dá)波左、右兩翼無明顯畸變，無虛假反射。而Mur二階吸收邊界條件在大偏移距處有強(qiáng)烈的虛假反射記錄；這是由于偏移距越大，波入射到匹配層的入射角度越大，掠射情況越嚴(yán)重。而在多個(gè)激勵(lì)源的情況下，在大偏移距下，虛假反射會更加明顯；這是由于加大了子波的時(shí)寬和介質(zhì)的不均一性導(dǎo)致的色散造成的。如圖9b箭頭所指示的位置，在直達(dá)波的左右翼都發(fā)生了明顯的波形畸變和虛假反射記錄。雖然在單個(gè)天線時(shí)的雷達(dá)數(shù)據(jù)記錄上，直達(dá)波左右兩翼畸變不明顯，但是也產(chǎn)生了明顯的虛假反射信號，如圖10b箭頭所指示的位置。

圖7 低洼模型不同時(shí)間步Ez波場快照Fig.7 Snapshots of the Ez wave field of the low-lying model at different time steps

a—遞歸卷積完全匹配層波場快照;b—Mur二階吸收邊界波場快照a—recursive convolution complete matching layer wave field snapshot;b—snapshot of Mur second-order absorbing boundary wave field圖8 低洼模型700時(shí)間步Ez波場快照Fig.8 Snapshot of the low-lying model at 700 time step Ez wave field

a—遞歸卷積完全匹配層數(shù)據(jù)記錄;b—Mur二階吸收邊界數(shù)據(jù)記錄a—recursive convolution complete matching layer data record;b—Mur second-order absorbing boundary data record圖9 多源并發(fā)雷達(dá)數(shù)據(jù)記錄Fig.9 Multi-source concurrent radar data recording

3 結(jié)論與展望

1)在均勻介質(zhì)下，Mur二階吸收邊界條件對多源并發(fā)下電磁波的吸收不完全，會在網(wǎng)格截?cái)嗵幃a(chǎn)生反射回波返回波場模擬區(qū)域。而基于非分裂的遞推卷積完全匹配層能很好地吸收多源并發(fā)下的電磁波，有效改善邊界位置處對干涉后的電磁波吸收效果。

2)在復(fù)雜地質(zhì)構(gòu)造下，電磁波在地下介質(zhì)中發(fā)生反射、繞射后以不同的入射角度到達(dá)網(wǎng)格截?cái)嗵帲瑐鹘y(tǒng)的Mur二階吸收邊界條件在網(wǎng)格截?cái)辔恢弥苯訉ξ⒎址匠踢M(jìn)行近似求解，不能很好地消除強(qiáng)電磁反射，在邊界網(wǎng)格截?cái)辔恢锰帟胁糠址瓷洳ǚ祷匮芯繀^(qū)域?qū)τ行Рㄔ斐啥胃蓴_，特別是在大偏移距下，會造成直達(dá)波波形畸變，形成虛假反射記錄；而基于非分裂遞歸卷積的完全匹配層作為吸收邊界條件能很好地抑制直達(dá)波波形畸變現(xiàn)象，消除虛假反射記錄。

3)基于非分裂遞歸卷積的完全匹配層作為吸收邊界條件在不分裂波場的情況下，又避免了直接對卷積的求解的復(fù)雜運(yùn)算，提升了計(jì)算的效率，改善了吸收的效果。

本文詳細(xì)介紹了非分裂遞歸卷積完全匹配層和Mur二階吸收邊界結(jié)合時(shí)域有限差分的實(shí)現(xiàn)方法，并將它們運(yùn)用到多個(gè)激勵(lì)源并發(fā)下探地雷達(dá)正演模擬，基于C語言和matlab平臺開發(fā)了相應(yīng)的數(shù)值計(jì)算程序，為陣列探地雷達(dá)的數(shù)值模擬研究提供了幫助。

致謝:感謝煤炭資源與安全開采國家重點(diǎn)實(shí)驗(yàn)室提供的計(jì)算設(shè)備，感謝評審的專家提出的寶貴意見。