基于吸收邊界條件的瞬變電磁法三維矢量有限元快速正演

張永超， 王光杰， 李宏杰， 廉玉廣， 李文，邱浩， 牟義

1 煤炭科學技術研究院有限公司安全分院， 北京 100013 2 煤炭資源高效開采與潔凈利用國家重點實驗室(煤炭科學研究總院)， 北京 100013 3 北京市煤礦安全工程技術研究中心, 北京 100013 4 中國科學院大學， 北京 100049 5 中國科學院地質與地球物理研究所 頁巖氣與地質工程重點實驗室， 北京 100029

0 引言

瞬變電磁法又稱時間域電磁法，是近年來發展迅速的一種地球物理勘探方法.它具有可在近區測量、旁側效應影響小、對低阻體敏感等優點(樸化榮，1990)，被廣泛地應用于水文地質調查、金屬礦勘探、工程勘查等領域.目前，實際工作中瞬變電磁法的反演解釋仍以一維方法為主(薛國強等，2008)，分辨率較低，在某些地質條件復雜的地區難以取得理想效果.Di等(2021)耦合彈性波和電磁波數據，解譯了華南陸塊深部-淺部的地質結構.為更好地指導瞬變電磁法的資料處理、反演解釋、進一步改善勘探效果，開發速度快、精度高、對復雜地電模型適用性好的三維正演技術已成必由之路.

目前常用的三維數值模擬方法主要有：積分方程法(Sanfilipo and Hohmann，1985；Zhdanov et al.，2006)、有限差分法(Wang and Hohmann，1993；Commer and Newman，2004；孫懷鳳，2018,2021)和有限元法(殷長春等，1994, 2013；李瑞雪等，2016a, 2016b；余翔等，2017；周建美等，2018；Zeng et al., 2019)，其中有限元法具有對復雜地質模型的適用性強、形成的總體系數矩陣是稀疏且對稱、易于并行等優點，成為電磁法領域解決正演問題應用最為廣泛的方法(Jin，2014).Kuo和Cho(1980)利用有限元法對低阻覆蓋層下的低阻礦體的瞬變電磁法響應進行了三維數值模擬，得到的曲線形狀和數量級與野外實測數據基本一致，證明了該方法的可行性.此后，Pridmore等(1981)也對有限元求解三維瞬變電磁法正演問題進行了探索性研究.Sugeng(1998)采用異常場和背景場分離的思路，使用六面體矢量單元對回線源瞬變電磁法進行了三維正演.Jin等(1999)基于蘭喬斯譜分解法(SLDM)，分別運用節點有限元法和矢量有限元法對時間域和頻率域的三維電磁場進行了數值模擬.Ralph-Uwe等(2008)提出了一種快速的基于頻域有限元法離散、Krylov子空間投影技術建立的模型簡化方法求得頻域響應，然后通過余弦變換求得時間域電磁場響應.Um等(2010)提出了用于電性源瞬變電磁法的時域矢量有限元三維正演算法，其模型剖分采用四面體網格，時間離散使用隱式后退歐拉法，邊界條件選用Dirichlet齊次邊界條件，模擬了海洋時間域電磁場.李建慧(2011，2013)，劉亞軍(2019)利用矢量有限元法對大回線源瞬變電磁法進行了三維正演，其算法從電場雙旋度方程出發，模型剖分采用六面體單元，邊界條件為電磁場切向分量為零的齊次邊界條件，先通過異常場/背景場法求得頻域響應，然后通過Gaver-Stehfest變換得到時間域響應.張博等(2016)采用矢量有限元法對帶有起伏地形的三維航空瞬變電磁法進行了數值模擬，其思路與李建慧的基本相同，不同之處在于為了更好地模擬地形起伏采用了非結構化的Delaunay四面體進行網格剖分，時頻變換則采用余弦變換進行.李賀(2016)也基于時間域電場雙旋度方程和電磁場切向分量為零的齊次邊界條件，對典型地電模型的瞬變電磁法響應進了三維正演，模型剖分采用六面體矢量單元，時間離散則采用向后差分方式.李建慧(2017)基于齊次邊界條件，采用將回線源視為多個電偶源的策略和非結構化四面體矢量單元正演了復雜形態回線源激發的瞬變電磁場分布.

綜上所述，針對瞬變電磁法的三維有限元正演，目前常采用齊次邊界條件作為模型的截斷邊界條件，這雖然簡化了有限元弱形式方程的推導和后續的單元分析，但為了滿足該條件，需要構建較大尺寸的模型，否則會產生由邊界引起的非物理反射，降低正演結果的精度.而較大尺寸的模型意味著更多的網格單元，這會降低正演問題的求解速度.吸收邊界條件能使邊界的非物理反射最小化，因此可以縮小模型的尺寸、減少單元數目、加快求速解速度，或者在模型大小相同時提高正演結果的精度(Jin，2014).鑒于此，本文在前人研究的基礎上引入吸收邊界條件，以加快瞬變電磁法三維有限元正演的求解速度.

1 基于吸收邊界條件的時間域矢量有限元法

1.1 矢量勢控制方程推導

考慮場源的情況下，時間域中微分形式的麥克斯韋方程組如下

(1)

(2)

(3)

(4)

此外，還有三個輔助的介質特性方程如下

B=μH,D=εE,J=σE,

(5)

式中μ、ε、σ分別為磁導率、介電常數和電導率.針對地學問題，常見礦物的磁導率與真空磁導率十分接近，故μ=μ0.

由(2)、(4)式可假設

(6)

(7)

式中，A稱為矢量勢，φ稱為標量勢.

根據(3)、(5)和(7)式可得

(8)

將(6)、(8)式和介質特性方程代入(1)式，可得

(9)

1.2 吸收邊界條件的選擇

與頻率域電磁法相比，瞬變電磁法的吸收邊界條件較為困難.針對這一難題，本文的解決思路為：對于瞬變電磁法，由模型邊界引起的非物理反射主要影響晚期數據，而晚期數據的頻段集中在低頻，針對這一特性，借鑒頻率域吸收邊界條件的思想，構建對低頻吸收效果好的時間域吸收邊界條件.

由于模型規模遠大于發射線圈尺寸，外表面處的電磁波可以近似為以發射線圈為中心的球面波.借鑒Webb和Kanellopoulos(1989)的研究成果，采用形如下式的對球面波吸收較好的一階吸收邊界條件來減小模型規模，提升正演速度.

(10)

式中n為模型外表面的單位法向量，r為發射線圈中心點到外表面上任意一點的距離，ε′=(ε+|ε*|)/2，ε*=ε-iσ/2πf，為復介電常數，與介電常數、電導率和指定的主吸收頻率f相關，f可以在n～n×10 Hz中選取.

1.3 弱形式推導

根據Galerkin加權余量法，結合矢量恒等式和散度定律，基于(9)式控制方程和(10)式邊界條件的電磁場問題的弱形式方程為

(11)

式中，N為矢量基函數，Ω為模型區域，Γ為其外表面.

1.4 單元分析

矢量單元(又稱棱邊單元)很好地解決了節點單元由于未強加散度條件或強加法向連續性邊界條件造成的“偽解”問題以及與結構相關的(如導體或介質的邊緣或銳角)奇異性問題(Jin，2014)，因此本文采用非結構化的一階四面體矢量單元進行單元分析，相比六面體單元不僅能更好地模擬地形起伏或傾斜界面等復雜地電模型，還能在場源和異常體附近細分網格，提高結果精度.一階四面體矢量單元的節點與棱邊間的編碼規則見圖1.

圖1 四面體矢量單元的節點和棱邊示意圖Fig.1 Sketch map of nodes and edges in a tetrahedral vector unit

單元e中，點(x,y,z)處的矢量勢A在某一時刻t的值可由六條棱邊的矢量基函數表示為(Jin，2014)

(12)

式中，上標e表示棱邊i等編碼均為基于單元的局部編碼.其中矢量基函數

i=1,2,3,4,5,6

(13)

矢量基函數具有如下性質：

①散度為零，即

(14)

將(13)式代入(11)式，即可對弱形式方程進行空間離散.首先對(11)式左端體積分的第一部分進行離散，可得

i,j=1,2,3,4,5,6

(15)

若單元e位于外表面Γ上則需要進一步計算(11)式左端的面積分部分：

(16)

對于(11)式的右端的計算，由于線電流的奇異性需將回線源看做多個較短導線元的組合，并在源附近采用較密的網格剖分.Ansari和Farquharson(2014)、李建慧等(2016)、李建慧等(2017)以該方法處理接地長導線源，取得了很好的效果.下面以圖2中的Ⅰ邊為例，詳述發射電流加載過程.

圖2 發射回線示意圖Fig.2 Sketch map of transmitter loop

設發射源位于z=0的平面上，四邊中通過的電流方向如圖2所示.網格剖分后，設Ⅰ邊分布于若干個單元中，且總是與這些單元的某一個棱邊重合，其源電流密度為

(17)

式中，I(t)為隨時間變化的電流，指定不同的I(t)即可實現不同發射電流波形的加載，δ為沖激函數(又稱狄拉克函數)，dl為導線元的長度，(0,0，0)為

其中心坐標，ex表示x方向的單位矢量.

令{Qe}為6×1的列向量，其元素為

(18)

(19)

令[Re]=[Ee]+[Ke]、[Pe]=[Fe]+[Me]，將上述單元矩陣后組裝起來后，最終得到形如下式的總體矩陣方程：

(20)

1.5 時間離散

(20)式中未知數A(t)的求解可采用時間逐步積分法，它包括以中心差分法為代表的時間步長受算法穩定性限制的顯式算法和以Newmark法為代表的時間步長只受精度限制的隱式算法(王長清和祝西里，2011).顯示算法不需要計算逆矩陣，但為滿足穩定性條件，時間步長必須很小，選擇顯式算法有極大的計算量.根據瞬變電磁場前期變化劇烈、后期變化則相對平緩的特點，本文選擇時間步長可逐漸增大的Newmark法進行求解.

該方法原理如下，設t時刻的A(t)已經求得，則t+Δt時刻

(21)

(22)

(23)

針對系數矩陣為對稱稀疏矩陣的特點，對(23)式的求解，采用針對稀疏矩陣的Intel MKL PARDISO求解器直接求解.

2 算法驗證

2.1 均勻半空間

均勻半空間模型的大小為6 km×6 km×6 km，其中空氣層厚度2 km，電導率10-12S·m-1，巖石厚度4 km，電導率0.01 S·m-1，發射線框為邊長100 m的正方形，水平置于巖石和空氣交界面(z=0)的中心，接收線圈面積100 m2，發射電流20 A，電流下降緣波形設為指數下降(如圖3所示)，關斷時間10-5s.

圖3 電流關斷波形Fig.3 Waveform of the turn-off current

模型的網格剖分采用約束Delaunay算法，共剖分3764個單元，起始時間步長設為10-10s，每5次將步長放大100.1倍，共計算336步，在CPU為i5-8250U，內存8GB的筆記本電腦上，計算時間為164 s.計算完成后通過插值提取了按對數等間距分布的10-5s至0.01 s共31個時間門的數據，并將其與解析解進行對比(如圖4所示)，除一個時間門的相對誤差稍大于3%外其余時間門的相對誤差的絕對值均小于3%，31門數據的均方相對誤差為1.37%，兩者吻合良好，證明了本文算法正確.

圖4 均勻半空間數值解及其相對誤差(σ=0.01 S·m-1, 6 km×6 km×6 km)Fig.4 Numerical solution and relative errors of homogeneous half-space (σ=0.01 S·m-1, 6 km×6 km×6 km)

其他設置不變的情況下，將吸收邊界條件替換為齊次邊界條件(即強制式(11)中的面積分部分為零)，此時的數值解與吸收邊界條件相比，早期的精度基本一致，但隨著時間的增加，電磁場逐漸向邊界擴散，誤差也逐漸增大，最大相對誤差接近-10%，在晚期的精度顯著低于吸收邊界條件.

圖5為均勻半空間模型尺寸為10 km×10 km×10 km(其中空氣層厚3 km)時，采用齊次邊界條件得到的數值解與解析解的對比，計算時間為231 s.31門數據的均方相對誤差為1.38%，與采用吸收邊界條件、模型大小為6 km×6 km×6 km時的精度相近.

圖5 均勻半空間數值解及其相對誤差(σ=0.01 S·m-1, 10 km×10 km×10 km)Fig.5 Numerical solution and relative errors of homogeneous half-space (σ=0.01 S·m-1, 10 km×10 km×10 km)

通過以上對比可見，采用吸收邊界條件確實可以縮小模型尺寸、加快正演速度，提高正演的精度.

2.2 H型地電斷面

將H型地電斷面模型、中心回線裝置的數值解與CR1Dmod一維正演程序的結果進行對比，進一步驗證本文算法的正確性.模型圍巖的電導率為0.01 S·m-1，其中有一頂部埋深80 m、厚20 m、電導率0.1 S·m-1的低阻層，其余設置與均勻半空間模型相同.網格剖分時在發射線框的正下方、低阻層的頂面設置一個正方形以控制該層的網格剖分，改善求解精度(任政勇和湯井田，2009).由于低阻層及其附近需要剖分較小尺寸的網格，因此單元數較均勻半空間大幅增長，共剖分16528個單元.在10-4s之前，時間步長的設置與均勻半空間模型相同，10-4s之后，可以認為受低阻層的影響，電磁場的擴散速度變慢，因此每5次將步長放大100.2倍，共計算294步，計算時間為3720 s.

由圖6可見，有限元解和CR1Dmod解在早期偏離較大，在0.01 ms時相對誤差接近25%，這是因為受發射電流“暫態效應”的影響，早期場值與關斷時間為零的理想條件(CR1Dmod即采用此種假設)相比會有明顯的增大(白登海和Maxwell，2001).此后，隨著時間的增加，暫態效應的影響降低，有限元解和CR1Dmod解的相對誤差也快速減小，在0.1 ms后，相對誤差絕對值小于2.5%，進一步證明了本文算法的正確性.

圖6 H型地電斷面的有限元解和CR1Dmod解及其相對誤差Fig.6 FEM solution and CR1Dmod solution of H type geoelectric section and their relative errors

3 算例分析

3.1 水平低阻異常體

如圖7所示，在電導率為0.01 S·m-1的圍巖中有一規模為600 m×300 m×20 m、頂部埋深80 m、電導率0.1 S·m-1的低阻異常體，異常體中心在地表的投影為坐標原點，其余設置與均勻半空間模型相同.正演時，在異常體的地表投影基礎上外擴200 m范圍內(即x：-350至350，y：-500至500)，以50 m×50 m的網度、移動中心回線的方式計算不同位置的瞬變電磁響應.根據回線位置的不同，模型剖分的單元數4635～5749個不等，時間步長的設置與均勻半空間模型相同，計算時間為961～1872 s.

圖7 水平低阻異常體模型示意圖Fig.7 Sketchmap of horizontal low resistivity anomalous body model

圖8為水平低阻異常體與前述的H型地電斷面、均勻半空間模型在坐標原點O的瞬變電磁響應對比.前兩者在0.4 ms之前的響應基本相同，與均勻半空間相比，在0.03～0.06 ms都出現了由于反射波干涉引起的電壓低于均勻半空間的Undershoot現象，之后感應電壓(V)逐漸高于均勻半空間.0.4 ms之后，水平低阻異常體的響應開始快速下降并最終趨于均勻半空間，而H型地電斷面受低阻屏蔽效應的影響，感應電壓始終高于均勻半空間.

圖8 O點不同模型的瞬變電磁響應Fig.8 TEM response of different models at point O

圖9為具有代表性時間門的瞬變電磁響應切片.電流完全關斷的時刻(0.01 ms)感應電壓隨位置有幅度不大的無規律波動，推測是由發射框移動引起的網格剖分的差異造成的.0.05 ms時刻，由于Undershoot現象出現了一片低電壓區域，其基本呈橢圓形，與異常體的地表投影相比范圍略有擴大，中心區域的電壓相比正常值下降約16%.0.4 ms時刻，受異常體影響出現的高電壓區域，范圍相比異常體的投影有所擴大，與正常值相比，測點在距異常體邊緣約50 m處感應電壓的增幅大于10%，在異常體邊緣電壓的增幅大于25%，進入異常體投影后電壓開始迅速增長并在達到某一值后趨于平穩，電壓增幅最大約為160%.10 ms時刻，隨著波場的擴散，異常體的影響減弱，電壓趨于均勻.

圖9 水平低阻異常體的瞬變電磁響應時間切片Fig.9 Time slice of TEM response of horizontal low resistivity anomalous body model

圖10 XZ平面0.05～0.06 ms電流密度變化(a) 均勻半空間； (b) 水平低阻異常體.Fig.10 Change of current density in XZ plane from 0.05 ms to 0.06 ms(a) Homogeneous half-space; (b) Horizontal low resistivity anomalous body model.

下面嘗試對Undershoot現象做出解釋.根據正演結果，抽取了發射框中心位于O點時，均勻半空間和水平低阻異常體模型在0.05 ms、0.06 ms時刻XZ平面處的電流密度，將兩者相減得到前后電流密度變化ΔJ，結果見圖10，圖示以外的部分尚未被感應渦流的擴散波及或電流密度變化極小.與均勻半空間相比，受低阻異常體的影響ΔJ主要集中于異常體(紅色虛線所示)內，上部的ΔJ減小.以10 m×10 m的網格剖分圖11中的區域，假設在每個網格內ΔJ是均勻的，并將ΔJ形成的電流環視為磁偶極子，則其在O點引起的磁感應強度變化可按下式計算

(24)

式中ds為網格面積，(y,z)為網格中心坐標.均勻半空間的計算結果為-3.99×10-7T，而低阻異常體為-3.19×10-7T.上面的估算是十分粗略的，因為很多時候并不滿足源點和場點之間的距離遠大于電流環半徑這一磁偶極子的假設條件，但仍在一定程度上說明之所以出現Undershoot現象是因為在某些時段低阻異常體對ΔJ的集中不足以補償其對上部ΔJ減少的影響.

圖11 傾斜低阻異常體模型示意圖Fig.11 Sketch map of sloping low resistivity anomalous body model

3.2 傾斜低阻異常體

在前節模型的基礎上，將地質異常體傾角調整為30°，維持異常體頂面的中心O′至原點O的距離仍為80 m.正演時，以50 m×50 m的網度計算x：-250至250、y：-400至400矩形區域內的瞬變電磁響應.根據回線位置的不同，模型剖分的單元數5349～8792個不等，時間步長的設置與均勻半空間模型相同，計算時間為1117～2616 s.

從圖12可見，傾斜和水平低阻異常體在O點的瞬變電磁響應基本相同，只是由于傾斜低阻異常體的垂向等效厚度增大，在0.05～0.2 ms 之間感應電壓有7%～18%的提升.

圖12 O點水平和傾斜低阻異常體模型的瞬變電磁響應Fig.12 TEM response of models of horizontal and sloping low resistivity anomalous body at point O

y=0線的多測道剖面見圖13，圖中只顯示了部分具有代表性的時間門.0.05 ms時刻，感應電壓隨x坐標的增長呈先降后升的波浪式起伏；0.1～0.8 ms之間，感應電壓呈不對稱的單波峰形態，波峰偏向于異常體埋深較淺的一側并隨時間的增加向中心移動且振幅逐漸減小；1.6 ms 時，感應電壓在異常體投影范圍內有一定的升高但開始趨于均勻.

圖13 y=0線瞬變電磁多測道剖面Fig.13 TEM multichannel profile of y=0 line

圖14為具有代表性時間門的瞬變電磁響應切片.0.01 ms時刻，在異常體埋深較淺的一側出現了一片低電壓條帶.0.05 ms時刻，在異常體投影范圍內，埋深較淺的一側呈現高電壓，較深一側呈現低電壓，且高壓升幅大于低壓降幅.0.4 ms時刻，感應電壓的變化與水平低阻異常體的特征基本相同，但是電壓最大增幅更大(約180%)且波峰偏向于異常體埋深較淺一側.10 ms時刻，電壓的分布基本均勻.

4 結論

(1)從時間域麥克斯韋方程組出發，推導了基于矢量勢的微分控制方程，在此基礎上結合一階球面波吸收邊界條件推導了弱形式方程，采用一階四面體矢量單元進行單元分析、帶“數值阻尼”的Newmark法進行時間離散、并將回線源視為多個電偶源的策略，實現了任意發射電流波形的瞬變電磁法響應的時間域快速求解.

(2)通過有限元數值解與均勻半空間模型的解析解、H型地電斷面模型的CR1Dmod解的對比，證明了本文算法的正確性.均勻半空間模型采用吸收邊界條件和齊次邊界條件的結果對比則表明：吸收邊界條件確實可以提高瞬變電磁法三維正演的精度或者縮小模型尺寸、加快計算速度.

(3)根據正演結果，其他條件相同時，水平有限體積的低阻異常體和H型地電斷面的瞬變電磁響應前期基本相同，但在晚期前者會快速趨于均勻半空間的響應.對比均勻半空間和水平低阻異常體模型在0.05 ms、0.06 ms時刻電流密度的變化，表明感應電壓的Undershoot現象之所以出現，是因為低阻異常體對電流的集中不足以補償其引起的上部電流減少的影響.對于傾斜低阻異常體，由于垂向上等效厚度增大，中心處的感應電壓會高于相同條件的水平低阻異常體，在沿傾向的多測道剖面上感應電壓的波峰會向異常體埋深較淺的一側偏移.

致謝中國地質大學(武漢)胡祥云教授、中國科學院地質與地球物理研究所安志國副研究員對本文提出了有益的建議，在此表示感謝.此外，對審稿人提出寶貴的修改意見和編輯們的辛苦工作表示感謝.