基于矩陣指數(shù)法的雙線性Z變換PML模擬FDTD地下探測與成像

馮乃星 張玉賢 鄭宏興 曾慶生 童美松

（1. 深圳大學，深圳 518060；2. 河北工業(yè)大學，天津 300401；3. 南京航空航天大學，南京 210016；4. 同濟大學，上海 200092）

引 言

在現(xiàn)代通信中，空中的無線電波可以為無線電、雷達、導航系統(tǒng)和其他應(yīng)用傳送信息[1-4]. 然而，在波長較短的無線電波中發(fā)現(xiàn)了局限性[5]，其信號在傳播越來越遠時變得越來越弱，甚至從不通過水傳播，并且容易被巖層阻止. 例如，GPS[5]信號不能穿透水、土壤或建筑物墻壁，甚至不能傳播到這些物體中. 因此，海底或地下活動，如礦產(chǎn)勘查，不能采用GPS信號. 此外，GPS和其他無線電信號在室內(nèi)或摩天大樓之間的通信中表現(xiàn)出相當?shù)偷墓ぷ餍阅?

為了克服上述問題，可以將波長較長的甚低頻(very low frequency, VLF)無線電視為一種有效的方法[6]. 眾所周知，VLF波具有遠距離傳播和強大的能量穿透能力的優(yōu)點[7]，它可以在水和地球表面上傳播和跨越數(shù)百英尺，在空氣中傳播和跨越數(shù)千英里. 目前，VLF信號已廣泛應(yīng)用于軍事和民用領(lǐng)域，如海底通信、遠程通信導航、無線心率檢測、地球物理研究和救災(zāi)等. 在所有的應(yīng)用中，時域有限差分(finitedifference time-domain, FDTD)數(shù)值模擬中，對于VLF地球物理勘探的研究很少，因為對于三維VLF地下傳感的完全匹配層(perfectly matched layer,PML)的吸收特性的研究在文獻中是相當少見的.

據(jù)我們所知，F(xiàn)DTD方法作為一種高效的求解工具[8-9]，在求解Maxwell方程時受到了廣泛的關(guān)注[10]，它可以基于一個簡單的公式來避免復雜的漸近或格林函數(shù). 目前，F(xiàn)DTD已被廣泛應(yīng)用于解決天線、波導、傳播、散射等各種電磁問題[11-12]. 然而，由于計算機資源的限制，F(xiàn)DTD方法需要采用吸收邊界條件(absorbing boundary condition, ABC)來模擬任意開域問題[13]. 作為最流行和最著名的ABC之一，Bérenger于1994年首次提出的PML是一種高效的吸收材料ABC，它可以終止由非均勻、色散、各向異性甚至非線性介質(zhì)組成的FDTD問題[14].

隨著Bérenger的開創(chuàng)性工作的演變，拉伸坐標PML(stretched coordinate PML, SC-PML)和各向異性PML(uniaxial anisotropic PML, UPML)作為兩種可供選擇的PML實現(xiàn)，是由Chew等人[15]和Sacks等人分別提出的[16]，且得到了非常廣泛的應(yīng)用. 然而，如文獻[17-18]所述，這些SC-PMLs和UPMLs在減少后期反射和衰減低頻隱失波方面效果非常差.

為了改善上述傳統(tǒng)PML的缺點[19]，文獻[20]提出了通過簡單地將頻率相關(guān)的極點從實軸移到復虛半平面來實現(xiàn)復頻率偏移PML(complex frequency shifted PML, CFS-PML)的方法，其在衰減低頻隱失波和減少后期反射等方面的高效性能，引起了人們的廣泛關(guān)注. 隨后，Berenger從理論上進一步解釋了CFS-PML能夠提高隱失波吸收的原因. 此外，通過在CFS-PML方法中采用遞歸卷積(recursive convolution, RC)技術(shù)，Roden等人開發(fā)了卷積PML(convolutional PML, CPML)，其具有與傳統(tǒng)PML相似的計算效率[21]. 然而，由于在CPML公式中使用了不適當?shù)臅r間步同步，所提出的CPML整體吸收能力在某些數(shù)值情況下有些退化. 為了解決這一問題，分別基于雙線性Z變換(bilinear Z-transform, BZT)方法[22-23]、零極點匹配Z變換(matched Z-transform,MZT)方法[24]、遞推積分[25]、輔助微分方程(auxiliary differential equation, ADE)法[26]、直接Z變換(direct Z-transform, DZT)方法[27]、分段線性RC(piecewise linear RC, PLRC)[25]、時間同步[28]，以及基于積分的指數(shù)時間差分(integrate exponential time difference,IETD)[29]的幾種具有較高吸收效果的改進型PMLs.

雖然上述PML的吸收性能得到了改善，然而迭代形式復雜、效率較低. 為了解決這一問題，文獻[30]提出了一種基于矩陣指數(shù)(matrix exponential, ME)方法的CFS-PML實現(xiàn)方案，該方案具有較高的精度，并為平面倒F天線(planar inverted-F antenna, PIFA)的應(yīng)用保持了良好的效率. 然而，此文獻中只有一個PIFA案例用于反映和分析其性能，所以我們不知道是否也可以計算其他問題并獲得良好的結(jié)果. 因此，為了展示基于ME的PML實現(xiàn)在不同應(yīng)用中的通用性，我們考慮了更多的情況.

本文提出了一種基于ME方法的CFS格式的高效、準確的BZT-PML方法，并將其應(yīng)用于VLF地下傳感問題. 另一方面，除了地下傳感之外，還利用甚高頻(very high frequency, VHF)地下成像問題來驗證所提公式的完整性和穩(wěn)健性. 該提案稱為BZTME-PML. 由于采用了BZT方法，頻域的乘法和時域的卷積均對應(yīng)于Z域的乘法. 該方案結(jié)合了ME和Z變換方法的優(yōu)點，既能保持較好的吸收精度，又能保持較好的效率，而且可以方便、直接地離散Maxwell方程組.

1 PML算法推導

在三維SC-PML區(qū)域中，頻域Maxwell旋度方程可以表示如下：

如式(1)與(2)所示，電通量D和磁通量B本構(gòu)關(guān)系式可以被應(yīng)用于更新電場與磁場，優(yōu)點在于可以避免考慮材料的具體屬性[31]，定義如下：

式中，εr(ω)和μr(ω)分別是FDTD計算域的相對介電常數(shù)和導磁率.

式(1)和(2)中的操作算子定義為

式中，Sη(η=x,y,z)是n階復數(shù)拉伸坐標變量，且有

為了節(jié)省更多的內(nèi)存，我們采用一階CFS-PMLSη(η=x,y,z)，表達式為

式中：κη是大于等于1的實數(shù)；ση和αη均被假定為正實數(shù). 接下來，重新對式(7)進行整理和簡化，如下：

接著，將式(8)的倒數(shù)轉(zhuǎn)換到s域，結(jié)合關(guān)系式j(luò)ω→s, 采用如下BZT方法：

可以得到

式中：至此，我們已經(jīng)取得了復數(shù)拉伸變量的Z域表達式.

接下來，以x場分量為例，在笛卡爾坐標系下展開式(1)如下：

因此，在Z域中有

將式(12)代入式(15)，并引入輔助變量ψ，有

經(jīng)過整理，我們得到：

為了更加簡易地采用ME方法[30,32]，需要進一步整理式(17)，并在式(18)和(19)等式左邊添加ψDxy?ψDxz，然后再在兩邊乘以1/Δt，有：

式中：

因此可以總結(jié)得到不同場量的表達通式為

式中，Θ和Φ為FDTD迭代過程的電場和磁場的通用符號. 當Θ=E時Φ=H，當Θ=H時Φ=E，我們可以取得所有的矢量方程，如下式：

式中：

假定初始條件為t0，我們可以取得

通過采用上一時間步nΔt取代t0，式(27)關(guān)于新的時間步(n+1)Δt可表達如下：

再次簡化式(28)，有

式中，

在FDTD計算中，在第nth時間步，PML區(qū)域的實現(xiàn)步驟如下：

第一步 邊界預(yù)存儲場

第二步 電通量場迭代(以x場分量為例)

第三步 更新輔助變量

目前為止，我們已經(jīng)推導得到了完整的數(shù)學表達式. 為了能夠闡述所提出算法BZT-ME-PML的有效性和高效性，另外七種方法也被實現(xiàn)用以對比，分別是ME-PML、Conventional PML、DZT-PML、CPML、ADE-PML、BZT-PML，以及MZT-PML.

為了驗證本文提出的方法在PML上的數(shù)值色散問題，在三維的數(shù)值色散關(guān)系中，必然滿足

通過對正弦函數(shù)泰勒展開可以得到

考慮到一般網(wǎng)格在各坐標軸的方向上都取相同的厚度δ，即Δx= Δy= Δz=δ，cΔt= 0.5δ，此時

假定各方向的數(shù)值波長為λdη=λ0δ?1(η=x,y,z)，將式(32)中的高階項消去可得

式中：

不妨令數(shù)值色散關(guān)系的特征矩陣為

發(fā)現(xiàn)在PML層中必然受到空間波矢(kx,ky,kz)的影響，在全局電磁空間的FDTD計算可以與矩陣K3構(gòu)成數(shù)值精度上的正相關(guān)關(guān)系，在側(cè)面上也驗證了在PML-ABC層下的FDTD方法與無吸收邊界的FDTD方法具有相同的數(shù)值色散關(guān)系. 因此，可以給出三維情形下的FDTD關(guān)于數(shù)值波長λdη(η=x,y,z)的評價函數(shù)為

由此得到表1的列表.

表1 λdz = 25時FDTD評價函數(shù)ζ3 (λdx, λdy, 25)的取值Tab. 1 The numerical value of evaluation function ζ3 (λdx, λdy, 25) of FDTD at λdz = 25

當各方向的數(shù)值波長均大于25時，PML-ABC下的三維FDTD方法可以保持數(shù)值誤差在2.867 6×10?3以下，也肯定了本文提出的PML算法能夠有效收斂且計算穩(wěn)定.

2 數(shù)值算例實驗

2.1 TE極化電磁波與完純導體片的相互作用問題

為了驗證本文所提出算法的有效性，我們采用8G內(nèi)存i7-8750H CPU的計算機進行算法驗證. 首先仿真一個二維數(shù)值算例. 如圖1所示，一個TE極化的電磁波與一個在z方向上無窮大而在x方向上有限長度的完純導體片的相互作用，空間步長是Δx=Δy=1 mm，時間步長是Δt=1.178 5 ps ≈0.5×CFL ，其中CFL表示FDTD方法離散時間長度的穩(wěn)定因子，一個寬度為100個元胞的PEC被自由空間所包圍，10個元胞厚度的PML截斷自由空間，在各個方向上PML與PEC的距離是3個元胞.

圖1 二維仿真中的FDTD模型Fig. 1 The FDTD model in 2D simulation

一個在z方向上無窮長的y方向極化的線電流源被放置在中心，并且被一個微分高斯脈沖所激勵[33-34]，其中tw=26.53 ps，t0=4tw. 觀察點“P”選取在PEC薄片的邊上. 參考解網(wǎng)格充分大以至于在1 500個時間步內(nèi)沒有來自外邊界的反射波. PML參數(shù)設(shè)定如下：對于傳統(tǒng)PML，κmax=9，σmax= 0.6σopt；對于CFS-PML，κmax= 10，σmax= 0.9σopt，αη= 0.05.

我們提出的方法與無吸收邊界的傳統(tǒng)FDTD進行對比得到的最大相對反射誤差(maximum relative reflection error, MRRE)公式為

式中：x為我們提出的方法所獲取到的場數(shù)據(jù)；xref為無吸收邊界的傳統(tǒng)FDTD的場數(shù)據(jù).

相對反射誤差(relative reflection error，RRE)曲線的計算結(jié)果如圖2所示. 另外，對不同算法中的MRRE對應(yīng)內(nèi)存和計算時間(時間步=20 000)進行統(tǒng)計，結(jié)果如表2所示. 圖2和表2表明，二維算例結(jié)果驗證了所提算法的有效性. 從圖2看出，所提出的算法比首次提出的基于輔助微分方程的ME-PML在吸收精度上有些許提高，比其他的PML方法在吸收精度上有明顯提高.

圖2 二維FDTD情況下不同時間步數(shù)下的RRE曲線Fig. 2 The RRE curve under the different time steps in 2D FDTD cases

表2 二維FDTD情況下不同PML算法實現(xiàn)的性能對比Tab. 2 The computation comparison in those different PML methods in 2D FDTD cases

2.2 金屬薄板嵌入在有耗電介質(zhì)中的三維模型

為了進一步驗證所提出算法的有效性，一個100 mm× 25 mm的金屬薄板(被仿真為完純導體)嵌入在本構(gòu)參數(shù)為ε和σ的有耗電介質(zhì)中的問題被采用[23]. 有耗電介質(zhì)的本構(gòu)參數(shù)被假設(shè)為εr=7.73和σ=0.273. 三維FDTD計算域(圖3)是106× 31× 6，空間步長和時間步長分別是Δx=Δy=Δz=1 mm和Δt=5.3 ps. 在每個方向上再用厚度為10個元胞的PML去截斷FDTD計算域. 在PML區(qū)域，ση和κη采用PML內(nèi)部中本構(gòu)參數(shù)分布，其中m=4；而αη是常數(shù)，即不會隨著空間位置發(fā)生變化. 一個在垂直方向上極化的激勵源Js被放置于金屬薄板中的左上角，其中，tw=83 ps，t0=4tw. 觀察點“P”放置于金屬薄板的對角線上.

圖3 三維仿真中的FDTD模型Fig. 3 The FDTD model in 3D simulation

PML參數(shù)設(shè)定如下：對于傳統(tǒng)PML，κmax= 14，σmax= 0.56σopt；對于CFS-PML，κmax=11，σmax= 0.4σopt，αη= 0.05.

由圖4和表3可知，在處理三維問題時，所提出的算法的有效性得以進一步驗證.

圖4 三維FDTD情況下不同時間步數(shù)下的RRE曲線Fig. 4 The RRE curve under the different time steps in 3D FDTD cases

表3 三維FDTD情況下不同PML算法實現(xiàn)的性能對比Tab. 3 The computation comparison in those different PML methods in 3D FDTD cases

2.3 三維VLF地下探測問題

目前為止，二維和三維算例已經(jīng)被用于驗證所提出算法的有效性和穩(wěn)定性. 接下來，為了進一步驗證所提出算法的優(yōu)勢：高精度(ME技術(shù))、易操作和粗空間網(wǎng)格設(shè)置(BZT變換)，兩個三維VLF地下探測和一個三維VHF地下成像問題被采用.

眾所周知，一個高效的ABC應(yīng)該被使用截斷帶有高損耗媒質(zhì)的FDTD計算區(qū)域，用以模擬開域空間問題，使得電磁波相互作用的解可以被精確計算.為此，一個三維色散半空間問題被實現(xiàn)，如圖5所示，可以看到三個礦體分布在不同的內(nèi)部區(qū)域.

圖5 三維情形下的內(nèi)嵌立方礦體的半空間幾何模型Fig. 5 The half-space geometry with three different cubic ores in 3-D cases

整個FDTD計算區(qū)域分別包含空氣、圍巖(電導率σrock=0.005 S/m)，以及礦體(電導率σore=4 S/m). 三個立方體礦體(大小為268 m× 268 m× 268 m)內(nèi)嵌在巖石區(qū)域，且它們的幾何中心位置坐標分別為(2 010, 1 608, 670) m、(2 010, 2 010, 804) m、(2 010, 2 412,670) m. 空間離散元胞大小為Δx= Δy= Δz= 67 m. 整個FDTD計算區(qū)域大小為4 020 m× 4 020 m×2 010 m. 激勵源采用的是磁點偶極子源，位于(2 010,2 010, 1 742) m的空氣中. 接收點位于(2 077, 2 077,1 809) m.如圖6，一個頻率為25 kHz的雙極性方波脈沖被采用，是因為該激勵是適合應(yīng)用于航空瞬變電磁系統(tǒng)問題. 如圖7所示，橫坐標為0.16 ms的時間窗被用于展示所提出算法與其他算法在接收點處的吸收性能.

圖6 頻率為25 kHz的雙極性方波脈沖(上升期是0.005 ～0.015 ms，下降期是0.025 ～ 0.035 ms)Fig. 6 Bipolar square wave pulse with frequency of 25 kHz(rising period at time = 0.005?0.015 ms, falling at time = 0.025?0.035 ms)

圖7 在內(nèi)嵌立方礦體的半空間幾何模型下的RRE曲線Fig. 7 The RRE curves under the half-space geometry with three different cubic ores

為了觀察比較MRRE，一個參考解需要被提供，即在原先計算模型的大小上向各個方向拓展100個元胞，最外層采用相同個數(shù)的PML截斷計算區(qū)域.觀察點處的MRRE可以通過式(30)計算得出.

由圖7可知，所提出算法與其他算法的MRRE分別是ADE-PML (?49.19 dB)、MZT-PML (?49.74 dB)、BZT-PML (?49.98 dB)、DZT-PML (?62.90 dB)、MEPML (?81.29 dB)，以及BZT-ME-PML (?82.09 dB).

為了更加符合真實環(huán)境中不規(guī)則的地形分布，地面不規(guī)則模型被模擬，其中三個礦體也是不規(guī)則地分別內(nèi)嵌在巖石層中，如圖8所示，整個計算區(qū)域大小為4 020 m× 4 020 m× 2 010 m，而地下巖石的探測區(qū)域為1 608 m. 激勵源采用的是電線源，坐落在巖石層表面，脈沖仍然是雙極性方波脈沖. 空間元胞大小為Δx= Δy= Δz= 67 m.

圖8 三維情形下的不規(guī)則的半航空瞬變電磁幾何模型Fig. 8 The electromagnetic geometry for irregular semiaeronautical transients in 3-D cases

采用同圖7一樣的時間窗進行觀察，結(jié)果如圖9所示. 可以看到，不同PML算法的MRRE分別是ADE-PML (?95.81 dB)、MZT-PML (?95.71 dB)、BZTPML (?95.16 dB)、DZT-PML (?109.2 dB)、ME-PML(?113.2 dB)、BZT-ME-PML (?114 dB).

圖9 在不規(guī)則半航空瞬變電磁幾何模型下的RRE曲線Fig. 9 The RRE curves under the electromagnetic geometry for irregular semi-aeronautical transients

2.4 三維VHF地下成像問題

為表明所提出算法的魯棒性和完整性，一個三維VHF地下成像問題被考慮. 據(jù)我們所知，嫦娥5號探測工程項目即將發(fā)射升空對月球表面進行探測，并在不久的將來采樣回地球. 因此，在地球上進行高分辨率地下成像是非常有必要的. 在過去的研究中，月球表面土壤的相對介電常數(shù)默認為εr= 3[35-37]. 為了模擬月球環(huán)境，我們采用相對介電常數(shù)為εr= 3的火山灰替代月球表面土壤作為嵌入目標體的背景，以此在地球上進行月球探測的實測實驗.

如圖10所示，來自實測月球探測系統(tǒng)嫦娥五號的兩個單站實驗數(shù)據(jù)被提取，填充了火山灰的地下實驗?zāi)Ｐ痛笮?.5 m × 2.0 m × 1.0 m.

圖10 嫦娥五號提供的兩個單站實驗室數(shù)據(jù)采集系統(tǒng)Fig. 10 The data-acquisition systems of two single station provided by Chang’e-5

在實測實驗中，激勵信號從嫦娥五號超寬帶雷達系統(tǒng)中發(fā)射，并在多道天線中記錄. 在嫦娥五號實測結(jié)束之后，所提出的BZT-ME-PML算法被應(yīng)用于實現(xiàn)地下成像，如圖11 所示. 因此，這也進一步證明了所提出算法的有效性和穩(wěn)定性.

圖11 嫦娥五號超寬帶系統(tǒng)實測數(shù)據(jù)的逆時偏移成像結(jié)果Fig. 11 Reverse time migration imaging results of Chang’ E-5 UWB system measured data

3 結(jié)束語

利用ME方法，所提出的緊湊型BZT-ME-PML算法能夠高效地求解VLF地下探測和VHF地下成像問題. 所提出的算法具有避免進行公式重整理、卷積操作和變量替換等優(yōu)勢，使得計算過程變得高效.此外，由于采用了BZT方法，無論是頻域相乘還是時域卷積都可以對應(yīng)到Z域相乘，使得所提出的算法具有易操作性. 因此，同時結(jié)合ME技術(shù)和Z變換方法，所提出的BZT-ME-PML算法不僅具有較好的吸收性能和高效率，而且易于直接離散Maxwell旋度方程. 最后，理論推導及算例表明該方法有效可行的.