柳尚斌，孫懷鳳,3*，李文翰，王震，楊洋,3

1 山東大學巖土與結(jié)構(gòu)工程研究中心，濟南 250061 2 山東大學地球電磁探測研究所，濟南 250061 3 山東省工業(yè)技術(shù)研究院先進勘探與透明城市協(xié)同創(chuàng)新中心，濟南 250061 4 山東大學北京研究院，北京 100086

0 引言

時域有限差分(FDTD)方法最早由Yee(1966)提出，通過在時間和空間上對電場、磁場交替采樣，用中心差分對Maxwell的兩個旋度方程離散，得到電磁場的顯式迭代方程.中心差分具有二階精度，但需要時間步長Δt嚴格滿足CFL條件才能保證數(shù)值穩(wěn)定性(Wand and Hohmann，1993).嚴苛的穩(wěn)定性條件使顯式FDTD需要大量的時間迭代步才能模擬晚期瞬變電磁響應，過多的迭代次數(shù)不僅會增加模擬時間，同時也會增加晚期的累積誤差.Crank-Nicolson時域有限差分法(CN-FDTD)(Fijany et al.，1995；Yang et al.，2006)是一種對任意時間步長都無條件穩(wěn)定的算法，由于時間步長不再受CFL條件限制，可以采用比ΔtCFL(顯式FDTD算法滿足CFL條件時的最大時間步長)大的多的時間步長，從而減少時間迭代步數(shù).但由于該算法的電場、磁場采樣時刻相同，在差分離散后的求解方程中電、磁場分量相互耦合在一起，每一時間步電磁場的計算都需要求解一個大型稀疏矩陣，不管采用直接法還是迭代法求解，都不可避免地會占用大量內(nèi)存或時間.因此，相較于原始FDTD算法，CN-FDTD計算效率的提升非常有限.之后，為了提高計算效率，一些學者相繼提出了一系列基于CN-FDTD算法的近似策略，例如用于二維的approximate-decoupling策略(Sun and Trueman，2004)和用于三維的direct-splitting(DS)策略(Sun and Trueman，2006)，以及交變方向隱式時域有限差分方法(ADI-FDTD)(Garcia et al.，2002)等.由于這些近似算法不再求解大型稀疏矩陣，而是求解低階且主對角占優(yōu)的三對角矩陣，因此可以用快速算法求解(Thomas，1949)，顯著的降低內(nèi)存消耗和計算時間.ADI-FDTD以及CNDS-FDTD算法的無條件穩(wěn)定性驗證在多篇文獻(Zheng and Chen，2001；Tan，2008；Jiang et al.，2019)中都有提及.

然而，Yang和Tan(2017)實驗發(fā)現(xiàn)，ADI-FDTD算法的穩(wěn)定性并不是完全沒有條件，該算法的穩(wěn)定前提是在均勻時間步長下，當時間步長不均勻時，ADI-FDTD算法不再穩(wěn)定.本文通過試驗發(fā)現(xiàn)，CNDS-FDTD算法也存在類似的問題，即在非均勻時間步長下，該算法變的和普通FDTD算法一樣，Δt需滿足CFL條件，否則，算法很快會發(fā)散.而在瞬變電磁三維正演中，時間步的大小并不是均勻的，在早期會采用較小的Δt，使得電磁響應能夠反映近地表的電性特征；隨著時間的迭代，電磁波中的高頻成分迅速衰減，低頻逐漸成為主要頻段，更大的波長使得沒有必要再采用小的Δt，因此會采用更大的Δt來節(jié)省模擬時間.根據(jù)上面的分析可知，瞬變電磁正演中的時間步長是非均勻的，這意味著如將CNDS-FDTD算法用于瞬變電磁三維正演，由于Δt無法突破CFL限制，正演的效率并不會提高.

Backward Euler(BE)差分格式同樣對任意時間步長都無條件穩(wěn)定，在地球物理數(shù)值模擬中有著廣泛的應用(Yang et al.，2014；惠哲劍等，2020；齊彥福等，2021).然而，BE差分后的Maxwell方程同樣是隱式求解格式，與CN-FDTD方法一樣，都面臨著求解大型稀疏矩陣的問題.為了提升求解效率，本文將DS近似策略引入到BE差分格式后的FDTD方程中，形成了BEDS-FDTD新方法.DS策略將電磁場分量解耦，并將大型稀疏矩陣降階和重構(gòu)為低階且主對角占優(yōu)的三對角矩陣，可以采用Thomas算法快速計算.BE方法在時間上是穩(wěn)定的，但加入DS策略是否會對其穩(wěn)定性產(chǎn)生影響，以及BEDS-FDTD在有耗介質(zhì)、非均勻時間步長下的穩(wěn)定性尚未見報道.借鑒Zheng和Chen(2001)對FDTD類算法無條件穩(wěn)定性證明的方法，采用von Neumann方法求幅值增長系數(shù)，判斷算法的穩(wěn)定性.在實踐中發(fā)現(xiàn)，在有耗介質(zhì)、非均勻時間步長下的場幅值增長系數(shù)的解析解非常復雜，對于這種未知量過多的符號類函數(shù)求解，是一個非常耗時的過程，作者在一臺CPU為Intel Core i5-7400的PC機上采用Matlab求解48 h，仍不收斂.因此，本文采用半解析的方式，首先減少未知量個數(shù)，將電導率、時間步長給定，然后采用大量隨機抽樣的方式來計算幅值增長系數(shù)，進而對BEDS-FDTD算法的穩(wěn)定性進行驗證.實驗結(jié)果驗證了全新算法BEDS-FDTD在有耗介質(zhì)、均勻和非均勻時間步長下都是穩(wěn)定的，這為將BEDS-FDTD算法用于瞬變電磁三維正演提供了理論前提.

1 理論

1.1 PML與普通介質(zhì)統(tǒng)一控制方程

電磁波的傳播本身是一個開域問題，但在數(shù)值模擬時，有限的計算資源無法模擬電磁波傳播的整個過程，因此需要在模型某處截斷，人為添加邊界條件.Dirichlet邊界條件是最常采用的邊界條件，其實質(zhì)是在模型最外側(cè)添加一層理想電導體(Perfect Electric Conductor，PEC)，電磁波遇到PEC會完全反射.該邊界條件施加方式簡單，只需在邊界網(wǎng)格處強制令電場的切向分量為0即可.然而，如果模型尺寸不夠大，由邊界反射產(chǎn)生的“偽反射波”會對觀測區(qū)域的計算結(jié)果產(chǎn)生嚴重影響，而采用過大的模型尺寸又會增加模擬時間.完全匹配層(Perfectly Matched Layer，PML)吸收邊界是為了解決電磁場遠距離傳播所設計的，電磁波在PML介質(zhì)內(nèi)傳播時會迅速衰減，而且在界面處無反射.通常情況下，PML介質(zhì)遠離發(fā)射源所在區(qū)域，那么在均勻、各向同性、非磁性媒質(zhì)中的頻率域Maxwell旋度方程(Chew and Weedon，1994)為

(1)

(2)

其中，η代表笛卡爾坐標系下坐標方向，且滿足循環(huán)移位規(guī)律，例如當η=x時，η-1和η+1分別代表z和y方向，ε為介電常數(shù)，μ為真空磁導率.S為復數(shù)形式的坐標伸縮因子，在CFS-PML介質(zhì)中，為了改善對低頻電磁波的吸收效果，將Sη定義(Roden and Gedney，2000)為

(3)

式中，κη為網(wǎng)格延拓因子，σp η是PML介質(zhì)中人為添加的電導率，αη是一個大于零的實數(shù).

將方程從頻率域轉(zhuǎn)換到時間域，往往通過傅里葉反變換，這不可避免的會引入卷積項，而卷積的計算通常費時費力.為了避免復雜的卷積運算，本文采用雙線性變換方法，首先利用關(guān)系jω→s將Sη變換到s域，則

(4)

其中，uη=αη/ε0，vη=uη+σp η/(κηε0).

然后采用雙線性變換關(guān)系s=(2/Δt)·(1-z-1)/(1+z-1)將(4)式變換到Z域：

(5)

其中，

將(5)式代入到方程(1)和(2)中，并將兩個方程等號左側(cè)的jω用其Z域的表達式(1-z-1)/Δt替換，可得到Maxwell旋度方程在Z域下的形式

(6)

(7)

式中，系數(shù)w、φ和θ下標中的e代表的是求解電場，h代表的是求解磁場，但計算公式是相同的.將方程(6)和(7)化簡：

(8)

(9)

其中，a0=ε/(ε+σΔt)，a1=Δt/(ε+σΔt)，b=Δt/μ.如果令

(10)

將式(10)展開并化簡，可得到Ψe(η,η+1)的計算公式：

(11)

同理，如果令

(12)

可得到Ψe(η,η-1)的計算公式：

(13)

將式(10)和(12)代入到方程(8)中，并考慮式(11)和(13)，整理得

Eη=a0z-1Eη+(Ψe(η,η+1)-θe(η+1)z-1Ψe(η,η+1))-(Ψe(η,η-1)-θe(η-1)z-1Ψe(η,η-1))

(14)

如果令reη=φeη-θeη，并將方程(14)在n+1時刻差分離散，同時考慮到Z變換的時移特性，經(jīng)整理，得到時域方程為

(15)

其中，D代表一階中心差分算子.

(16)

(17)

同理，方程(9)經(jīng)過相似的變換，可得到下式

(18)

其中，

(19)

(20)

經(jīng)過上述的處理過程，方程(15)和(18)在時間上為BE差分格式，下面進行簡單證明.將時間域Maxwell旋度方程中場對時間的一階導數(shù)用BE差分近似，可得

(21)

(22)

將上面兩個方程(21)和(22)寫成分量形式，并化簡：

(23)

(24)

其中，η、a0、a1和b與上文中的定義一致.

分析方程(15)、(18)，當w和Ψ分別取值1和0后，即將PML介質(zhì)內(nèi)的迭代方程退化到普通介質(zhì)內(nèi)的迭代格式時，它們與方程(23)、(24)完全一致.因此，我們通過選擇合適的變換關(guān)系得到的方程(15)和(18)是在時間上滿足BE差分格式的電、磁場方程.

1.2 引入DS策略的近似算法

若要對方程組(15)和(18)求解，可將(18)代入(15)中，使得新方程的左側(cè)不再包含n+1時刻的H，化簡后，可得

(25)

其中，D2代表二階中心差分算子.

由方程(25)可以看到，等號的左側(cè)是n+1時刻的場值(待求量)，等號的右側(cè)是n時刻的場值(已知量)，可通過求解線性方程組來計算待求場值.但在等號左側(cè)，由于電場的不同方向、不同空間位置的分量相互耦合在一起，使得需要處理的系數(shù)矩陣是一個大型稀疏矩陣，這種類型的矩陣無論是采用直接法還是迭代法求解，都不可避免地會占用大量內(nèi)存和時間.為了提升求解效率，本文將DS策略引入到求解方程中，將方程(15)和(18)寫成下列矩陣的形式：

(26)

在方程(26)的左側(cè)加入高階誤差項AB(Wn+1-Wn)，整理后得

注意到方程(27)等號左側(cè)可做因式分解，因此求解過程可分為以下兩個子步驟：

(28)

(I6-B)Wn+1=W*-BWn，

(29)

其中，場值的上標“*”代表非真實的采樣時刻，是為了方便計算而引入的中間過程.將算子矩陣A、B和C代入方程(28)和(29)中并化簡，可得

(30)

(31)

En+1=a1MTHn+1+E*-a1MTHn，

(32)

Hn+1=bMEn+1+H*-bMEn.

(33)

通過一系列變換可使得方程(30)、(32)和(33)中等號的右側(cè)不包含H*和Hn+1，然后將矩陣M代入并整理，得

(34)

(35)

(36)

式中，c=a1b.

方程(34)—(36)即為BEDS-FDTD算法在PML介質(zhì)中的控制方程，其求解過程為：首先通過方程(34)計算E*，再通過方程(35)計算n+1時刻的電場，最后由方程(36)顯式迭代計算n+1時刻的磁場.需要說明的是，這里計算的磁場不是真實的物理場，而是一個中間量，若要得到真實的磁場，需將電場值代入式(18)中求解.BEDS-FDTD算法在普通介質(zhì)中的控制方程和在PML介質(zhì)中的控制方程，均可寫成(34)—(36)的形式，當w和Ψ分別取值1和0后，該控制方程就退化到普通介質(zhì)內(nèi)的格式.

1.3 系數(shù)矩陣快速求解

以方程(34)為例，當η=x時，通過z方向的坐標逐行變化(x和y坐標不變)，方程可寫為

(37)

(38)

(39)

(40)

其中，2≤k≤n.那么系數(shù)矩陣P1為

矩陣P是滿秩矩陣，且只有三個主對角上有非零元素，考慮到α、β和γ的表達式滿足不等式條件：|β2|>|γ2|>0，|βk|≥|αk|+|γk|，αk、γk≠0，k=3,4,…,n-1；|βn|>|αn|>0.因此，矩陣P是主對角占優(yōu)的三對角矩陣，可以采用Thomas算法快速求解.

2 穩(wěn)定性驗證

當在瞬變電磁正演中對復雜地形模擬時，地空邊界采用向上延拓不再適用，此時不可避免地需要將空氣層引入到計算內(nèi)，理想的空氣是無耗介質(zhì)(σ=0 S·m-1).但在實際的算法設計中，不會讓空氣的電導率真正的為0，而是給一個足夠小的值(例如，σ=10-6S·m-1)，此時，空氣由無耗介質(zhì)轉(zhuǎn)為有耗介質(zhì)；除了空氣之外，大地本身同樣為有耗介質(zhì).當算法用于瞬變電磁三維正演時，需要考慮這種使用場景下的特殊條件，包括時間步長非均勻(參考引言中的內(nèi)容)、介質(zhì)有耗，所以下文主要討論這種特殊情況下算法的數(shù)值穩(wěn)定性.

本文采用von Neumann方法求幅值增長因子來判斷算法的數(shù)值穩(wěn)定性(陳娟等，2016)，如果計算的場幅值增長系數(shù)均小于或等于1，那么說明該算法是穩(wěn)定的；否則，說明該算法是不穩(wěn)定的.由于在瞬變電磁正演中的時間步長不均勻，只考慮從n時刻到n+1時刻的場值增長情況并不能反映算法在非均勻時間步長下的穩(wěn)定性，需要增大觀測時間，考慮從n時刻到n+2時刻的場值增長情況.BEDS-FDTD算法的最根本求解方程為方程(27)(之后的一系列數(shù)學變換都是基于該方程)，為了簡化計算，下面僅考慮普通介質(zhì)中的方程

(41)

在線性、各向同性介質(zhì)中傳播的平面波方程為

(42)

(43)

當η=x時，將函數(shù)ψ對空間的導數(shù)進行中心差分離散：

(44)

同理，當η=y,z時，也有類似的表達式：

Dyψ(x,y,z)=τyψ(x,y,z)，

(45)

Dzψ(x,y,z)=τzψ(x,y,z)，

(46)

將式子(42)—(46)代入方程(41)中，得到BEDS-FDTD從n時刻到n+1時刻的場值計算公式

(47)

同理，BEDS-FDTD算法從n+1時刻到n+2時刻的場值計算公式可寫做

(48)

式中，

綜合考慮公式(47)和(48)，從而得到BEDS-FDTD算法從n時刻到n+2時刻的迭代公式：

(49)

式中，右上角的“-1”代表矩陣求逆運算.考慮平面波方程(42)和(43)，可以得到Wn+2=ξ2Wn；那么，方程(49)經(jīng)整理、化簡為

(50)

若方程(50)具有非零解，系數(shù)矩陣的行列式必須為零，所以：

(51)

通過求解式(51)可獲得幅值增長系數(shù)ξ的解析解，然而，由于式(51)中的未知量過多，很難求得ξ解析解的表達式.本文采用半解析的方式，用大量隨機抽樣的方式計算ξ.首先減少未知量的個數(shù)，令空間網(wǎng)格的尺寸為10 m，空氣電導率σ=10-6S·m-1.若Δtn=αΔtCFL，Δtn+1=βΔtCFL，那么可記作CFLN(n,n+2)=(α,β)；在Matlab R2020上對kx、ky和kz各隨機抽樣取值1000次，將計算的|ξ|繪制在復平面的上半平面上.為了方便對比，將采用同樣方法計算的CN-FDTD和CNDS-FDTD算法的結(jié)果同樣繪制在圖中.圖1a、1b展示了非均勻時間步長下的結(jié)果，在圖1a中非均勻時間步長為CFLN(n,n+2)=(2,8)，從圖1a可以看到，原始CN-FDTD算法的|ξ|位于單位圓的圓周上，這意味著CN-FDTD算法是穩(wěn)定的，但近似算法CNDS-FDTD卻在某些波數(shù)時失去了穩(wěn)定性.考慮到相鄰Δt變化的過于劇烈，在圖1b中展示了當CFLN(n,n+2)=(2,3)時的結(jié)果，圖1b的結(jié)果與圖1a中的現(xiàn)象是一致的，CNDS-FDTD算法同樣發(fā)散了.作為對比，本文提出的全新算法BEDS-FDTD在CFLN(n,n+2)=(2,8)和CFLN(n,n+2)=(2,3)時的|ξ|都位于單位圓內(nèi)，這說明BEDS-FDTD在不均勻時間步長下是穩(wěn)定的.圖1c和1d中展示了均勻時間步長下三種算法的結(jié)果，可以看到，三種算法在均勻時間步長下都是穩(wěn)定的.

圖1 從n時刻到n+2時刻的場值增長系數(shù)ξ

從上面的分析可以看出，盡管原始CN-FDTD算法在均勻和非均勻時間步長下都是穩(wěn)定的，但近似后的CNDS-FDTD算法卻只能在均勻時間步長下維持穩(wěn)定，在非均勻時間步長下不再穩(wěn)定.而本文提出的新算法BEDS-FDTD，不管是在均勻時間步長還是非均勻時間步長下均能保持數(shù)值穩(wěn)定.上面的實驗結(jié)果為新算法BEDS-FDTD應用到瞬變電磁三維正演中提供了理論前提.

3 模型試驗

3.1 BEDS-FDTD算法計算效率與精度測試

為了檢驗新算法BEDS-FDTD用于瞬變電磁三維正演的適用性和效率，本文將其用于典型層狀模型的模擬中，選擇的第一個測試模型為H型模型，模型參數(shù)細節(jié)如表1所示.空氣的電導率為σ=10-6S·m-1，大地與空氣的相對介電常數(shù)均為εr=7.8.110 m×110 m的矩形回線敷設在空氣與大地之間作為發(fā)射源，電流為1 A，發(fā)射波形采用梯形波(孫懷鳳等，2013)，波形1的具體參數(shù)如圖2所示，除此之外，本文還設計了波形2，與波形1的區(qū)別是Turn-off時間由10-6s改為10-7s.模型(包括PML介質(zhì)區(qū)域)采用均勻網(wǎng)格剖分，網(wǎng)格尺寸采用兩種方案，策略1是網(wǎng)格在x、y和z方向上的尺寸為22 m×22 m×20 m；策略2是網(wǎng)格在x、y和z方向上的尺寸為10 m×10 m×10 m.PML介質(zhì)的參數(shù)按照下式選?。?/p>

表1 模型參數(shù)

圖2 激勵波形示意圖

(52)

(53)

(54)

其中，d是PML介質(zhì)的厚度；x是PML介質(zhì)中的某點距離PML與正常介質(zhì)的分界面的距離，m和ma為多項式的階數(shù)，影響著參數(shù)的變化速度，通常情況下ma=1，m=3；本文經(jīng)過大量的模擬測試后發(fā)現(xiàn)，在參數(shù)σmax=1.0×10-2S·m-1，αmax=1.0×10-4S·m-1和κmax=1.0時，CFS-PML邊界能取得較好的吸收效果.

BEDS-FDTD算法求解過程中，最耗時的部分為時間步迭代過程，其中包含了大量的多重循環(huán)，這種結(jié)構(gòu)非常適合并行加速.本文算法的代碼基于GPU并行計算的OpenACC框架編寫，實驗中的GPU采用專門用于高性能計算(High performance computing，HPC)的專業(yè)卡Tesla A100.由于BEDS-FDTD算法使用了吸收邊界，通常不需要考慮太多的網(wǎng)格數(shù)目，但更多的網(wǎng)格可以模擬更加精細的模型，因此本次實驗將網(wǎng)格數(shù)量作為變量測試算法的效率，所有模型的網(wǎng)格數(shù)均包括外部10層的CFS-PML介質(zhì)，模型剖分采用策略1.圖3給出了BEDS-FDTD算法和常規(guī)FDTD算法在不同網(wǎng)格數(shù)量時的計算效率對比.可以看到BEDS-FDTD算法在模擬網(wǎng)格數(shù)為50×50×50(圖3中計作503)的模型(計算域是30×30×30)時，整個計算過程共需要16000個迭代時間步，用時僅10 s，而即使網(wǎng)格數(shù)增加到200×200×200(圖3中計作2003)，用時也僅224 s，計算效率非常高.作為對比，同樣采用了CFS-PML邊界以及GPU并行的FDTD算法的整個計算過程共需要233000個時間迭代步，需要注意的是，該FDTD算法已經(jīng)考慮了虛擬介電常數(shù)(孫懷鳳等，2013)，即在確保傳導電流占據(jù)主導時，時間步長相較于普通FDTD已經(jīng)有了很大的增長.但由于迭代時間步數(shù)仍遠超BEDS-FDTD，F(xiàn)DTD在計算不同尺寸的模型時，平均用時要比BEDS-FDTD多9倍左右.

圖3 BEDS-FDTD和FDTD算法在計算不同數(shù)量網(wǎng)格時的效率

為了驗證BEDS-FDTD算法的精度，我們將采用不同波形、不同網(wǎng)格剖分策略的H模型(網(wǎng)格數(shù)都為50×50×50，計算域為30×30×30)的數(shù)值解與一維數(shù)字濾波解(李貅，2002)做了對比，觀測時間為電流關(guān)斷后10-5～10-2s，計算結(jié)果如圖4所示.圖4a中是數(shù)值解與數(shù)字濾波解的對比，圖4b為相對誤差曲線.可以看到，當BEDS-FDTD算法采用網(wǎng)格剖分策略1時，采用波形2的解比采用波形1的解在早期精度更高；而在晚期，兩者的精度是一致的，這和文獻(Zeng et al.，2019)中的結(jié)論一致，即更短的關(guān)斷時間會使得模擬結(jié)果的早期精度更高；而采用網(wǎng)格剖分策略2和波形2的BEDS-FDTD的解精度最高，在觀測時間內(nèi)，與數(shù)字濾波解的相對誤差都在5%以下，滿足模擬計算要求.在圖4中還繪制了同樣采用網(wǎng)格剖分策略2和波形2的FDTD的解，可以看出，整體的誤差和BEDS-FDTD算法的誤差基本上處在同一水平，但值的關(guān)注的是，晚期(電流關(guān)斷3 ms之后)FDTD的誤差有繼續(xù)擴大的趨勢，而且其晚期的誤差值要大于BEDS-FDTD的誤差.本次實驗表明，不僅電流的關(guān)斷時間會對算法的早期精度產(chǎn)生明顯的影響，網(wǎng)格的大小同樣會對計算結(jié)果的早期精度產(chǎn)生顯著影響.而且，實驗結(jié)果也證明了本文針對全新算法BEDS-FDTD提出的CFS-PML邊界非常有效，在較小的模型尺度下(采用網(wǎng)格剖分策略2時，包括PML介質(zhì)的整個模型尺寸為500 m×500 m×500 m)也能獲得精度非常高的結(jié)果.

圖4 H型模型的BEDS-FDTD數(shù)值解與一維數(shù)字濾波解的對比

為了進一步驗證計算精度，本文還將BEDS-FDTD算法用于模擬K型層狀模型，參數(shù)如表1所示，模型(包括PML介質(zhì))采用第二種剖分策略(網(wǎng)格大小為10 m)，網(wǎng)格數(shù)量為50×50×50，其中PML介質(zhì)在各個方向上占據(jù)10個網(wǎng)格，因此，計算域大小為30×30×30.發(fā)射源的參數(shù)與H型模型的一致，發(fā)射電流波形采用波形2，計算結(jié)果展示在圖5中.從圖中可以看出，除了早期外，數(shù)值解與一維數(shù)字濾波解均吻合的較好，而早期差異較大的原因是由于BEDS-FDTD發(fā)射波形2雖然對階躍波進行了模擬，但實際上只是修改后的梯形波，關(guān)斷效應導致了與采用階躍波響應的數(shù)字濾波解在早期差異較大，而隨著關(guān)斷效應的逐漸減弱，在大部分觀測時間內(nèi)，二者的相對誤差都不超過5%.

圖5 K型層狀模型的BEDS-FDTD數(shù)值解與一維線性數(shù)字濾波解的對比

在上面的精度驗證模型中的激勵源都是回線源，下面的算例中用A型層狀模型驗證算法在電性源激勵下的計算精度.模型的具體參數(shù)如表1所示，采用網(wǎng)格剖分策略2剖分，模型在x、y和z方向上的網(wǎng)格數(shù)為50×50×50，各個方向最外面的10層為PML介質(zhì)(PML介質(zhì)的參數(shù)與上面的模型一致)，因此，計算域的網(wǎng)格數(shù)目為30×30×30.發(fā)射源為110 m長的接地線源，發(fā)射電流為1 A，發(fā)射波形采用波形2.接收點距離線源的中心110 m，圖6a繪制了接收點處的感應電動勢隨時間的變化曲線，圖6b展示了數(shù)值解與數(shù)字濾波解的相對誤差.可以看到，模擬結(jié)果與數(shù)字濾波解吻合的較好，誤差在大部分時間內(nèi)均小于5%，僅在極早期大于5%，具體原因在K型模型中進行了詳細說明，這里不再贅述.

圖6 A型層狀模型的BEDS-FDTD數(shù)值解與一維線性數(shù)字濾波解的對比

3.2 復雜模型實驗

為了測試BEDS-FDTD算法的適用性，本文將其用于復雜三維模型的計算，模型的具體細節(jié)如圖7所示.空氣的電導率為σ=10-6S·m-1，大地與空氣的相對介電常數(shù)均為εr=7.8.大地表面是一厚50 m，電導率為0.1 S·m-1的覆層，薄層下面存在著一個形狀非常復雜的三維低阻體(在z方向上呈階梯狀分布)，電導率為1 S·m-1，該低阻體將大地分成電導率為0.01 S·m-1和0.003333 S·m-1的兩部分.100 m長的導線敷設在如圖7所示位置，發(fā)射電流為1 A，發(fā)射波形采用波形2.PML介質(zhì)的參數(shù)與層狀模型選取的一致.整個模型在x、y和z方向上為900 m×2050 m×1200 m，采用均勻網(wǎng)格剖分，網(wǎng)格在x、y和z方向的尺寸為10 m×12.5 m×10 m.計算域為70Δx×144Δy×100Δz，最外側(cè)包裹了10層的PML介質(zhì).

圖7 復雜三維模型示意圖

坐標的中心位于發(fā)射線源的中心，圖8繪制了三個觀測點處的感應電動勢的結(jié)果.圖中曲線突變的地方為數(shù)據(jù)從正值轉(zhuǎn)為負值，從圖中可以看到，BEDS-FDTD的結(jié)果與FDTD(Commer et al.，2004)和MFVTD(Liu et al.，2019)的結(jié)果一致性非常好，在具體細節(jié)上，本文提出的BEDS-FDTD的結(jié)果與FDTD的結(jié)果吻合的程度更高，在整個觀測時間內(nèi)，基本看不出明顯的偏差.

圖8 復雜模型的計算結(jié)果

4 結(jié)論

本文提出了瞬變電磁三維正演新算法BEDS-FDTD，該算法在有耗介質(zhì)、均勻和非均勻時間步長下均能維持無條件穩(wěn)定性，而且其時間步長不再受CFL穩(wěn)定性條件限制，可有效減少迭代次數(shù).DS策略的引入將大型稀疏矩陣轉(zhuǎn)變?yōu)橐幌盗械碗A、主對角占優(yōu)的三對角矩陣，可快速求解.實驗發(fā)現(xiàn)，采用GPU并行后，用Tesla A100模擬50×50×50規(guī)模的模型，精度能滿足計算要求的情況下僅用時10 s.即使模型規(guī)模擴大到200×200×200，用時也不超過4 min.之后，將新算法BEDS-FDTD用于不同復雜程度的模型計算中，計算結(jié)果驗證了該算法具有較高的精度.模型實驗證明了BEDS-FDTD算法是穩(wěn)定、高效的，可為瞬變電磁三維快速反演提供正演基礎.

致謝本文在孫懷鳳等(2013)開發(fā)的TEM3DFDTD程序基礎上完成.作者在研究過程中就時間步長的選擇與算法精度的問題與長安大學的齊彥福老師進行了深入討論；審稿專家對本文寫作的提高提出了很多建設性的意見，在此一并表示感謝.