孫懷鳳， 柳尚斌， 楊洋

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

0 引言

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

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

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

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

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

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

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

1 CN-FDTD方法

1.1 CN-FDTD時(shí)空采樣

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

1.2 控制方程與差分離散

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

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

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

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

(2a)

(2b)

(2c)

(2d)

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

(3)

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

(4a)

(4b)

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

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

(5)

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

(6)

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

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

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

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

(7a)

(7b)

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

(8)

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

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

(9a)

(9b)

(9c)

(9d)

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

(10a)

(10b)

(10c)

(10d)

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

(11a)

(11b)

(11c)

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

(12)

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

(13a)

(13b)

(14a)

(14b)

(14c)

(14d)

(14e)

(14f)

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

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

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

本節(jié)計(jì)算模型統(tǒng)一使用中心回線裝置，發(fā)射回線為100 m×100 m，激發(fā)電流為1 A，發(fā)射波形采用文獻(xiàn)(孫懷鳳等, 2013)中開(kāi)關(guān)函數(shù)處理后的梯形波，上升沿與下降沿均為1 μs，電流持續(xù)時(shí)間為30 ms，觀測(cè)時(shí)間為30 ms，占空比為1∶1，相當(dāng)于實(shí)際觀測(cè)中雙極性矩形波條件下瞬變電磁基頻為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)對(duì)比時(shí)間范圍為關(guān)斷后1 μs～30 ms，圖3b給出了三維CN-FDTD計(jì)算結(jié)果與階躍波解析解的感應(yīng)電動(dòng)勢(shì)衰減曲線對(duì)比.可以看到，在關(guān)斷后的早期二者差異較大，這是由于三維計(jì)算中考慮了關(guān)斷時(shí)間而不是階躍關(guān)斷引起的(孫懷鳳等, 2013)，而隨著關(guān)斷效應(yīng)的消失，誤差逐漸減小，10 μs之后兩者吻合較好.

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

2.2 復(fù)雜模型對(duì)比與精度驗(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，二者界面上有一個(gè)形狀復(fù)雜的接觸帶，其電阻率為1 Ωm，接觸帶長(zhǎng)400 m.發(fā)射回線邊長(zhǎng)為100 m，設(shè)置四個(gè)接收點(diǎn)，以最左側(cè)的源邊為相對(duì)坐標(biāo)起點(diǎn)，四個(gè)接收點(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與解析解對(duì)比 (a) 模型示意圖; (b) CN-FDTD數(shù)值解與解析解的感應(yīng)電動(dòng)勢(shì)衰減曲線對(duì)比.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ù)字濾波解對(duì)比 (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é)果對(duì)比Fig.6 The numerical results of three-dimensional complex model

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

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

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

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

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

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

4 結(jié)論

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

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