瞬變電磁Crank-Nicolson FDTD三維正演

孫懷鳳， 柳尚斌， 楊洋

1 山東大學 巖土與結構工程研究中心， 濟南 250061 2 山東大學 地球電磁探測研究所， 濟南 250061

0 引言

瞬變電磁法(TEM)利用不接地回線或接地長導線向大地發(fā)射一次場，在電流關斷間隙，通過采集到的感應電動勢分析地下電性分布，在工程勘察、金屬礦產(chǎn)資源勘探、環(huán)境與水文地質(zhì)調(diào)查、采空區(qū)探測等領域得到廣泛應用(Chang et al., 2019a; Chang et al., 2019b; 薛國強等,2007,2008).瞬變電磁三維正演是進行勘探可行性研究、方案設計、觀測數(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方程組，獲得隨時間衰減的電磁響應曲線.盡管求解方程的過程與數(shù)學方法各不相同，但瞬變電磁的基本方程是相同的.

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

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

然而，常規(guī)的時域有限差分方法為了保證數(shù)值計算的穩(wěn)定，時間步長Δt需要嚴格滿足Courant-Friedrich-Lewy(CFL)限制條件，這使得迭代計算過程中時間采樣過密、計算時間過長.為了克服上述問題， Fijany等(1995)首次提出采用Crank-Nicolson差分格式的有限差分方法求解麥克斯韋方程.Crank-Nicolson方法是一種隱式差分方法，具有無條件穩(wěn)定的特點，由于Δt不再受CFL的限制，可以有效的減少迭代次數(shù).然而，該方法的缺點也非常明顯，即在每次迭代中都需要求解大型稀疏矩陣，當矩陣階數(shù)較大時會占用大量計算資源和時間，這制約了該方法的推廣和應用.直到近年來，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)定性，求解速度得到了巨大提升，同時，計算精度遠高于ADI等隱式方法(Garcia et al., 2002).

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

由于計算資源是有限的，無法模擬電磁波在開域情況下的傳播過程，因此需要將計算模型在某處截斷.常規(guī)的瞬變電磁數(shù)值模擬中，邊界條件多采用Dirichlet邊界，這是因為該邊界條件理論簡單、容易實現(xiàn)，但如果想要邊界反射對目標區(qū)域的擴散場影響足夠小，則需要將邊界設置的足夠遠，采用更大尺度的模型，不可避免的會增大計算量.因此，一些學者開始研究其他的邊界條件來代替?zhèn)鹘y(tǒng)的Dirichlet邊界，如Mur吸收邊界(Mur, 1981)、 Liao吸收邊界(Liao et al., 1984)和完全匹配層邊界(Berenger, 1994; Berenger, 1996)等.由于復頻率參數(shù)完全匹配層(CFS-PML)(Kuzuoglu et al., 1996)對低頻波有較好的吸收效果，李展輝和黃清華(2014)以及余翔和王緒本等(2017)將CFS-PML應用于瞬變電磁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)差分格式提出的時域有限差分方法，該方法是一種既能夠滿足時間步長放大又能夠滿足計算精度且無條件穩(wěn)定的數(shù)值計算方法(Fijany et al., 1995; Sun et al., 2003; Yang et al., 2006)，在微波、毫米波以及光學領域等高頻電磁波問題中已經(jīng)得到應用.與常規(guī)FDTD一樣，CN-FDTD仍然采用Yee元胞作為基本空間離散單元，如圖1所示，電場在網(wǎng)格棱邊中心采樣，磁場在網(wǎng)格各面的中心采樣.然而，在時間采樣上，CN-FDTD與常規(guī)FDTD存在較大差別，圖2給出了CN-FDTD與常規(guī)FDTD時間步進采樣對比圖，與常規(guī)FDTD的電場和磁場在n和n+1/2時刻交替采樣不同，CN-FDTD的電場和磁場同時在整數(shù)時刻采樣.

1.2 控制方程與差分離散

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

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

其中，H和E分別為磁場強度和電場強度，ε、σ和μ分別是介電常數(shù)、電導率和磁導率，x、y、z構成直角坐標系.

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

(2a)

(2b)

(2c)

(2d)

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

(3)

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

(4a)

(4b)

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

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

(5)

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

(6)

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

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

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

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

(7a)

(7b)

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

(8)

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

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

(9a)

(9b)

(9c)

(9d)

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

(10a)

(10b)

(10c)

(10d)

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

(11a)

(11b)

(11c)

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

(12)

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

(13a)

(13b)

(14a)

(14b)

(14c)

(14d)

(14e)

(14f)

2 精度驗證與復雜模型對比

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

2.1 簡單模型對比與精度驗證

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

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

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

2.2 復雜模型對比與精度驗證

采用三維垂直接觸帶復雜模型(Commer et al., 2004)進行進一步驗證.模型參數(shù)如圖5所示，地下介質(zhì)可以分為四部分，表層厚50 m，電阻率為10 Ωm，左側與右側部分的電阻率分別為100 Ωm 和300 Ωm，二者界面上有一個形狀復雜的接觸帶，其電阻率為1 Ωm，接觸帶長400 m.發(fā)射回線邊長為100 m，設置四個接收點，以最左側的源邊為相對坐標起點，四個接收點坐標分別為(0,50,0)、(0,150,0)、(0,450,0)和(0,1050,0).模型采用均勻網(wǎng)格進行剖分，網(wǎng)格大小為25 m，計算區(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ù)值解與解析解的感應電動勢衰減曲線對比.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 三維垂直接觸帶復雜模型 (Li et al., 2017)Fig.5 Three-dimensional complex model with a vertical contact zone

圖6 三維垂直接觸帶復雜模型計算結果對比Fig.6 The numerical results of three-dimensional complex model

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

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

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

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

表2 不同CFLNs下CN-FDTD算法內(nèi)存開銷及計算時間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 結論

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

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