孫懷鳳， 柳尚斌， 楊洋

1 山東大學(xué) 巖土與結(jié)構(gòu)工程研究中心， 濟(jì)南 250061 2 山東大學(xué) 地球電磁探測研究所， 濟(jì)南 250061

0 引言

瞬變電磁法(TEM)利用不接地回線或接地長導(dǎo)線向大地發(fā)射一次場，在電流關(guān)斷間隙，通過采集到的感應(yīng)電動勢分析地下電性分布，在工程勘察、金屬礦產(chǎn)資源勘探、環(huán)境與水文地質(zhì)調(diào)查、采空區(qū)探測等領(lǐng)域得到廣泛應(yīng)用(Chang et al., 2019a; Chang et al., 2019b; 薛國強(qiáng)等,2007,2008).瞬變電磁三維正演是進(jìn)行勘探可行性研究、方案設(shè)計(jì)、觀測數(shù)據(jù)反演與解釋的必備工具.國際上，瞬變電磁或時間域電磁三維正演研究主要有五類方法，即體積分方程法(VIE)(Cox et al., 2010; Xiong, 1992; Zhdanov et al., 2006)、有限單元法(FEM)(Sugeng, 1998; Li et al., 2017; Li et al., 2011)、時域有限差分法(FDTD)(孫懷鳳等, 2013; 許洋鋮等, 2012; Commer et al., 2004; Wang et al., 1993)、有限體積法(FV)(Zhou et al., 2018; Ren et al., 2018; Oldenburg et al., 2013)以及譜Lanczos分解法(Spectral Lanczos Decomposition Method, SLDM)(Druskin et al., 1999).瞬變電磁正演的過程就是給定地電模型通過數(shù)值方法求解Maxwell方程組，獲得隨時間衰減的電磁響應(yīng)曲線.盡管求解方程的過程與數(shù)學(xué)方法各不相同，但瞬變電磁的基本方程是相同的.

在時域有限差分瞬變電磁正演方面：Oristaglio和Hohmann(1984),閆述等(2002)采用Du Fort-Frankel有限差分方法研究了二維地電斷面平行導(dǎo)線電流源產(chǎn)生的瞬變電磁場，并分析了均勻半空間中包含二維異常體時的瞬變電磁場分布特征.Wang和Hohmann(1993)采用改進(jìn)的Du Fort-Frankel方法首次給出了通過求解一階Maxwell方程組進(jìn)行電磁探測建模的三維時域有限差分算法.之后，宋維琪和仝兆歧(2000)針對電偶源瞬變電磁進(jìn)行了三維有限差分法正演計(jì)算.Commer和Newman(2004;2005)針對電性源長偏移瞬變電磁(LOTEM)建模問題建立了求解非因果場的三維時域有限差分并行算法.針對井下瞬變電磁探測的問題，岳建華和楊海燕等(2008)進(jìn)行了礦井瞬變電磁探測三維時域有限差分正演研究.孟慶鑫、潘和平(2012)采用三維時域有限差分研究了地井和井中瞬變電磁響應(yīng).

自2012年以來，瞬變電磁時域有限差分三維正演成為國內(nèi)的研究熱點(diǎn)，許洋鋮和林君等(2012)的航空時間域電磁響應(yīng)三維有限差分?jǐn)?shù)值計(jì)算，邱稚鵬和黃清華等(2013)的起伏地形下瞬變電磁場的三維建模方法，李建慧和胡祥云等(2013)基于電場Helmholtz方程的回線源瞬變電磁三維正演方法，孫懷鳳和李貅等(2013)考慮關(guān)段時間的回線源瞬變電磁三維正演方法，余翔和王緒本等(2017)的三維時域有限差分的CPML邊界以及趙越和李貅等(2017)針對航空電磁開展的復(fù)雜模型三維正演研究.之后，孫懷鳳等(2018)還開發(fā)了針對瞬變電磁三維時域有限差分網(wǎng)格剖分的多尺度網(wǎng)格方法，用于解決小目標(biāo)建模問題并提高正演速度.

然而，常規(guī)的時域有限差分方法為了保證數(shù)值計(jì)算的穩(wěn)定，時間步長Δt需要嚴(yán)格滿足Courant-Friedrich-Lewy(CFL)限制條件，這使得迭代計(jì)算過程中時間采樣過密、計(jì)算時間過長.為了克服上述問題， Fijany等(1995)首次提出采用Crank-Nicolson差分格式的有限差分方法求解麥克斯韋方程.Crank-Nicolson方法是一種隱式差分方法，具有無條件穩(wěn)定的特點(diǎn)，由于Δt不再受CFL的限制，可以有效的減少迭代次數(shù).然而，該方法的缺點(diǎn)也非常明顯，即在每次迭代中都需要求解大型稀疏矩陣，當(dāng)矩陣階數(shù)較大時會占用大量計(jì)算資源和時間，這制約了該方法的推廣和應(yīng)用.直到近年來，CN-FDTD近似快速算法陸續(xù)出現(xiàn)，才使得CN-FDTD方法得到逐步推廣，主要的近似求解算法有：CNDS-FDTD(Sun et al., 2003)，CNAFS-FDTD(Sun et al., 2004)，CNCSU-FDTD(Sun et al., 2006)等，這些方法保留了CN-FDTD方法的無條件穩(wěn)定性，求解速度得到了巨大提升，同時，計(jì)算精度遠(yuǎn)高于ADI等隱式方法(Garcia et al., 2002).

本研究基于Crank-Nicolson差分方法對Maxwell方程組重新離散，空間網(wǎng)格仍然采用Yee元胞，時間步進(jìn)采用整時間步電場、磁場同時采樣，建立無條件穩(wěn)定的FDTD格式.在時間采樣上，與常規(guī)FDTD交替采樣相比，CN-FDTD電場、磁場同時采樣，構(gòu)成了隱式差分格式，需要求解稀疏矩陣方程組.采用Crank-Nicolson-Cycle-Sweep-Uniform(CNCSU)近似求解方法，在保證精度的同時，計(jì)算效率大幅提高，且內(nèi)存占用小.

由于計(jì)算資源是有限的，無法模擬電磁波在開域情況下的傳播過程，因此需要將計(jì)算模型在某處截?cái)?常規(guī)的瞬變電磁數(shù)值模擬中，邊界條件多采用Dirichlet邊界，這是因?yàn)樵撨吔鐥l件理論簡單、容易實(shí)現(xiàn)，但如果想要邊界反射對目標(biāo)區(qū)域的擴(kuò)散場影響足夠小，則需要將邊界設(shè)置的足夠遠(yuǎn)，采用更大尺度的模型，不可避免的會增大計(jì)算量.因此，一些學(xué)者開始研究其他的邊界條件來代替?zhèn)鹘y(tǒng)的Dirichlet邊界，如Mur吸收邊界(Mur, 1981)、 Liao吸收邊界(Liao et al., 1984)和完全匹配層邊界(Berenger, 1994; Berenger, 1996)等.由于復(fù)頻率參數(shù)完全匹配層(CFS-PML)(Kuzuoglu et al., 1996)對低頻波有較好的吸收效果，李展輝和黃清華(2014)以及余翔和王緒本等(2017)將CFS-PML應(yīng)用于瞬變電磁FDTD的數(shù)值模擬中.本文采用雙線性變換方法將CFS-PML邊界條件施加于無條件穩(wěn)定的CN-FDTD方法，并將算法用于模擬三維瞬變電磁場的傳播.

圖1 CN-FDTD使用的Yee網(wǎng)格Fig.1 Yee grid

1 CN-FDTD方法

1.1 CN-FDTD時空采樣

圖2 常規(guī)FDTD與CN-FDTD電場、 磁場時間采樣分布對比 (a) 常規(guī)FDTD時間采樣； (b) CN-FDTD時間采樣.Fig.2 Comparison of original FDTD and CN-FDTD sampling in time for electric and magnetic fields (a) Original FDTD； (b) CN-FDTD.

CN-FDTD是基于克蘭克-尼科爾森(Crank-Nicolson)差分格式提出的時域有限差分方法，該方法是一種既能夠滿足時間步長放大又能夠滿足計(jì)算精度且無條件穩(wěn)定的數(shù)值計(jì)算方法(Fijany et al., 1995; Sun et al., 2003; Yang et al., 2006)，在微波、毫米波以及光學(xué)領(lǐng)域等高頻電磁波問題中已經(jīng)得到應(yīng)用.與常規(guī)FDTD一樣，CN-FDTD仍然采用Yee元胞作為基本空間離散單元，如圖1所示，電場在網(wǎng)格棱邊中心采樣，磁場在網(wǎng)格各面的中心采樣.然而，在時間采樣上，CN-FDTD與常規(guī)FDTD存在較大差別，圖2給出了CN-FDTD與常規(guī)FDTD時間步進(jìn)采樣對比圖，與常規(guī)FDTD的電場和磁場在n和n+1/2時刻交替采樣不同，CN-FDTD的電場和磁場同時在整數(shù)時刻采樣.

1.2 控制方程與差分離散

在均勻、各向同性、非磁性、無源媒質(zhì)中， Maxwell方程組旋度方程可以寫成如下分量形式(葛德彪和閆玉波，2005)：

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

其中，H和E分別為磁場強(qiáng)度和電場強(qiáng)度，ε、σ和μ分別是介電常數(shù)、電導(dǎo)率和磁導(dǎo)率，x、y、z構(gòu)成直角坐標(biāo)系.

Crank-Nicolson差分策略在空間域使用中心差分，以式(1a)為例，用n和n+1時刻的平均值代替式中出現(xiàn)的n+1/2時刻的場值,可得：

(2a)

(2b)

(2c)

(2d)

將(2a—2d)代入(1a)并化簡可得：

(3)

類似地，可以得到其他離散方程，經(jīng)過進(jìn)一步壓縮與簡化可以概化為

(4a)

(4b)

其中，η(η=x,y,z)為坐標(biāo)方向，滿足循環(huán)移位規(guī)律，即若η=x, 則η-1和η+1分別對應(yīng)z和y；Dη為沿η軸方向的一階中心差分算子，系數(shù)a2=Δt/2μ.

將(4b)代入(4a) 并代換n+1時刻的磁場分量，可得：

(5)

式中，D2η為沿著η軸方向的二階中心差分算子，例如η=x時，公式(5)中的二階差分項(xiàng)可以表示為

(6)

從(5)式可以看出，等號的左側(cè)為n+1時刻的電場值，為待求的未知量，等號右側(cè)全部為n時刻的已知場值，通過求解3個方程組就能得到待求電場.然而，等號左側(cè)的三個電場分量相互耦合在一起，求解每一時間步上的電場分量都需要求解一個大型稀疏矩陣，這將會消耗大量的計(jì)算資源與時間，計(jì)算效率較低(Feng et al., 2018; Sun et al., 2004).

1.3 CNCSU近似及CFS-PML吸收邊界

事實(shí)上，在CNCSU-FDTD近似公式的推導(dǎo)過程中，發(fā)現(xiàn)迭代公式在CFS-PML介質(zhì)內(nèi)和普通介質(zhì)內(nèi)具有類似的形式，因而可以推導(dǎo)更一般的CFS-PML介質(zhì)內(nèi)的迭代公式，將相應(yīng)的系數(shù)替換就可以得到CNCSU-FDTD在普通介質(zhì)內(nèi)的求解迭代公式.

為了獲得復(fù)頻率參數(shù)下的PML邊界條件，在均勻、各向同性、非磁性無源媒質(zhì)中頻率域麥克斯韋旋度方程寫為(Chew et al., 1994)

(7a)

(7b)

其中S為坐標(biāo)伸縮因子，在CFS-PML介質(zhì)中，Roden等(2000)將S定義為

(8)

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

通常情況下，F(xiàn)DTD計(jì)算會將公式(7)、(8)轉(zhuǎn)換到時間域來求解，然而，從頻率域通過傅里葉反變換到時間域會出現(xiàn)卷積項(xiàng)，為避免三維模型中的復(fù)雜卷積運(yùn)算，可將公式(7)、(8)進(jìn)行雙線性變換(Ramadan and Oztoprak, 2002)，將其轉(zhuǎn)換到Z域，整理得

(9a)

(9b)

(9c)

(9d)

考慮Z變換的時移特性，可將(9)式變換為

(10a)

(10b)

(10c)

(10d)

將(10b)代入(10a)并替換n+1時刻的磁場，在η=x,y,z坐標(biāo)方向上可以表示為

(11a)

(11b)

(11c)

公式(11)可以簡化表示為：

(12)

合并整理并進(jìn)行因式分解，可得下述子步驟公式：

(13a)

(13b)

(14a)

(14b)

(14c)

(14d)

(14e)

(14f)

2 精度驗(yàn)證與復(fù)雜模型對比

為了驗(yàn)證CN-FDTD瞬變電磁三維正演方法的計(jì)算精度以及對復(fù)雜模型的適用性，與均勻半空間模型解析解、三層模型線性數(shù)字濾波解進(jìn)行了對比，之后，選擇三維垂直接觸帶復(fù)雜模型并與三維矢量有限元、三維有限體積以及常規(guī)時域有限差分的計(jì)算結(jié)果進(jìn)行對比.

2.1 簡單模型對比與精度驗(yàn)證

本節(jié)計(jì)算模型統(tǒng)一使用中心回線裝置，發(fā)射回線為100 m×100 m，激發(fā)電流為1 A，發(fā)射波形采用文獻(xiàn)(孫懷鳳等, 2013)中開關(guān)函數(shù)處理后的梯形波，上升沿與下降沿均為1 μs，電流持續(xù)時間為30 ms，觀測時間為30 ms，占空比為1∶1，相當(dāng)于實(shí)際觀測中雙極性矩形波條件下瞬變電磁基頻為8.33 Hz.如圖3a所示，模型分為計(jì)算區(qū)域與PML邊界區(qū)域，計(jì)算區(qū)域采用20 m均勻網(wǎng)格剖分，網(wǎng)格數(shù)為61×61×30，計(jì)算區(qū)域模型尺寸為1220 m×1220 m×600 m，其中，空氣層采用6層網(wǎng)格，大地采用24層網(wǎng)格；PML邊界區(qū)域采用15層，包裹在計(jì)算區(qū)域之外.

均勻半空間模型大地電阻率為100 Ωm，空氣電阻率設(shè)為106Ωm，數(shù)值模擬結(jié)果與解析解(李貅, 2002; Anderson, 1979)對比時間范圍為關(guān)斷后1 μs～30 ms，圖3b給出了三維CN-FDTD計(jì)算結(jié)果與階躍波解析解的感應(yīng)電動勢衰減曲線對比.可以看到，在關(guān)斷后的早期二者差異較大，這是由于三維計(jì)算中考慮了關(guān)斷時間而不是階躍關(guān)斷引起的(孫懷鳳等, 2013)，而隨著關(guān)斷效應(yīng)的消失，誤差逐漸減小，10 μs之后兩者吻合較好.

為了進(jìn)一步測試，采用典型的三層模型，并與線性數(shù)字濾波解(李貅, 2002; Anderson, 1979)進(jìn)行對比，層狀模型的層厚和電阻率如表1所示，計(jì)算時所使用的其他參數(shù)與均勻半空間模型一致.圖4給出了四種層狀模型的感應(yīng)電動勢衰減曲線的對比結(jié)果.從圖中可以看出，層狀模型的對比結(jié)果與均勻半空間模型類似，均為早期誤差較大，這是由于關(guān)斷時間的影響，隨著時間的推移，關(guān)斷時間的影響逐漸減弱，誤差逐漸減小，在晚期，數(shù)值解與解析解曲線吻合較好，二者幾乎重合.

2.2 復(fù)雜模型對比與精度驗(yàn)證

采用三維垂直接觸帶復(fù)雜模型(Commer et al., 2004)進(jìn)行進(jìn)一步驗(yàn)證.模型參數(shù)如圖5所示，地下介質(zhì)可以分為四部分，表層厚50 m，電阻率為10 Ωm，左側(cè)與右側(cè)部分的電阻率分別為100 Ωm 和300 Ωm，二者界面上有一個形狀復(fù)雜的接觸帶，其電阻率為1 Ωm，接觸帶長400 m.發(fā)射回線邊長為100 m，設(shè)置四個接收點(diǎn)，以最左側(cè)的源邊為相對坐標(biāo)起點(diǎn)，四個接收點(diǎn)坐標(biāo)分別為(0,50,0)、(0,150,0)、(0,450,0)和(0,1050,0).模型采用均勻網(wǎng)格進(jìn)行剖分，網(wǎng)格大小為25 m，計(jì)算區(qū)域尺寸為1500 m×1500 m×750 m(網(wǎng)格數(shù)60×60×30)，其中，空氣層層數(shù)為6，大地層數(shù)為24，CFS-PML介質(zhì)層數(shù)為15.

表1 層狀模型參數(shù)Table 1 Parameters of the layered models

圖3 均勻半空間模型CN-FDTD與解析解對比 (a) 模型示意圖; (b) CN-FDTD數(shù)值解與解析解的感應(yīng)電動勢衰減曲線對比.Fig.3 Comparison of CN-FDTD numerical solution and analytical solution for Homogeneous half space model (a) Schematic diagram of the model; (b) Decay curve comparison of the CN-FDTD and analytic solutions.

圖4 典型三層模型CN-FDTD數(shù)值解與線性數(shù)字濾波解對比 (a) A型； (b) H型； (c) K型； (d) Q型.Fig.4 Comparison of digital filter solution and CN-FDTD for typical there-layered models (a) A type Model; (b) H type Model; (c) K type Model; (d) Q type Model.

圖5 三維垂直接觸帶復(fù)雜模型 (Li et al., 2017)Fig.5 Three-dimensional complex model with a vertical contact zone

圖6 三維垂直接觸帶復(fù)雜模型計(jì)算結(jié)果對比Fig.6 The numerical results of three-dimensional complex model

圖6給出了CN-FDTD計(jì)算的感應(yīng)電動勢絕對值隨時間的變化曲線，在曲線的尖點(diǎn)處數(shù)據(jù)符號發(fā)生改變，在圖中給出了標(biāo)示，并同時給出了和矢量有限元法(Li et al., 2017)、有限體積法(Zhou et al., 2018)、以及多尺度網(wǎng)格FDTD(孫懷鳳等，2018)方法的對比結(jié)果，發(fā)現(xiàn)僅在感應(yīng)電動勢變號處有明顯的差異，其他區(qū)域吻合較好.

3 CN-FDTD時間步放大對計(jì)算精度的影響

圖7 CN-FDTD在不同CFLNs下結(jié)果的對比Fig.7 Comparison of CN-FDTD numerical results under different CFLNs

對比試驗(yàn)在一臺CPU為Intel Core i5-7300 HQ的筆記本電腦上完成，且所有的對比試驗(yàn)均使用單線程計(jì)算，未啟用任何并行策略.表2給出了不同CFLN時的計(jì)算時間和內(nèi)存開銷，可以看出，隨著CFLN的增加，迭代次數(shù)和計(jì)算時間總體減小，但由于模擬過程考慮了發(fā)射波形，激勵源采用梯形波，上升沿以及下降沿非常陡峭，因此在計(jì)算過程中，為了保證結(jié)果穩(wěn)定性，上述階段的Δt采用逐漸增大和逐漸減小的操作(見圖8).經(jīng)過測試，在梯形波持續(xù)時間的前半段以及關(guān)斷后Δt的增長速率都各有一個閾值，不能增長過快，如果超過這個閾值，模擬結(jié)果與解析解會產(chǎn)生較大的誤差，這使得即使CFLN成倍數(shù)增長，而迭代次數(shù)以及計(jì)算時間并不是同比例的減小.從圖8可以看出，由于關(guān)斷后Δt增長較慢，當(dāng)CFLN≥800，計(jì)算結(jié)束時，Δt還沒有增長到預(yù)設(shè)的大小，因此，對整體迭代步數(shù)貢獻(xiàn)主要集中在on-time階段.如果不考慮關(guān)斷時間，而是采用(Wang and Hohmann, 1993)給出的設(shè)定初始場的方法施加激勵源，則計(jì)算時間會大幅減小.表2顯示CFLN=1600時，計(jì)算時間減少到了50 min，后續(xù)如果考慮FDTD的天然并行性，上述模型的計(jì)算時間有望減少到3 min以內(nèi)，這將為后續(xù)的三維反演工作帶來便利條件.而且，針對所設(shè)計(jì)的模型所需的內(nèi)存消耗僅為192 Mb，在普通PC就可以輕松完成.即使使用較多的網(wǎng)格模擬非常精細(xì)的地下目標(biāo)時，也能夠保持較小的內(nèi)存開銷.

表2 不同CFLNs下CN-FDTD算法內(nèi)存開銷及計(jì)算時間Table 2 Memory overhead and computing time of CN-FDTD algorithm under different CFLNs

圖8 CN-FDTD在不同CFLNs下Δt的變化對比Fig.8 Changes of Δt in the CN-FDTD under different CFLNs

4 結(jié)論

常規(guī)FDTD進(jìn)行瞬變電磁三維正演時，迭代時間步長受CFL穩(wěn)定性條件限制，造成迭代次數(shù)過多、計(jì)算時間過長、晚期累積誤差增大.本研究提出CN-FDTD隱式差分方法用于瞬變電磁三維正演，突破了CFL條件的限制，顯著減少迭代步數(shù)，通過CNCSU-FDTD近似算法實(shí)現(xiàn)方程組快速求解，通過雙線性變換方法施加CFS-PML吸收邊界條件，降低計(jì)算時間成本.通過與均勻半空間解析解、層狀模型的線性數(shù)字濾波解，以及三維垂直接觸帶復(fù)雜模型的矢量有限元、有限體積法、常規(guī)FDTD解進(jìn)行了對比與精度驗(yàn)證，結(jié)果均吻合較好.本研究給出的方法內(nèi)存消耗非常小，可以使用大量網(wǎng)格模擬非常復(fù)雜的模型，計(jì)算過程可以在普通PC上快速完成.

致謝研究在孫懷鳳等(2013)開發(fā)的tem3dfdtd程序基礎(chǔ)上完成.作者在研究過程中與長安大學(xué)李貅教授團(tuán)隊(duì)魯凱亮博士進(jìn)行了深入討論，在此一并表示感謝.