戴世坤，陳輕蕊*，凌嘉宣，李昆，趙東東，張瑩，張錢江

1 中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室，長沙 410083 2 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083 3 西南石油大學(xué)，成都 610500 4 桂林電子科技大學(xué)電子工程與自動(dòng)化學(xué)院，桂林 541004 5 桂林理工大學(xué)，桂林 541004

0 引言

大地電磁法利用天然電磁場進(jìn)行地下電性結(jié)構(gòu)的勘探，由于其野外工作簡單、成本低、受野外地質(zhì)條件約束小和探測深度廣等優(yōu)點(diǎn)被廣泛應(yīng)用于礦產(chǎn)資源勘查、深部地質(zhì)構(gòu)造研究和工程和環(huán)境勘查等領(lǐng)域（王家映,1997;Pedersen et al.,2006;閆永利等,2007;Tezkan and Saraev,2008;Becken et al.,2011;徐光晶等,2015）.正演是反演的基礎(chǔ)，高效高精度的正演是提高反演速度和效果的關(guān)鍵，也是準(zhǔn)確進(jìn)行實(shí)測數(shù)據(jù)資料處理和地質(zhì)地球物理解釋的重要工具.

空間-波數(shù)域算法，是一種結(jié)合傅里葉變換算法和空間域方程的方法，有基于積分方程和微分方程兩條思路，在重磁位場（李昆等，2019；Dai et al.，2019）、直流電正演（Dai et al.,2021）計(jì)算中已成功應(yīng)用，文獻(xiàn)（戴世坤等，2022）基于Coulomb規(guī)范的微分方程實(shí)現(xiàn)了空間-波數(shù)域三維電磁場數(shù)值模擬，體現(xiàn)了這種方法的低內(nèi)存占用、高效高精度的特性.本文基于Lorenz規(guī)范實(shí)現(xiàn)了空間-波數(shù)域三維大地電磁場數(shù)值模擬.用二次場方法，利用沿水平方向的二維傅里葉變換，將Lorenz矢量位滿足的三維控制方程轉(zhuǎn)換為多個(gè)波數(shù)之間相互獨(dú)立的常微分方程組，用有限單元法求解一維常微分方程.引入壓縮算子，構(gòu)成了穩(wěn)定收斂的迭代格式.設(shè)計(jì)模型驗(yàn)證了算法的正確性、收斂性和適應(yīng)性，試驗(yàn)結(jié)果表明，新方法具有計(jì)算精度高、占用內(nèi)存少和計(jì)算效率高的特性，相比基于Coulomb規(guī)范的空間域數(shù)值模擬方法（戴世坤等，2022），計(jì)算效率進(jìn)一步提高.

1 理論方法

設(shè)時(shí)諧因子為eiω t，無源頻率域Maxwell方程組可寫為

（1）

（2）

（3）

（4）

（5）

將式（5）代入式（1）可得

（6）

由于標(biāo)量的梯度的旋度為0，引入標(biāo)量位Φ，可以將E重新寫為

（7）

（8）

式（8）即為矢量位和標(biāo)量位滿足的微分方程（Weiss，2013；陳輝等，2016）.

基于二次場方法進(jìn)行三維大地電磁數(shù)值模擬，將場分為背景場和二次場，背景介質(zhì)設(shè)為均勻?qū)訝罱橘|(zhì)，此時(shí)，背景場和總場滿足的方程可寫為

（9）

（10）

（11）

（12）

式（12）即為二次場矢量位和標(biāo)量位滿足的微分方程.

1.1 基于Lorenz規(guī)范的空間域矢量位控制方程

引入二次場矢量位與二次場標(biāo)量位之間的Lorenz規(guī)范，表達(dá)式為

（13）

式（12）可改寫為

（14）

（15）

（16）

（17）

利用電磁場的邊界條件可得三個(gè)矢量位的邊界條件，將式（16）滿足的偏微分方程加載相應(yīng)的邊界條件，即可對(duì)三維大地電磁場進(jìn)行求解.

采用傳統(tǒng)空間域數(shù)值模擬算法求解式（16）的偏微分方程時(shí)，將構(gòu)成一個(gè)大型稀疏線性方程組，計(jì)算量大，存儲(chǔ)要求高，且規(guī)模越大，計(jì)算效率越低.空間-波數(shù)域方法對(duì)偏微分方程進(jìn)行水平方向二維傅里葉變換，將一個(gè)三維問題分解成不同波數(shù)下相對(duì)獨(dú)立的一維常微分方程問題，一維方程簡單且計(jì)算量小，能大大減少計(jì)算量和內(nèi)存需求，提高三維數(shù)值模擬的效率.

1.2 空間-波數(shù)域矢量位控制方程

（18）

上邊界

（19）

下邊界

（20）

Lorenz規(guī)范下，利用電磁場水平分量在邊界連續(xù)，垂直方向電流密度連續(xù)的特點(diǎn)，當(dāng)背景為突變介質(zhì)時(shí)，突變位置的邊界條件如下：

（21）

綜合式（18）—（21），即為基于Lorenz規(guī)范下的空間-波數(shù)域三維電磁場矢量位的邊值問題.

1.3 有限單元法

采用基于二次插值的一維有限單元法對(duì)式（18）—（21）組成的邊值問題進(jìn)行求解，利用伽遼金方法可轉(zhuǎn)化為有限元方程：

（22）

（23）

1.4 空間-波數(shù)域電磁場計(jì)算

利用有限單元法求得空間-波數(shù)域矢量位后，利用矢量位與電磁場之間的關(guān)系式（17）在空間-波數(shù)域中的表達(dá)式可求得空間-波數(shù)域二次電磁場：

（24）

（25）

式（24）和（25）中z方向的求導(dǎo)均采用差分求導(dǎo)，將空間-波數(shù)域電磁場進(jìn)行水平方向二維反傅里葉變換，即可求得到空間域電磁場.

1.5 迭代求解格式

將式（22）右端項(xiàng)中的總場采用背景場替代，得到的解為born近似解，在異常體電導(dǎo)率與背景場電導(dǎo)率差異很?。?倍以內(nèi)）時(shí)，能近似作為真解.當(dāng)電導(dǎo)率對(duì)比度差異大時(shí)，born近似解誤差大，本文采用迭代法逐次逼近真解，其核心思想是構(gòu)造一個(gè)穩(wěn)定收斂的迭代算子.在基于積分方程法的三維電磁場數(shù)值模擬中，Singer（1995）、Pankratov等（1997），Zhdanov和Fang（1997）和Avdeev等（2002）構(gòu)建一種收斂的波恩級(jí)數(shù)，把電磁場積分方程中格林算子G，利用坡印廷定理和能量不等式（Hursán and Zhdanov，2002）構(gòu)造了格林算子的線性變換，獲得了一個(gè)范數(shù)小于1的修正格林算子‖Gm‖<1.構(gòu)造了滿足任意有損耗介質(zhì)的電磁場迭代計(jì)算穩(wěn)定收斂的格式（Gao and Torres-Verdin,2006）：

E（n）=αE（n）+βE（n-1），

（26）

式（26）中，左端項(xiàng)E（n）是通過壓縮算子更新的第n次正演總場值，右端E（n）和E（n-1）分別表示第n和n-1次正演計(jì)算得到的總場，n=1時(shí)，E0=Eb，式中α，β的表達(dá)式與背景電導(dǎo)率σb、異常體電導(dǎo)率與背景電導(dǎo)率的差Δσ有關(guān)，寫為

（27）

（28）

若相鄰兩次迭代的計(jì)算結(jié)果滿足給定的精度要求，則停止迭代，輸出最新一次的正演總場值；否則將利用式（26）更新的第n次正演總場值作為第n+1迭代的初始總場，代入式（22）中右端項(xiàng)進(jìn)行新一輪的正演求解，直到滿足精度要求為止.

基于積分法與微分法的統(tǒng)一性，本文沿用積分方程法中構(gòu)造的迭代格式（27）—（28），測試結(jié)果證明，對(duì)于任意有損耗的介質(zhì)，算法的迭代過程穩(wěn)定收斂.

綜上，本文提出的基于Lorenz規(guī)范的空間-波數(shù)域三維電磁場數(shù)值模擬方法的要點(diǎn)主要有以下四點(diǎn)：

（1） 從Maxwell方程組出發(fā)，引入Lorenz規(guī)范，將電磁場滿足的偏微分方程轉(zhuǎn)化為關(guān)于Lorenz矢量位的亥姆霍茲方程組，基于二次場的數(shù)值模擬方法，得到關(guān)于二次場Lorenz矢量位的亥姆霍茲方程組；

（2） 利用水平方向二維傅里葉變換，將二次場矢量位滿足的三維空間域偏微分方程轉(zhuǎn)化為多個(gè)一維空間-波數(shù)域常微分方程組，一維方程組計(jì)算量小，求解簡單，大大減少了計(jì)算量和存儲(chǔ)需求；

（3） 常微分方程組采用有限單元法求解，單元采用二次插值形函數(shù)，構(gòu)成三個(gè)五對(duì)角方程，采用追趕法求解，且不同波數(shù)之間的方程組求解相互獨(dú)立，并行性好，垂向網(wǎng)格剖分靈活，進(jìn)一步提升了算法的計(jì)算精度和計(jì)算效率；

（4） 采用迭代法逐次逼近進(jìn)行求解，借用積分方程法中構(gòu)造的壓縮算子，算法穩(wěn)定收斂，占用內(nèi)存小，結(jié)合空間-波數(shù)域三維電磁場數(shù)值模擬算法，進(jìn)一步提高了算法的計(jì)算效率.

2 算法與算例

圖1所示為本文算法的流程圖.圖中ε〈〉為迭代終止的誤差計(jì)算公式，εmin為設(shè)定的最大迭代誤差，本文設(shè)置的迭代終止條件為：相鄰兩次迭代所有節(jié)點(diǎn)電場模的總和的相對(duì)誤差小于εmin=10-4，表達(dá)式為

（29）

式中|En|表示第n次迭代的總場的模，|En+1|表示第n+1次迭代的總場的模.

傅里葉變換的精度受波數(shù)選取、網(wǎng)格剖分的影響，本文模型水平方向網(wǎng)格均勻剖分，波數(shù)選取規(guī)律滿足采樣定理.測試表明，綜合模型復(fù)雜程度和頻率趨膚深度的影響選擇合適的網(wǎng)格剖分間距，并利用采樣定理選取的離散波數(shù)進(jìn)行數(shù)值模擬，數(shù)值解基本能滿足精度要求.具體的波數(shù)選取規(guī)則詳見文獻(xiàn)（陳龍偉等，2016），本文不再贅述.

本節(jié)設(shè)計(jì)模型進(jìn)行大地電磁場數(shù)值模擬，驗(yàn)證了算法的正確性，研究了算法的收斂性，并通過對(duì)比其他算法分析了算法效率.本文模型采用擴(kuò)邊FFT進(jìn)行計(jì)算，介質(zhì)磁導(dǎo)率μ均設(shè)為真空磁導(dǎo)率μ0，介電常數(shù)ε均設(shè)為真空介電常數(shù)ε0.本文測試的計(jì)算機(jī)為Intel（R） Core（TM） i9-7980XE CPU 主頻為2.60 GHz，內(nèi)存為64 GB.

2.1 正確性驗(yàn)證

用猶他大學(xué)開發(fā)的基于積分方程法（Integral Equation algorithm,IE）的三維正演軟件INTEM3D的計(jì)算結(jié)果為參照，驗(yàn)證本文算法（Space-wavenumber domain finite element method,SWFEM）在不同頻率和不同異常體電導(dǎo)率情況下視電阻率和相位的計(jì)算精度.模型為半空間內(nèi)的長方體，如圖2所示，空氣、下半空間介質(zhì)的電導(dǎo)率分別是10-12S·m-1和0.01 S·m-1，地面z=0，模型計(jì)算范圍：x方向-500～500 m，y方向-500～500 m，z方向0～500 m.異常體范圍：x方向-100～100 m，y方向-100～100 m，z方向50～250 m.剖分網(wǎng)格節(jié)點(diǎn)個(gè)數(shù)101×101×101，三個(gè)方向均勻剖分，Δx、Δy均為10 m，Δz為5 m.觀測面z=0，測點(diǎn)數(shù)據(jù)101×101個(gè).

表1 不同頻率積分方程法與本文算法的視電阻率和相位均方根誤差Table 1 The RRMS of apparent resistivity and phase between SWFEM algorithm and IE algorithm at different frequencies

圖1 基于Lorenz規(guī)范空間-波數(shù)域大地電磁三維正演流程圖Fig.1 Flow chart of the 3D MT modeling based on Lorenz gauge in the space-wavenumber domain

圖2 長方體模型Fig.2 The sketch of prism model

圖3 不同頻率積分方程法與本文算法視電阻率和相位對(duì)比曲線Fig.3 Apparent resistivity and phase of SWFEM algorithm and IE algorithm at different frequencies

圖4 不同頻率積分方程法與本文算法視電阻率和相位的相對(duì)誤差曲線Fig.4 Relative errors of apparent resistivity and phase between SWFEM algorithm and IE algorithm at different frequencies

測試分為兩組，第一組測試給定異常體電導(dǎo)率為0.1 S·m-1，頻率分別設(shè)為0.01 Hz、1 Hz、100 Hz和10000 Hz，地面y=0測線上視電阻率和相位的計(jì)算結(jié)果及其誤差如表1、圖3和圖4所示.由圖表可知，不同頻率本文算法與積分方程算法視電阻率和相位數(shù)據(jù)曲線擬合程度好，相對(duì)誤差均小于1%，均方根誤差均小于0.1%，計(jì)算精度高，驗(yàn)證了算法的正確性，且表明算法在不同頻率下均能計(jì)算正確.

第二組測試給定頻率為10 Hz，異常體電導(dǎo)率分別設(shè)為0.0001 S·m-1、0.001 S·m-1、0.1 S·m-1和1 S·m-1，地面y=0測線上視電阻率和相位的計(jì)算結(jié)果及其誤差如表2、圖5和圖6所示.綜合圖表可知，不同電導(dǎo)率異常體模型本文算法與積分方程算法結(jié)果吻合，當(dāng)異常體電導(dǎo)率為1 S·m-1時(shí)，與背景電導(dǎo)率對(duì)比度為100，測線上視電阻率最大相對(duì)誤差小于2%，相位小于0.3%；其他異常體模型兩種數(shù)值解的相對(duì)誤差均小于1%.兩種數(shù)值解地面視電阻率和相位的均方根誤差均小于0.15%，表明算法對(duì)于不同電導(dǎo)率對(duì)比度異常體模型的適應(yīng)性好.

2.2 收斂性分析

本節(jié)研究算法的收斂性，測試了不同電導(dǎo)率異常體迭代計(jì)算的收斂速度，記錄了電場三個(gè)分量的迭代誤差隨迭代次數(shù)的變化.模型為半空間內(nèi)有一正方體異常體，如圖7所示，空氣、下半空間電導(dǎo)率分別為10-12S·m-1和0.01 S·m-1，計(jì)算頻率為10 Hz.模型計(jì)算范圍：x方向-1000～1000 m，y方向-1000～1000 m，z方向0～1000 m，地面z=0.剖分網(wǎng)格節(jié)點(diǎn)個(gè)數(shù)101×101×101，三個(gè)方向均勻剖分，Δx、Δy均為20 m，Δz為10 m.異常體大小400×400×400 m3，范圍：x方向-200～200 m，y方向-200～200 m，z方向200～600 m.將異常體電導(dǎo)率分別設(shè)為高阻和低阻，與背景電導(dǎo)率對(duì)比度均分別為5、10、20、50和100倍，高阻異常體電導(dǎo)率分別為0.0001 S·m-1、0.0002 S·m-1、0.0005 S·m-1、0.001 S·m-1和0.005 S·m-1；低阻異常體電導(dǎo)率分別為0.05 S·m-1、0.1 S·m-1、0.2 S·m-1、0.5 S·m-1和1 S·m-1.進(jìn)行大地電磁正演，以x極化模式下電場迭代收斂情況為例，研究電導(dǎo)率不同對(duì)比度情況下電場三分量的收斂速度.

圖6 不同異常體電導(dǎo)率積分方程法與本文算法計(jì)算的視電阻率和相位的相對(duì)誤差曲線Fig.6 Relative errors of apparent resistivity and phase between SWFEM algorithm and IE algorithm at different conductivity of anomalies.

表2 不同電導(dǎo)率異常體積分方程法與本文算法視電阻率和相位均方根誤差Table 2 The RRMS of apparent resistivity and phase between SWFEM algorithm and IE algorithm at different conductivity of anomalies

圖8和圖9分別為不同高阻異常體和低阻異常體電場三分量的迭代誤差變化，圖中黑色橫線為本文設(shè)置的迭代終止誤差，εmin=10-4.由圖中可知，隨著迭代的增加，電場三分量的誤差均逐漸減小，且與背景電導(dǎo)率差異越大，達(dá)到計(jì)算精度所需迭代次數(shù)越多；相比高阻，低阻電場Ex和Ey分量達(dá)到計(jì)算精度所需的迭代次數(shù)更多，Ez分量在電導(dǎo)率對(duì)比度更大時(shí)迭代收斂曲線的震蕩更劇烈，總體來說，低阻異常體比高阻異常體的收斂慢，但算法均能穩(wěn)定收斂，表明壓縮算子適用于本文算法.

2.3 Dublin模型測試

采用Dublin（DTM1）模型測試方法的效率.均勻半空間中有三個(gè)長方體異常，模型計(jì)算范圍x方向-37.5～37.5 km；y方向-25～25 km；z方向0～55 km，三個(gè)方向均采用均勻剖分.半空間地下電導(dǎo)率為σ0=0.01 S·m-1，空氣電導(dǎo)率為10-12S·m-1.模型如圖10所示，異常體的范圍和電導(dǎo)率參數(shù)如表3所示.計(jì)算頻率為0.01 Hz.

文獻(xiàn)（戴世坤等，2022）基于Coulomb規(guī)范實(shí)現(xiàn)了空間-波數(shù)域三維電磁場數(shù)值模擬，本文采用Lorenz規(guī)范空間-波數(shù)域方法進(jìn)行三維電磁場數(shù)值模擬.二者采用規(guī)范不同導(dǎo)致最后得到的空間-波數(shù)域控制方程不同，Coulomb規(guī)范有4個(gè)未知量（三個(gè)矢量位、一個(gè)標(biāo)量位），各未知量滿足的方程相互耦合；Lorenz規(guī)范有3個(gè)未知量（三個(gè)矢量位），在各向同性介質(zhì)中，未知量方程能相對(duì)獨(dú)立.將DTM1 模型采用不同網(wǎng)格剖分，對(duì)比Coulomb規(guī)范和Lorenz規(guī)范空間-波數(shù)域算法三維大地電磁正演效率，表4為統(tǒng)計(jì)結(jié)果，記錄了兩種算法不同網(wǎng)格剖分情況下單次迭代時(shí)間、正演總時(shí)間和內(nèi)存占用，圖11和12為相應(yīng)的變化曲線.

圖7 正方體模型Fig.7 The sketch of cube model

圖8 不同高阻異常體的電場三分量迭代誤差曲線Fig.8 Iterative error curves of three electric fields for high resistance anomaly

圖9 不同低阻異常體的電場三分量迭代誤差曲線Fig.9 Iterativeerror curves of three electric fields for low resistance anomaly

圖10 DTM1 模型示意圖Fig.10 Schematic diagram of DTM1 model

表3 DTM1 模型三個(gè)長方體異常的分布范圍和電導(dǎo)率Table 3 Range and conductivity of three cuboid anomalies in DTM1 model

綜合圖表可知，兩種空間-波數(shù)域算法耗時(shí)和內(nèi)存占用隨計(jì)算節(jié)點(diǎn)數(shù)增多均近似線性增長，且相同網(wǎng)格剖分情況下，本文基于Lorenz規(guī)范下的空間-波數(shù)域正演方法比Coulomb規(guī)范的耗時(shí)更短、內(nèi)存需求更小.對(duì)于872289個(gè)節(jié)點(diǎn)數(shù)模型，本文算法單次迭代時(shí)間約0.47 s，基于Coulomb規(guī)范算法約1.91 s；本文算法總耗時(shí)約240.09 s，基于Coulomb規(guī)范算法約962.43 s；本文算法內(nèi)存占用約458.5 MB，基于Coulomb規(guī)范算法588.3 MB，兩種算法內(nèi)存需求差別小，但本文算法耗時(shí)減少近3倍.

DTM1模型是一個(gè)復(fù)雜模型，其計(jì)算規(guī)模大、異常體范圍大，異常體電導(dǎo)率與背景電導(dǎo)率對(duì)比度差異大.而電導(dǎo)率差異大會(huì)導(dǎo)致本文算法迭代終止時(shí)迭代次數(shù)較多，影響效率.文獻(xiàn)（Long and Farquharson,2019）采用基于徑向基函數(shù)的無網(wǎng)格單元法進(jìn)行三維大地電磁數(shù)值模擬，計(jì)算DTM1模型范圍與本文相同，在16核Intel Xeon E5-2670 CPU、主頻為2.6 GHz的電腦上，計(jì)算120598個(gè)離散節(jié)點(diǎn)，正演一次耗時(shí)約797.7 s.在節(jié)點(diǎn)數(shù)為222589時(shí)，Coulomb規(guī)范的空間-波數(shù)域算法總耗時(shí)312.10 s，本文算法總耗時(shí)約69.82 s.節(jié)點(diǎn)數(shù)比文獻(xiàn)（Long and Farquharson,2019）多，在配置更差的電腦上，空間-波數(shù)域算法效率更高，且本文算法的優(yōu)勢(shì)更明顯.

圖13為節(jié)點(diǎn)數(shù)為61×41×89的模型本文算法正演結(jié)果與Siripunvaraporn等（2002）的對(duì)比圖，數(shù)據(jù)來源于文獻(xiàn)（Miensopust et al.,2013），測線y=0 m上視電阻率的相對(duì)誤差最大為3%，平均相對(duì)誤差為1.2%，相位的相對(duì)誤差均在1%以下，滿足精度要求，表明本文算法對(duì)復(fù)雜模型適應(yīng)性好，且算法效率高.

表4 Dublin模型不同網(wǎng)格剖分計(jì)算效率Table 4 Calculation efficiency of Dublin model at different grids

圖11 不同網(wǎng)格剖分Coulomb規(guī)范和Lorenz規(guī)范算法的計(jì)算時(shí)間Fig.11 Calculation time of Coulomb gauge algorithm and Lorenz gauge algorithm of space-wavenumber domain at different grids

圖12 不同網(wǎng)格剖分Coulomb規(guī)范和Lorenz規(guī)范算法的內(nèi)存占用Fig.12 Memory of Coulomb gauge algorithm and Lorenz gauge algorithm of space-wavenumber domain at different grids

3 結(jié)論

本文提出一種基于Lorenz規(guī)范的空間-波數(shù)域三維大地電磁數(shù)數(shù)值模擬方法.引入Lorenz規(guī)范，將Maxwell方程轉(zhuǎn)換成矢量位的偏微分方程組.采用二次場方法，將總場拆分為背景場和二次場，處理二次場時(shí)，利用水平方向二維傅里葉變換，將二次場矢量位滿足的偏微分方程轉(zhuǎn)換為多個(gè)波數(shù)的常微分方程，常微分方程采用追趕法求解，計(jì)算量小、存儲(chǔ)需求小、計(jì)算速度快.最后利用壓縮算子，采用迭代法求解空間域電磁場.取得以下結(jié)論：

（1）設(shè)計(jì)不同頻率和不同電導(dǎo)率的異常體驗(yàn)證了算法的正確性和適應(yīng)性；

（2）分析異常體電導(dǎo)率與背景電導(dǎo)率不同對(duì)比度時(shí)算法的迭代收斂性，對(duì)比度越大，收斂越慢，相同對(duì)比度的高阻異常體比低阻異常體的迭代收斂速度更快；

（3）采用DTM1 模型對(duì)比了基于Coulomb規(guī)范和Lorenz規(guī)范的空間-波數(shù)域三維大地電磁正演效率，隨著計(jì)算規(guī)模的增大，兩種算法耗時(shí)與存儲(chǔ)均呈近似線性增長，但Lorenz規(guī)范的算法耗時(shí)更短、內(nèi)存需求更少.相比傳統(tǒng)算法，空間-波數(shù)域算法效率高，且本文算法優(yōu)勢(shì)更明顯，非常適合大規(guī)模三維大地電磁正演.

圖13 地面y=0 m測線上視電阻率和相位Fig.13 The apparent resistivity and phase in line y=0 m on surface

本文僅研究大地電磁場三維數(shù)值模擬，若加上外加源項(xiàng)，即可進(jìn)行帶源的電磁場三維數(shù)值模擬，當(dāng)頻率為0時(shí)，可以退化為直流電模擬.為大規(guī)模三維電磁法數(shù)值模擬的快速計(jì)算提供了一種高效的新方法，為精細(xì)、高效反演成像提供新工具.下一步將研究帶地形的處理和各向異性介質(zhì)的數(shù)值模擬.

致謝感謝審稿專家和編輯提出的寶貴意見．

附錄A 單元積分形式

式（22）種存在三種形函數(shù)積分，查表（徐世浙,1994）可得三種積分的單元積分表達(dá)式.

（1）第一類積分

（A1）

（2）第二類積分

（A2）

（3）第三類積分

（A3）

其中，

以上三種類型的積分涵蓋了式（22）全部積分類型，其他積分項(xiàng)均可以在這三種積分類型中找到相應(yīng)的單元積分.