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

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

0 引言

大地電磁法利用天然電磁場進行地下電性結(jié)構(gòu)的勘探，由于其野外工作簡單、成本低、受野外地質(zhì)條件約束小和探測深度廣等優(yōu)點被廣泛應(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)鍵，也是準確進行實測數(shù)據(jù)資料處理和地質(zhì)地球物理解釋的重要工具.

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

1 理論方法

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

（1）

（2）

（3）

（4）

（5）

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

（6）

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

（7）

（8）

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

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

（9）

（10）

（11）

（12）

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

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

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

（13）

式（12）可改寫為

（14）

（15）

（16）

（17）

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

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

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

（18）

上邊界

（19）

下邊界

（20）

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

（21）

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

1.3 有限單元法

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

（22）

（23）

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

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

（24）

（25）

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

1.5 迭代求解格式

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

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

（26）

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

（27）

（28）

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

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

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

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

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

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

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

2 算法與算例

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

（29）

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

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

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

2.1 正確性驗證

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

表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 不同頻率積分方程法與本文算法視電阻率和相位對比曲線Fig.3 Apparent resistivity and phase of SWFEM algorithm and IE algorithm at different frequencies

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

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

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

2.2 收斂性分析

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

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

表2 不同電導率異常體積分方程法與本文算法視電阻率和相位均方根誤差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.由圖中可知，隨著迭代的增加，電場三分量的誤差均逐漸減小，且與背景電導率差異越大，達到計算精度所需迭代次數(shù)越多；相比高阻，低阻電場Ex和Ey分量達到計算精度所需的迭代次數(shù)更多，Ez分量在電導率對比度更大時迭代收斂曲線的震蕩更劇烈，總體來說，低阻異常體比高阻異常體的收斂慢，但算法均能穩(wěn)定收斂，表明壓縮算子適用于本文算法.

2.3 Dublin模型測試

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

文獻（戴世坤等，2022）基于Coulomb規(guī)范實現(xiàn)了空間-波數(shù)域三維電磁場數(shù)值模擬，本文采用Lorenz規(guī)范空間-波數(shù)域方法進行三維電磁場數(shù)值模擬.二者采用規(guī)范不同導致最后得到的空間-波數(shù)域控制方程不同，Coulomb規(guī)范有4個未知量（三個矢量位、一個標量位），各未知量滿足的方程相互耦合；Lorenz規(guī)范有3個未知量（三個矢量位），在各向同性介質(zhì)中，未知量方程能相對獨立.將DTM1 模型采用不同網(wǎng)格剖分，對比Coulomb規(guī)范和Lorenz規(guī)范空間-波數(shù)域算法三維大地電磁正演效率，表4為統(tǒng)計結(jié)果，記錄了兩種算法不同網(wǎng)格剖分情況下單次迭代時間、正演總時間和內(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 模型三個長方體異常的分布范圍和電導率Table 3 Range and conductivity of three cuboid anomalies in DTM1 model

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

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

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

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

圖11 不同網(wǎng)格剖分Coulomb規(guī)范和Lorenz規(guī)范算法的計算時間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)換成矢量位的偏微分方程組.采用二次場方法，將總場拆分為背景場和二次場，處理二次場時，利用水平方向二維傅里葉變換，將二次場矢量位滿足的偏微分方程轉(zhuǎn)換為多個波數(shù)的常微分方程，常微分方程采用追趕法求解，計算量小、存儲需求小、計算速度快.最后利用壓縮算子，采用迭代法求解空間域電磁場.取得以下結(jié)論：

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

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

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

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

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

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

附錄A 單元積分形式

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

（1）第一類積分

（A1）

（2）第二類積分

（A2）

（3）第三類積分

（A3）

其中，

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