湯文武 柳建新 葉益信 張 華(東華理工大學地球物理與測控技術學院，江西南昌 330013； 中南大學地球科學與信息物理學院，湖南長沙 410083)

1 引言

電磁法是根據電磁感應原理研究天然電磁場或人工(可控)場源在大地中激勵的交變電磁場分布， 并由觀測到的電磁場分布研究地下目標體的電性及地質特征的一種地球物理方法[1]，已被廣泛地應用于金屬礦勘探[2，3]、地下水資源調查[4]、油氣田勘探[5]、油氣動態(tài)監(jiān)測[6]等領域。根據場源的類型可分為天然源及可控源兩類，現今可控源電磁法發(fā)展迅速，但場源的引入給其資料處理和解釋帶來了諸多困難。作為資料解釋的重要手段，正反演技術的發(fā)展對進一步推動可控源電磁法的廣泛應用具有重要作用，尤其是三維正、反演方法，已成為研究熱點。

針對三維可控源電磁法的正演問題，近年來國內外學者開展了大量研究。Newman等[7]采用基于二次電場的Helmholtz 方程，將全空間或層狀半空間電磁場響應作為一次場引入到方程右端，避免了直接模擬電偶源或磁偶源的困難。利用交錯網格有限差分法對Helmholtz方程進行離散，得到的復系數線性方程組通過雅可比預處理的QMR(quasi-minimal residual method)迭代法進行求解，提高了正演效率。Weiss等[8]利用有限差分法對三維各向異性介質的電磁感應問題進行了正演模擬，并采用全張量表示電導率參數，更有利于解決實際問題。Streich[9]利用MUMPS直接求解器對有限差分法離散頻率域可控源電磁方程得到的線性方程組進行求解，對于中等規(guī)模的求解問題精度高，且適用于多激發(fā)源的正演計算。Badea等[10]利用有限單元法計算了三維可控源的電磁響應，通過二次耦合勢的引入消除了源的奇異性，并通過局部網格加密提高了求解精度。Stalnaker[11]利用基于非結構化四面體網格的有限元對三維可控源電磁進行正演模擬，采用基于二次耦合勢的方式避免了直接模擬場源?；诜墙Y構化的四面體網格能夠適應更復雜的地形及地質體結構，且模擬精度很高。Farquharson等[12]采用矢量有限元進行大地電磁三維正演，并在各個單元接觸面上施加散度校正以提高正演計算速度。Puzyrev等[13]利用節(jié)點型有限單元法求解基于二次耦合勢的三維可控源電磁方程，并考慮了電導率的各向異性。Tong等[14]利用有限單元法對三維大地電磁法進行了正演模擬，通過將散度校正引入到有限元方程中，避免了偽解的出現。徐志峰等[15]基于磁矢量勢及電標量勢將Maxwell方程轉化為位勢的類Helmholtz方程，通過引入罰項及穩(wěn)定化方法避免了偽解及數值不穩(wěn)定，采用有限元方法計算了三維可控源的電磁響應。另外，付長民等[16]、殷長春等[17]、王超[18]、韓波等[19]、湯文武等[20]、嚴波等[21]、李建慧等[22]分別對三維可控源電磁正演問題開展了研究。

有限單元法由于其對地形及不規(guī)則異常體具有很好的適應性受到廣泛關注[23]。對于有限單元法在三維可控源電磁中的應用，早期主要采用節(jié)點型有限元，其概念及實現相對簡單，但可能存在偽解問題。近年來，矢量有限元法逐步獲得更多關注，究其原因在于矢量有限元在電導率分界面上自動滿足電磁場的邊界條件，可消除偽解。但目前并未見兩種正演方法實際應用的詳細對比。因此，基于前人研究成果，本文開展關于節(jié)點有限元及矢量有限元在三維頻率域可控源電磁正演計算的應用對比，旨在為兩種方法的合理選用提供參考。

2 理論基礎

頻率域可控源電磁法的三維正演可基于電場開展，亦可依據耦合勢進行，文中的三維正演是基于二次電場的有限元正演，可避免源的奇異性，與總場正演方式相比，同等條件下具有可采用更小網格的優(yōu)勢?；诠?jié)點有限元與矢量有限元的三維正演，其基本方程是一致的，但具體處理邊值問題時會有所差異。文中矢量有限元三維正演實現是基于筆者已發(fā)表的成果[20]，為了敘述完整，以下簡要討論其基本原理。

2.1 基于二次電場的偏微分方程

電磁場理論是可控源電磁法的基本理論。頻率域中，在準靜態(tài)條件下(忽略位移電流)，電磁場滿足以下麥克斯韋方程組

(1)

(2)

式中：E為電場矢量(單位為V/m)；H為磁場矢量(單位為A/m)；Jc為外電流密度矢量(單位為A/m2)；μ0為真空中的磁導率；σ為地下電性介質的電導率。在式(1)和式(2)中，時諧因子設為e-iω t。聯立方程式(1)和式(2)，可得到以下關于電場矢量的雙旋度方程

(3)

為了消除場源引起的奇異性，常將總場分解為一次場和二次場。通過預先選定某個簡單參考模型，利用計算得到的一次場Ep求解二次場Es，并最終得到總場，即

E=Ep+Es

(4)

對于一次場Ep，同樣可得以下雙旋度方程

(5)

式中σp為預先選定的參考模型的電導率分布，一般可選為全空間、均勻半空間或水平層狀介質模型以簡化一次場的計算。用式(3)減去式(5)，并化簡，得到

(6)

式中σa為異常電導率分布，且有σa=σ-σp。

2.2 邊界條件

采用一次場、二次場分離的方式計算總場時，邊界條件可以取得比較簡單。當邊界距離目標研究區(qū)域足夠遠時，二次場可近似為零。因此，采用如下狄利克雷邊界條件

n×Es=0

(7)

式(6)、式(7)便構成了可控源電磁三維正演的邊值問題。這里，選定均勻半空間模型計算一次場，其水平電偶極的諧波場存在解析積分公式[24]，具體表達式見附錄A。積分表達式中都是含有第一類貝塞爾函數J0或J1項的漢克爾變換，可利用數字線性濾波器快速求解，本文采用的是Guptasarma等[25]提出的數值濾波算法。

3 有限元分析

對于上節(jié)建立的邊值問題，可采用節(jié)點有限元或矢量有限元方法求解。為簡單起見，這里僅以長方體剖分單元為例展開討論。當采用節(jié)點有限元時，電磁場分量置于節(jié)點上采樣，而矢量有限元則將電磁場分量置于各個單元的棱邊上采樣(圖1)。圖1a為節(jié)點有限元中各個電場分量的分布，每個節(jié)點上都布置有電場三分量，因此每個單元有24個自由度；而圖2b為矢量有限元各個電場分量的布置，每條棱邊上都有與之平行的電場分量，每個單元的自由度為12。

圖1 電場分量在節(jié)點有限元(a)和矢量有限元(b)的剖分單元中的布置示意圖

對于以上邊值問題，基于變分原理可推導相應的變分問題如下

iωμ0σaEs·Ep]dV

(8)

(9)

式中λ是調節(jié)散度項在以上泛函表達式的比重因子，當取值為1時可獲得較滿意結果。

當采用矢量有限元時，由于矢量元在棱邊上取樣，在各個單元內自動滿足散度條件，因此不需要在各個剖分單元內施加散度條件。但進一步研究發(fā)現，基于矢量有限元求解時獲得的有限元方程的剛度矩陣條件數很大，采用迭代法收斂很慢，原因是在不同單元接觸面上可能不滿足散度條件，因此有必要進一步校正以加速迭代法求解過程。

4 模型實例對比

為了對比節(jié)點有限元與矢量有限元在三維可控源電磁中的模擬效果，設計了兩個模型開展正演計算，包括一個存在解析解的水平層狀介質模型及一個無解析解的低阻體模型。所有模擬結果都是設定兩種方法在同一剖分網格及線性方程組求解器的情況下獲得的。采用的計算平臺(laptop)為：Windows 7；CPU，Intel i7-3610qm；8GB RAM。

4.1 水平層狀介質模型

水平層狀介質模型的電磁響應具有解析表達式，這里選用如圖2所示的三層水平介質模型對正演計算結果進行對比驗證。場源取為單位電偶源并置于原點處，分別計算1～8192Hz按對數等間隔分布的14個頻率的電磁響應。采用的網格為91×76×50，其中z方向包含20個空氣層單元和30個地下介質層單元。圖3展示了坐標為(0，-2300m，0)的測深點處卡尼亞視電阻率及相位的解析解、節(jié)點有限元及矢量有限元卡尼亞視電阻率及相位的數值解，其中圖3a為卡尼亞視電阻率曲線的對比，圖3b為相位曲線的對比。

圖2 水平三層介質模型

由圖3a定性分析可知：矢量有限元獲得的卡尼亞視電阻率曲線與解析解吻合很好，而節(jié)點有限元獲得的卡尼亞視電阻率數值大體上與解析解吻合，但在個別低頻及高頻端數據誤差較大。進一步由計算結果統(tǒng)計可知，對于卡尼亞視電阻率而言，矢量有限元數值解所有頻點數據誤差在1%以下，平均誤差為0.17%；節(jié)點有限元數值解高頻、低頻段數據誤差大，中間頻段數據誤差小，基本在3%以下。對圖3b做類似分析可知，矢量有限元數值解整體上與解析解吻合較好，平均誤差為4.4%；而節(jié)點有限元在低頻、高頻段誤差大，中間頻段數據誤差小，基本在5%以下。

進一步通過模擬實驗可知，縮小中間網格剖分間距可以改善節(jié)點有限元在高頻段誤差大的現象，因此高頻段的誤差大主要跟網格剖分有關，相比矢量有限元而言需要更密的網格才能達到合理的精度。但低頻段的誤差通過網格剖分控制基本無法進一步改善，可能的原因應該是與節(jié)點有限元方法本身有關。

由于本文算法是基于二次場開展的，因此，為了更準確地對比矢量有限元與節(jié)點有限元的計算精度，有必要直接針對二次場進行對比。圖4所示是坐標為(0，-2300m，0)的測深點處的二次電場Exs的對比圖，其中二次場的解析解由圖2中三層介質模型及另一個電阻率為100Ω·m的均勻半空間模型的解析解之差獲得。圖4a代表二次電場分量Exs實部絕對值(以|Re(Exs)|表示)的對比，圖4b代表二次電場分量Exs虛部絕對值(以|Im(Exs)|表示)的對比。由圖可知，矢量有限元及節(jié)點有限元求解獲得的二次場與解析解基本吻合，但相比而言，矢量有限元精度更高。

圖3 卡尼亞視電阻率(a)及相位曲線(b)對比圖

圖4 二次電場實部分量(a)和虛部分量(b)曲線對比圖

4.2 單一低阻體模型

以一個單一低阻體模型為實例，分別進行基于節(jié)點有限元及矢量有限元法的三維可控源電磁響應模擬，并作對比。圖5是其斷面(圖5a)及平面(圖5b)示意圖。該模型描述了一個置于電阻率為100Ω·m的均勻半空間的長方體低阻體，激發(fā)源為單位電偶源，距離低阻體的最大水平距離為8.8km。該低阻體長、寬、高分別為0.6、0.5、0.2km，電阻率為10Ω·m。

統(tǒng)一采用91×116×50的網格對該模型進行正演計算，圖6分別給出了采用節(jié)點有限元及矢量有限元法的模擬結果。圖6a為節(jié)點有限元視電阻率不同頻率下的三維等值線切片圖； 圖6b為矢量有限元法正演結果對應的切片。5個不同切片分別對應64、32、16、8、4Hz。由頻率測深原理可知，5個不同切片分別對應由淺至深的地下介質的電阻率分布及變化情況。由于采用相同色標，因此顏色的變化即體現了視電阻率的變化。從圖中等值線分布情況來看，兩種正演方法得到的不同頻率下的視電阻率分布形態(tài)基本一致，因此相互驗證了各自正演算法的正確性。

具體分析各張切片圖：64Hz切片，盡管從等值線分布看兩種方法差別較大，但實際上視電阻率基本等于均勻大地的電阻率(100Ω·m)，此時趨膚深度最小，低阻體的影響基本可忽略； 從32Hz逐步過渡到8Hz，隨著趨膚深度增大，低阻體影響顯現出來，切片圖上可見藍色低值區(qū)，視電阻率最低約為86Ω·m； 4Hz對應的切片未見明顯低值區(qū)，視電阻率普遍高于100Ω·m，此時已基本進入近區(qū)。因此，分析不同頻率下切片圖的變化也在一定程度上驗證了正演算法的正確性。

對于高維電磁正演計算，有限元分析得到的往往是一個大型、稀疏、對稱復系數線性方程組，其求解是正演計算重點考慮的一個問題。采用預處理的QMR迭代法求解線性方程組時，統(tǒng)計了兩種不同算法在不同頻率的迭代殘差。如圖7所示，圖7a為節(jié)點有限元方程在不同頻率下的正演收斂曲線，圖7b為矢量有限元方程在不同頻率下的正演收斂曲線。

圖5 低阻體模型斷面(a)和平面(b)示意圖

圖6 節(jié)點有限元法(a)和矢量有限元法(b)模擬的不同頻率下視電阻率三維等值線切片圖

圖7 不同頻率下的迭代收斂曲線圖(a)節(jié)點有限元法； (b)矢量有限元法

分析圖7a可知，不同頻率的迭代收斂曲線基本一致，大約80次迭代可使得相對誤差降到10-6；分析圖7b可知，不同頻率下線性方程組求解的收斂曲線區(qū)別較大，總體規(guī)律是：頻率越低，迭代收斂越慢，達到同等精度需要更多迭代次數。因此，當采用迭代法求解有限元方程時，矢量有限元正演計算速度受頻率影響較大，且迭代收斂較節(jié)點有限元明顯要慢。表1給出了兩種正演計算在各個頻點的計算時間，總體而言，節(jié)點有限元在不同頻率下計算時間基本一致，而矢量有限元低頻時正演計算時間長。值得注意的是，盡管從自由度的角度來分析，矢量有限元的自由度小于節(jié)點有限元的自由度，同等條件下節(jié)點有限元離散得到的剛度矩陣更大。但由于節(jié)點有限元線性方程組迭代求解收斂較矢量有限元快，矢量有限元總體正演時間較節(jié)點有限元長，平均到各個頻率基本上矢量有限元計算速度較節(jié)點有限元慢一半。

表1 各頻率下正演計算時間對比

5 結論

本文對節(jié)點有限元和矢量有限元可控源電磁三維正演問題進行對比、分析，在控制并確保網格剖分一致并采用相同線性方程組求解器的情況下，得出以下認識和結論：

(1)同等情況下，矢量有限元求解精度較節(jié)點有限元精度高，在全頻率都能取得滿意的計算結果；而節(jié)點有限元在低頻段正演計算精度不高；

(2)采用迭代法求解有限元方程時，矢量有限元正演計算速度跟頻率有較大關系，通常來說，頻率越低，正演速度越慢；而強加散度條件的節(jié)點有限元正演計算速度與頻率關系不大；

(3)同等情況下，節(jié)點有限元正演計算速度較矢量有限元法快1倍以上。

總之，在計算效率優(yōu)先的情況下可考慮選擇節(jié)點有限元開展正演計算；當計算精度是第一指標時應當選用矢量有限元開展正演計算。

附錄A 水平電偶極源的諧波場

在水平電偶極源激發(fā)下，均勻半空間模型在地面上的電磁場響應可利用已有的積分公式計算，但由于基于二次場的三維正演模擬要利用地下空間的一次場，因此還需另外推導相應的積分公式。

考夫曼等[24]通過頻率域和時間域電磁測深法推導了均勻半空間模型上水平電偶極源在地面上的電磁場響應?；诖?，本文進一步推導了地下空間的電磁場響應的積分表達式。

圖A-1為一個均勻半空間模型，x軸沿正北方向，y軸沿正東方向，z軸豎直向下。時諧電流為I，長為dl的水平電偶極子位于坐標原點。均勻半空間的電導率為σb，設電磁場待求點為P，P與原點的距離為r，r與x軸的夾角為φ。假設時諧因子為e-iω t。均勻半空間電磁場矢量可由一個矢勢及一個標勢來表示(不考慮位移電流作用)

(A-1)

(A-2)

圖A-1 水平電偶極激發(fā)下的均勻半空間模型

(A-3)

(A-4)

且電場矢量可用矢勢單獨表示

(A-5)

依據文獻[24]的推導，假設Ay=0并求解式(A-4)，可得地下空間中矢勢的其他二分量積分表達式

(A-6)

(A-7)

(A-8)

采用圓柱坐標系，將式(A-6)～式(A-8)代入式(A-5)并整理，可得電場三分量表達式

(A-9)

(A-10)

(A-11)

同理，結合式(A-1)、式(A-6)和式(A-7)亦可得到磁場三分量表達式。

參考文獻

[1] 湯井田，何繼善．可控源音頻大地電磁法及其應用．湖南長沙：中南大學出版社，2005.

Tang Jingtian，He Jishan．The Controlled-Source Audio-frequency Magnetotellurics and Its Application．Central South University Press，Changsha，Hunan，2005.

[2] Nabighian M，Corbett J．Electromagnetic Methods in Applied Geophysics：Theory．SEG，1988.

[3] 楊懷杰，潘和平，孟慶鑫等．導電圍巖對井中三維瞬變電磁響應的影響規(guī)律研究．石油物探，2016，55(2)：288-293.

Yang Huaijie，Pan Heping，Meng Qingxin et al．Influence laws of conductive host on borehole 3D transient electromagnetic responses．GPP，2016，55(2)：288-293.

[4] 柳建新，麻昌英，孫麗影等．可控源音頻大地電磁測深法在地熱勘探中的應用．工程地球物理學報，2014，11(3)：319-325.

Liu Jianxin，Ma Changying，Sun Liying et al．The application of CSAMT to the geothermal exploration．Chinese Journal of Engineering Geophysics，2014，11(3)：319-325.

[5] He Z，Liu X，Qiu W et al．Mapping reservoir boundary by borehole-surface TFEM：Two case studies．The Leading Edge，2005，24(9)：896-900.

[6] Wirianto M，Mulder W A，Slob E C．A feasibility study of land CSEM reservoir monitoring in a complex 3-D model．Geophysical Journal of the Royal Astronomical Society，2010，181(2)：741-755.

[7] Newman G A，Alumbaugh D L．Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences．Geophysical Prospecting，1995，43(8)：1021-1042.

[8] Weiss C，Newman G．Electromagnetic induction in a generalized 3D anisotropic earth，Part 2：The LIN preconditioner．Geophysics，2003，68(3)：922-930.

[9] Streich R．3D finite-difference frequency-domain mode-ling of controlled-source electromagnetic data：Direct solution and optimization for high accuracy．Geophy-sics，2009，74(5)：F95-F105.

[10] Badea E A，Everett M E，Newman G A et al．Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials．Geophy-sics，2001，66(3)：786-799.

[11] Stalnaker J L．A Finite Element Approach to the 3D CSEM Modeling Problem and Applications to the Study of the Effect of Target Interaction and Topo-graphy [D]．Texas A & M University，2004.

[12] Farquharson C G，Miensopust M P．Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction．Journal of Applied Geophysics，2011，75(4)：699-710.

[13] Puzyrev V，Koldan J，Puente Jdl et al．A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling．Geophysical Journal International，2013，193(2)：678-693.

[14] Tong X Z，Liu J X，Xie W et al．Three-dimensional forward modeling for magnetotelluric sounding by finite element method．Journal of Central South University of Technology，2009，16(1)：136-142.

[15] 徐志鋒，吳小平．可控源電磁三維頻率域有限元模擬．地球物理學報，2010，53(8)：1931-1939.

Xu Zhifeng，Wu Xiaoping．Controlled source electromagnetic 3-D modeling in frequency domain for finite element method．Chinese Journal of Geophysics，2010，53(8)：1931-1939.

[16] 付長民，底青云，王妙月．海洋可控源電磁法三維數值模擬．石油地球物理勘探，2009，44(3)：358-363.

Fu Changmin，Di Qingyun，Wang Miaoyue．The 3D numerical modeling for MCSEM．OGP，2009，44(3)：358-363.

[17] 殷長春，賁放，劉云鶴等．三維任意各向異性介質中海洋可控源電磁法正演研究．地球物理學報，2014，57(12)：4110-4122.

Yin Changchun，Ben Fang，Liu Yunhe et al．MCSEM 3D modeling for arbitrarily anisotropic media．Chinese Journal of Geophysics，2014，57(12)：4110-4122.

[18] 王超．海洋可控源電磁法三維正演研究[學位論文]．北京：中國地質大學(北京)，2014.

Wang Chao．Marine CSEM 3D Forward Modeling [D]．China University of Geosciences (Beijing)，Beijing，2014.

[19] 韓波，胡祥云，Schultz A等．復雜場源形態(tài)的海洋可控源電磁三維正演．地球物理學報，2015，58(3)：1059-1071.

Han bo，Hu Xiangyun，Schultz A et al．Three dimensional forward modeling of marine controlled source electromagnetic field of complex source geometries．Chinese Journal of Geophysics，2015，58(3)：1059-1071.

[20] 湯文武，李耀國，柳建新等．基于二次電場的可控源電磁法三維矢量有限元正演模擬．石油物探，2015，54(6)：665-673.

Tang Wenwu，Li Yaoguo，Liu Jianxin et al．Three-dimensional controlled-source electromagnetic forward modeling by edge-based finite element using secondary electrical field．GPP，2015，54(6)：665-673.

[21] 嚴波，李予國，韓波等．任意方位電偶源的MCSEM電磁場三維正演．石油地球物理勘探，2017，52(4)：859-868.

Yanbo，Li Yuguo，Hanbo et al．The 3D forward mode-ling of the EM field for MCSEM using the electrical dipole source with arbitrary azimuth．OGP，2017，52(4)：859-868.

[22] 李建慧，胡祥云，陳斌等．復雜形態(tài)回線源激發(fā)電磁場的矢量有限元解．石油地球物理勘探，2017，52(6)：1324-1332.

Li Jianhui，Hu Xiangyun，Chen Bin et al．The vector finite element solution for the EM field induced by the loop source with complex form．OGP，2017，52(6)：1324-1332.

[23] 趙慧，劉穎，李予國．自適應有限元海洋大地電磁場二維正演模擬．石油地球物理勘探，2014，49(3)：578-585.

Zhao Hui，Liu Ying，Li Yuguo．The adaptive finite ele-ment forward modeling for 2D marine MT．OGP，2014，49(3)：578-585.

[24] 考夫曼 A A，凱勒 G V．頻率域和時間域電磁測深．北京：地質出版社，1987.

Kaufman A A，Keller G V．Frequency and Transient Sounding．Geological Publishing House，Beijing，1987.

[25] Guptasarma D，Singh B．New digital linear filters for Hankel J0and J1transforms．Geophysical Prospecting，1997，45(5)：745-762.