基于非結(jié)構(gòu)化網(wǎng)格的人工源頻率域三維電磁場數(shù)值模擬

熊彬，蔡紅柱，羅天涯，李長偉，丁彥禮，趙建國

(1.桂林理工大學(xué)a.地球科學(xué)學(xué)院;b.廣西礦冶與環(huán)境科學(xué)實驗中心，廣西桂林541004;2.College of Mines＆Earth Sciences，University of Utah，Salt Lake City，USA，84112;3.中國石油大學(xué)地球物理與信息工程學(xué)院，北京102249)

在減少數(shù)據(jù)解釋歧義、降低勘探風(fēng)險方面，可控源電磁法(CSEM)勘查一直以來都被認為是一種重要的地球物理手段[1-2]。不可否認的是，CSEM勘探方法效果的優(yōu)劣，將直接取決于人們對野外地質(zhì)環(huán)境電磁響應(yīng)的深入理解和認識[3]。隨著數(shù)值技術(shù)及計算機性能的快速發(fā)展，現(xiàn)在普通的臺式工作站已經(jīng)能容易地處理CSEM三維正演問題。就當(dāng)前而言，在電磁問題的正演模擬中，應(yīng)用最為廣泛的數(shù)值方法就是有限單元法和有限差分法[4]，盡管其他一些方法也伴隨著得到了長足的提高和完善，譬如積分方程法[5-7]、有限體積法[8]、對模型電導(dǎo)率施加某些額外約束的近似求解方法[9]、Neumann級數(shù)法[10]、人工神經(jīng)網(wǎng)絡(luò)法[11]，以及基于低對比度假設(shè)的Born-Rytov型近似求解[12]，等等[13]。

完全非結(jié)構(gòu)化三角形網(wǎng)格或四面體網(wǎng)格以Delaunay法或推進波陣面法為基礎(chǔ)，采用三角形(四面體)來填充二維或三維空間，它消除了結(jié)構(gòu)化網(wǎng)格中節(jié)點的結(jié)構(gòu)性限制，節(jié)點和單元的可控性好，能更好地處理邊界，適用于模擬不規(guī)則計算區(qū)域[14]，在刻畫復(fù)雜地質(zhì)構(gòu)造時有明顯優(yōu)勢。非結(jié)構(gòu)化網(wǎng)格在其生成過程中采用一定的準(zhǔn)則進行優(yōu)化判斷，能生成高質(zhì)量的網(wǎng)格，很容易控制網(wǎng)格的大小和節(jié)點的密度;此外，它采用隨機的數(shù)據(jù)結(jié)構(gòu)也有利于進行網(wǎng)格自適應(yīng)。一旦在邊界上指定網(wǎng)格的分布，邊界之間就可以自動生成網(wǎng)格，無需分塊或用戶的干預(yù)，而且不需要在子域之間傳遞信息。與此同時，有限單元法能實施完全非結(jié)構(gòu)化網(wǎng)格，使其單元邊界能與地下非均質(zhì)體的特有幾何形態(tài)達到一致，如偏離的鉆孔、流體浸入帶以及石油測井中經(jīng)常遇到的傾斜地層，正是這種自然支持非結(jié)構(gòu)化網(wǎng)格的特性，使得在處理具有一定程度幾何復(fù)雜性問題時，有限單元法成為了普遍性的首選數(shù)值方法[4，15]。

用有限單元法求解電磁感應(yīng)問題，可基于耦合的矢量位、標(biāo)量位(A，Ψ)，進而通過數(shù)值微分得到電磁場各分量，或直接求解電場與磁場矢量(E，H)，兩種方案各有優(yōu)勢[3，16-17]，本文將考慮基于位的CSEM三維有限元正演方案。

1 磁矢量位和電標(biāo)量位所滿足的微分方程

在Coulomb規(guī)范下，磁矢量位A和電標(biāo)量位Ψ在求解模型的邊界上連續(xù)，且A-Ψ方程對稱、簡單、通用。限于篇幅，下面直接列出基于Coulomb規(guī)范(▽·A=0)下磁矢量位A及電標(biāo)量位Ψ所滿足的二階微分方程[3-4]

及

式中:采用正弦諧變型時間因子e-iωt;ω為角頻率;μ0為自由空間磁導(dǎo)率;Js為場源電流分布;σ為介質(zhì)電導(dǎo)率。

為了避免式(1)和(2)中顯含場源項帶來的奇異性[18]，不妨定義關(guān)于電磁場二次位的表達式

由此不難獲得關(guān)于電磁場二次矢量位與標(biāo)量位的控制方程

這里

σp為背景電導(dǎo)率。另外，設(shè)定區(qū)域邊界遠離發(fā)射場源，以至邊界上的電磁場值可以忽略不計;此時，在外邊界上施加均勻Dirichlet邊界條件[19]

上述方程(4)、(5)和(7)即構(gòu)成了磁矢量位及電標(biāo)量位所滿足的三維邊值問題。

2 有限元分析

2.1 Galerkin線性方程組

很明顯，對于均勻全空間介質(zhì)中的偶極源與回線源，Ψp=0恒成立。在笛卡爾坐標(biāo)系中，將方程(4)和(5)沿三軸向依次展開

應(yīng)用Galerkin方法[20]，同時引用普通恒等式和散度定理，并注意到外邊界條件，容易得到關(guān)于標(biāo)量函數(shù)(Asx、Asy、Asz、Ψs)的體積分方程組

試驗函數(shù)N(r)可以選擇為基于節(jié)點的插值函數(shù)[20－21]。

2.2 總體系數(shù)矩陣的形成

用四面體單元離散整個計算區(qū)域，在各個單元中，電導(dǎo)率為常數(shù)，磁矢量位As及電標(biāo)量位Ψs呈線性變化

將上述兩式代入方程(12)～(15)，很明顯，這里的單元矩陣是一個分塊矩陣[4，21]，由4×4個子矩陣Kst(s，t=1，…，4)構(gòu)成

這里

如此逐單元分析，并將單元的系數(shù)矩陣擴展成由全體節(jié)點組成的系數(shù)矩陣，而后，將各個單元系數(shù)矩陣中的元加在總體系數(shù)矩陣數(shù)組中的相應(yīng)位置，即K=∑，由此，二階微分方程經(jīng)由離散最終得到如下線性代數(shù)方程組

b為一次位所形成的右端項列向量。

2.3 稀疏線性方程組的求解

方程組(20)中的K是一個大型稀疏對稱的非正定矩陣。鑒于目前計算機存貯空間及計算量方面的客觀需求，為了求解方程組(20)，即便采用現(xiàn)今最為高效的直接求解方法依然會是一個不小的挑戰(zhàn);再者，為便于算法日后的并行化，文中采用基于復(fù)系數(shù)Krylov子空間的迭代方法(QMR)進行大型稀疏線性方程組的求解[22]。

這里，diag()即為取總體系數(shù)矩陣對角線上的值。

3 均勻介質(zhì)中的一次位

根據(jù)文獻[25]，這里直接列出電導(dǎo)率為σ0的均勻介質(zhì)中一次磁矢量位及電標(biāo)量位的解析計算式

這里

式中:I為供電電流;d為偶極長度;xs、ys、zs為偶極源中心位置坐標(biāo)。

4 位的導(dǎo)數(shù)

上述有限單元法求解得到的是關(guān)于電磁場的二次位(As、Ψs)，為了進一步得到物理意義上的電磁場(Es、Hs)各分量，必須借助有效的數(shù)值微分方法來獲取

在規(guī)則化網(wǎng)格中，常規(guī)做法就是在節(jié)點周圍的單元中，分別計算場的導(dǎo)數(shù)，再通過面積加權(quán)平均或體積加權(quán)平均來獲得節(jié)點上場的導(dǎo)數(shù)[1];或者直接利用節(jié)點附近沿某一坐標(biāo)軸方向若干點上的場值，通過等距差商的辦法來獲取節(jié)點上場的偏微分值[21]。

針對完全非結(jié)構(gòu)化網(wǎng)格上的二次位，則采用滑動最小二乘插值的辦法[26-27]來獲取所期望節(jié)點t處電磁場各二次分量。不失一般性，設(shè)單元中場的線性描述

顯然，式中的系數(shù)pi(i=1，…，4)即為所需的場沿各坐標(biāo)軸的導(dǎo)數(shù)。略去中間推導(dǎo)過程，這里直接列出最終求取系數(shù)pi的方程

其中:

F即為期望節(jié)點t附近N個網(wǎng)格節(jié)點上的二次位的值。

5 數(shù)值算例

為了驗證文中所提方案的可靠性和有效性，利用COMSOL Multiphysics生成四面體單元網(wǎng)格，以3例典型地電模型從不同角度進行分析和闡述。具體計算平臺參數(shù)描述:Dell Precision T7610 CTO Base工作站，雙英特爾至強處理器E5-2697 v2(12核HT，2.7GHzTurbo，30 MB)。

5.1 均勻介質(zhì)中賦存單一目標(biāo)體

如圖1所示，均勻背景電導(dǎo)率為0.1 S/m，沿軸向布設(shè)水平電偶極子源，電偶極矩為100 Am，位于z=0界面以上100 m的(0，0，-100)處，發(fā)射頻率為10 Hz。目標(biāo)體電導(dǎo)率為0.01 S/m，大小為500 m×500 m×400 m，其幾何中心位于(0，0，800)處。

計算區(qū)域設(shè)置為Ω={-2 000 m≤x，y≤2 000 m;-1 000 m≤z≤2 000 m}，以美國猶他大學(xué)電磁法正反演科研組CEMI的積分方程法結(jié)果為參照，本文給出了對比結(jié)果，詳見圖2(為充分體現(xiàn)各分量的細節(jié)，這里分別考察了二次場的實部和虛部)。結(jié)果表明，兩者吻合情況令人滿意?；诖?，圖3、圖4分別給出了模型Ⅰ中二次電場和二次磁場實、虛部的平面分布。其中，由于受趨膚深度及網(wǎng)格剖分方面的影響，磁場分量Hsx出現(xiàn)一定程度的誤差。

圖2 有限單元解與積分方程解的對比Fig.2 Comparison of FE solution against IE solution

圖3 模型Ⅰ二次電場三分量實部與虛部平面圖Fig.3 Real and imaginary parts of secondary electric field for ModelⅠ

圖4 模型Ⅰ二次磁場三分量實部與虛部平面圖Fig.4 Real and imaginary parts of secondary magnetic field for ModelⅠ

5.2 多個目標(biāo)體的電磁響應(yīng)

圖5給出了一均勻背景中賦存有2個矩形異常體的情況，背景電導(dǎo)率值為0.01 S/m，同上，在z=0界面以上100 m的(0，0，100)處，沿軸向布設(shè)水平電偶極子源，電偶極矩為100 Am，發(fā)射頻率為10 Hz。兩目標(biāo)體電導(dǎo)率分別為0.1、0.001 S/m，大小均為400 m×400 m×200 m，其幾何中心分別位于(-800，0，1 000)、(800，0，600)處。

計算區(qū)域的設(shè)置及有限元網(wǎng)格的剖分均同上，仍以二次場各分量的實部和虛部為研究對象，考察各分量的固有特征及變化規(guī)律，模型Ⅱ中二次電場和二次磁場實、虛部的平面分布見圖6、圖7。

圖5 模型Ⅱ示意圖Fig.5 Map of ModelⅡ

圖6 模型Ⅱ二次電場三分量實部與虛部平面圖Fig.6 Real and imaginary parts of secondary electric field for ModelⅡ

圖7 模型Ⅱ二次磁場三分量實部與虛部平面圖Fig.7 Real and imaginary parts of secondary magnetic field for ModelⅡ

5.3 海底地形模型

圖8所示為儲藏在海丘之下的一塊藏油區(qū)，最大海水深度1 km，最淺500 m，海水電導(dǎo)率取3.3 S/m，沉積層水平和垂直方向電導(dǎo)率分別為1、0.8 S/m，油藏范圍:xy方向-1 000～1 000 m，z方向2 000～2 100 m，油藏區(qū)域水平和垂直方向電導(dǎo)率分別為0.05、0.005 S/m，場源位于(-3 000，0，500)m處，沿軸向布設(shè)水平電偶極子源，電偶極矩為100 Am，發(fā)射頻率為1 Hz。

如圖9所示，Ex二次場的實部在海丘處顯著增大，與5.1節(jié)相比較可知，海丘對Ex二次場的實部的影響遠大于藏油區(qū)的影響。

圖8 模型Ⅲ示意圖Fig.8 Map of ModelⅢ

圖9 模型Ⅲ二次場的實部Fig.9 Real part of secondary electric field for ModelⅢ

6 結(jié)束語

具有良好的通用性，適用于井中電磁、海洋電磁、環(huán)境地球物理等非均勻介質(zhì)中的電磁感應(yīng)基礎(chǔ)研究。在后續(xù)研究中，本課題組將在本文所提正演方案的基礎(chǔ)上，致力于解決人工源電磁法的三維反演問題。

美國猶他大學(xué)電磁法正反演科研組(CEMI)提供了優(yōu)良的科研環(huán)境，審稿人對本文提出了寶貴意見，在此一并致謝。

采用COMSOL生成完全非結(jié)構(gòu)化網(wǎng)格，本文建立了基于磁矢量位、電標(biāo)量位的人工源頻率域電磁法三維正演的流程與算法，通過與國際同行研究結(jié)果進行對比，文中正演算法是可靠和有效的。

由于Coulomb規(guī)范下磁矢量位、電標(biāo)量位在求解區(qū)域內(nèi)連續(xù)，滿足節(jié)點型有限元法的最基本假設(shè)，避免了直接求解時對場的強制連續(xù)性要求。采用非結(jié)構(gòu)化網(wǎng)格進行區(qū)域剖分，一方面可以在不影響解的精度的前提下，它能使得單元諸多屬性參數(shù)達到最大程度的優(yōu)化;另一方面，也是其不足之處，它會降低方程求解的收斂速度，針對這一問題，文中采用預(yù)條件因子改善稀疏矩陣條件數(shù)，大大提高了方程組求解收斂效率。本文算法具有良好的移植性，對方程組右端項稍加修改，即可適用于磁性源CSEM的三維正演。該方法也

[1]Nabighian M N.Electromagnetic Methods in Applied Geophysics—Theory(Volume 1)[M].Tulsa:Soc.Expl.Geophys.，1988.

[2]Nabighian M N.Electromagnetic Methods in Applied Geophysics:Volume 2，Application，Parts A and B[M].Tulsa:Soc.Expl.Geophys.，1991.

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

[4]Puzyrev V，Koldan J，de la Puente J，et al.A parallel finiteelement method for three-dimensional controlled-source electromagnetic forward modelling[J].Geophysical Journal International，2013，193(2):678-693.

[5]Zhdanov M S，F(xiàn)ang S.Quasi-linear series in three dimensional electromagnetic modeling[J].Radio Science，1997，32(6):2167-2188.

[6]Zhdanov M S.Geophysical Electromagnetic Theory and Methods[M].Oxford:Elsevier，2009.

[7]何展翔，王志剛，孟翠賢，等.基于三維模擬的海洋CSEM資料處理[J].地球物理學(xué)報，2009，52(8):2165-2173.

[8]Weiss C J，Constable S.Mapping thin resistors and hydrocarbons with marine EM methods，Part II—Modeling and analysis in 3D[J].Geophysics，2006，71(6):G321-G332.

[9]McKirdy D McA，Weaver J T，Dawson T W.Induction in a thin sheet of variable conductance at the surface of a stratified earth-II.Three-dimensional theory[J].Geophysical Journal International，1985，80(1):177-194.

[10]Avdeev D B，Kuvshinov A V，Pankratov O V，et al.Highperformance three-dimensional electromagnetic modelling using modified Neumann series.Wide band numerical solution and examples[J].Journal of Geomagnetism and Geoelectricity，1997，49:1519-1539.

[11]Spichak V，Popova I.Artificial neural network inversion of magnetotelluric data in terms of three-dimensional earth macroparameters[J].Geophys.J.Int.，2000，142:15-26.

[12]Song L P，Liu Q H.Fast three-dimensional electromagnetic nonlinear inversion in layered media with a novel scattering approximation[J].Inv.Prob.，2004，20(6):S171-S194.

[13]Avdeev D B.Three-dimensional electromagnetic modelling and inversion from theory to application[J].Surveys in Geophysics，2005，26:767-799.

[14]陳炎，曹樹良，梁開洪，等.結(jié)合前沿推進的Delaunay三角化網(wǎng)格生成及應(yīng)用[J].計算物理，2009，26(4):527-533.

[15]雷林林，劉四新，傅磊，等.海洋電磁法勘探油氣的三維正演模擬[J].世界地質(zhì)，2012，31(4):797-802.

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

[17]湯井田，任政勇，化希瑞.Coulomb規(guī)范下地電磁場的自適應(yīng)有限元模擬的理論分析[J].地球物理學(xué)報，2007，50(5):1584-1594.

[18]林昌洪，譚捍東，舒晴，等.可控源音頻大地電磁三維共軛梯度反演研究[J].地球物理學(xué)報，2012，55(11):3829-3839.

[19]Bгró O，Preis K.On the use of the magnetic vector potential in the finite element analysis of three-dimensional eddy currents[J].IEEE Trans.Magn.，1989，25(4):3145-3159.

[20]Jin J M.Finite Element Method in Electromagnetics[M].2nd ed.New York:John Wiley＆Sons，2002.

[21]徐世浙.地球物理中的有限單元法[M].北京:科學(xué)出版社，1994.

[22]Saad Y.Iterative Methods for Sparse Linear Systems[M].2nd Edition.Philadelphia:Society for Industrial and Applied Mathematics，2003.

[23]吳建平，王正華李曉梅.稀疏線性方程組的高效求解與并行計算[M].長沙:湖南科學(xué)技術(shù)出版社，2004.

[24]Ascher U M，Chen Greif.A First Course in Numerical Methods[M].Philadelphia:Society for Industrial＆Applied Mathematics，2011.

[25]Liu C S，Everett M E，Lin J，et al.Modeling of seafloor exploration using electric source frequency domain CSEM and the analysis of water depth effect[J].Chinese J.Geophys.，2010，53(4):669-683.

[26]Tabbara M，Blacker T，Belytschko T.Finite element derivative recovery by moving least square interpolants[J].Comp.Meth.Appl.Mech.Eng.，1994，117:211-223.

[27]Alexa M，Behr J，Cohen-Or D，et al.Computing and rendering point set surfaces[J].IEEE Trans.Vis.Comput.Graphics，2003，9(1):3-15.

[28]Schwarzbach C，B?rner R-U，Spitzer K.Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—A marine CSEM example[J].Geophysical Journal International，2011，187:63-74.