萬(wàn)文武, 李長(zhǎng)偉, 熊 彬, 高 磊, 宋衛(wèi)文

(桂林理工大學(xué) a.地球科學(xué)學(xué)院; b.廣西隱伏金屬礦產(chǎn)勘查重點(diǎn)實(shí)驗(yàn)室, 廣西 桂林 541006)

0 引 言

在過(guò)去的幾十年里, 地球物理電磁法取得了快速發(fā)展[1], 出現(xiàn)了大量的研究成果。傳統(tǒng)電磁正演主要包含了積分方程法(IEM)、 有限差分法(FDM)和有限元法(FEM)[2-6]。與FDM和FEM相比較, IEM的效率非常高, 主要是其離散的區(qū)域和方程組都非常小, IEM解的精度與對(duì)稱矩陣的稠密性息息相關(guān), 難以達(dá)到高精度。FDM使用交錯(cuò)有限差分法離散控制方程, 應(yīng)用于電磁場(chǎng)中, 計(jì)算簡(jiǎn)單、 編程容易實(shí)現(xiàn), 被廣泛應(yīng)用于電磁正演模擬中，其缺點(diǎn)在于無(wú)法使用非結(jié)構(gòu)化網(wǎng)格, 難以模擬復(fù)雜地質(zhì)體。FEM具有較好的靈活性, 能夠使用結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格。李長(zhǎng)偉等[7]利用仿射坐標(biāo)變換技術(shù), 應(yīng)用于平行多面體單元的單元分析, 極大地提高計(jì)算效率, 成為有效的電磁正演模擬方法, 特別是在復(fù)雜地形的正演研究中。熊彬等[8]借助COMSOL 生成完全非結(jié)構(gòu)化三維網(wǎng)格, 引入Coulomb 規(guī)范和適當(dāng)?shù)倪吔鐥l件保證矢量位唯一。在有限元方法中, 為提高求解精度需要對(duì)網(wǎng)格加密, 會(huì)導(dǎo)致計(jì)算量增加；另一種提高精度的策略是提高插值多項(xiàng)式的階數(shù), 即高階有限元法。Schwarzbach等[9]和Grayver等[10]運(yùn)用高階有限元結(jié)合面向目標(biāo)的局部網(wǎng)格細(xì)化, 可以達(dá)到較高的正演精度。

譜元法(spectral element method, SEM)最早是Patera在1984年提出的求解流體Navier-Stokes方程的一種數(shù)值方法[11]。 近年來(lái), SEM在計(jì)算地球物理中得到廣泛應(yīng)用, 如Seriani等 、 Komatitsch等提出基于Chebyshev和Legendre多項(xiàng)式的譜元法模擬地震波[12-13], 能夠模擬物性不連續(xù)介質(zhì)中的體波和面波的傳播特征。隨著開源包SPECFEM的出現(xiàn), 在地震研究中SEM成為解決局部和區(qū)域問(wèn)題的一種流行方法[14-16]。如Magnoni等[17]采用了SEM對(duì)2009年意大利中部發(fā)生的6.3級(jí)地震產(chǎn)生的復(fù)雜波場(chǎng)進(jìn)行模擬；王秀明等[18]引入基于預(yù)條件下共軛梯度法的元算法；劉有山等[19]和李琳等[20]基于Fekete和Cohen節(jié)點(diǎn)實(shí)現(xiàn)了譜元法采用三角網(wǎng)格剖分的地震模擬。在地球物理電磁法領(lǐng)域, 一些研究者將SEM應(yīng)用于航空電磁和可控源電磁法的正演模擬， 如Lee等[21-22]提出了基于Gauss-Lobatto-Legendre(GLL)多項(xiàng)式矢量電磁波方程和三維瞬變電磁問(wèn)題， Ren等[23]提出了一種新的時(shí)間域伽遼金SEM, 減少了未知數(shù)個(gè)數(shù)； Liu等[24]利用高階混合階SEM分析了各向異性和波導(dǎo), 在邊界上使用了旋度一致的GLL多項(xiàng)式和節(jié)點(diǎn)基于GLL的標(biāo)量基函數(shù)來(lái)計(jì)算正向和切向分量得到場(chǎng)值； Huang等[25]利用基于GLL插值多項(xiàng)式基函數(shù)的譜元法實(shí)現(xiàn)了頻域航空正演模擬； Zhou等[26]將基于GLL多項(xiàng)式的譜元法應(yīng)用到井地電磁中解決了低頻電磁場(chǎng)問(wèn)題； 劉玲等[27]實(shí)現(xiàn)了基于GLC插值基函數(shù)的譜元法在頻域三維海洋可控源電磁正演模擬研究。

與有限元相似, SEM也是將區(qū)域剖分成互不重疊的子單元, 在每個(gè)子單元內(nèi)基于譜方法進(jìn)行插值。 SEM能夠用結(jié)構(gòu)化或非結(jié)構(gòu)化網(wǎng)格模擬復(fù)雜的幾何形狀, 但與傳統(tǒng)的有限元相比, 譜元法具有指數(shù)收斂特性。SEM結(jié)合了有限元的靈活性和譜方法的高精度, 減少了內(nèi)存和CPU的需求。為提高大地電磁法正演計(jì)算的效率及精度, 本文將高精度的譜元法應(yīng)用到低頻平面波電磁場(chǎng)數(shù)值計(jì)算中。 在研究區(qū)域使用稀疏網(wǎng)格, 采用GLL多項(xiàng)式構(gòu)建插值基函數(shù), 利用GLL多項(xiàng)式的完全正交性, 形成對(duì)稱的質(zhì)量矩陣, 減少了計(jì)算時(shí)間, 并通過(guò)高階插值基函數(shù)來(lái)提高計(jì)算精度。與國(guó)際標(biāo)準(zhǔn)模型COMMEMI2D-0、 COMMEMI2D-1計(jì)算比對(duì), 驗(yàn)證了本文算法的可行性和有效性, 還設(shè)計(jì)了幾個(gè)含有構(gòu)造特征的典型地電模型, 分析其大地電磁響應(yīng)特征, 為實(shí)際野外地質(zhì)工作提供參考。

1 譜元法基礎(chǔ)理論

1.1 控制方程

采用時(shí)變因子e-iwt, Maxwell方程可寫為

(1)

其中:E為電場(chǎng)強(qiáng)度(V/m)；H為磁場(chǎng)強(qiáng)度(A/m)；B是磁感應(yīng)強(qiáng)度矢量(T)；ω為角頻率;σ為介質(zhì)的電導(dǎo)率;ε為介電常數(shù);μ為介質(zhì)的磁導(dǎo)率;ε、μ分別取真空中的介電常數(shù)和介質(zhì)磁導(dǎo)率, 即ε=ε0=1/36π×10-9(F/m)、μ=μ0=4π×10-7(H/m), 對(duì)Maxwell方程組(1)中的第一個(gè)方程取旋度, 得到E滿足的矢量Helmholtz方程

(2)

1.2 SEM基函數(shù)

采用高斯-洛巴托-勒讓德(Gauss-Lobatto-Legendre)多項(xiàng)式作為插值基函數(shù), 一維參考域的N階GLL插值多項(xiàng)式[21]的表達(dá)式

(3)

其中:LN(ξj)為L(zhǎng)egendre多項(xiàng)式;ξj為GLL插值節(jié)點(diǎn);j=0, 1,…,N。

對(duì)于二維等參單元, 其N階GLL插值基函數(shù)為

(4)

(5)

1.3 GLL插值基函數(shù)節(jié)點(diǎn)和階數(shù)的選擇

在SEM中, 選擇基函數(shù)和用于物理量插值的配置點(diǎn)是決定SEM法能否成功地關(guān)鍵。在標(biāo)準(zhǔn)單元中, 選擇不同函數(shù)展開形式對(duì)計(jì)算效率和數(shù)值精度都有很大影響。由于Legendre正交多項(xiàng)式的SEM, 在單元數(shù)值積分節(jié)點(diǎn)和函數(shù)插值節(jié)點(diǎn)選擇的是同一點(diǎn), 在數(shù)值計(jì)算過(guò)程中產(chǎn)生的質(zhì)量矩陣是對(duì)角陣, 具有簡(jiǎn)化數(shù)值計(jì)算和提高計(jì)算效率等優(yōu)點(diǎn), 從而本文選擇GLL積分。根據(jù)劉玲等[27]和Yin等[28]的研究, 當(dāng)插值階數(shù)小于4階時(shí), 必須通過(guò)增加單元個(gè)數(shù)才能保證SEM的高精度求解, 避免數(shù)值計(jì)算誤差較大問(wèn)題。理論上, 當(dāng)多項(xiàng)式插值階數(shù)無(wú)限增大時(shí), SEM計(jì)算的結(jié)果越逼近精確解, 但計(jì)算機(jī)內(nèi)存和計(jì)算時(shí)間相應(yīng)增加。通常采用4～8階多項(xiàng)式插值基函數(shù), 考慮到精度和計(jì)算量的平衡, 本文選擇了4階的GLL插值多項(xiàng)式。

提高SEM的精度可以通過(guò)網(wǎng)格單元數(shù)量和增加每個(gè)參考單元的插值節(jié)點(diǎn)或提高多項(xiàng)式基函數(shù)的階數(shù), 也就是說(shuō), 網(wǎng)格大小和插值基函數(shù)的階數(shù)直接決定了SEM的精度和耗時(shí)。在選取網(wǎng)格大小和基函數(shù)的階數(shù)時(shí), 網(wǎng)格尺寸越小和基函數(shù)階數(shù)越高, 隨之精度越高, 越接近數(shù)值解,但需要計(jì)算機(jī)的內(nèi)存越大和計(jì)算時(shí)間越長(zhǎng)。網(wǎng)格剖分在SEM計(jì)算中是一個(gè)極其關(guān)鍵的步驟, 需要不斷重復(fù)調(diào)試適合的網(wǎng)格單元(包括網(wǎng)格形狀、 網(wǎng)格數(shù)量、 網(wǎng)格尺寸)。

1.4 變分方程離散化

將計(jì)算區(qū)域Ω離散成一組互不重疊的子區(qū)域Ωe。對(duì)于任意四面體單元Ωe, 利用一個(gè)等參變換將其轉(zhuǎn)換成標(biāo)準(zhǔn)參考單元, 即設(shè)物理單元坐標(biāo)和標(biāo)準(zhǔn)參考單元坐標(biāo)之間的映射函數(shù)為Fe:Λ?Ω, 對(duì)剛度矩陣和質(zhì)量矩陣使用同樣的映射函數(shù)轉(zhuǎn)化到參考域中計(jì)算。物理單元與參考單元之間的轉(zhuǎn)換如圖1所示, 右邊為物理坐標(biāo)(x,y)下的四邊形單元, 左邊是標(biāo)準(zhǔn)參考(ξ,η)下的四邊形單元,ξ,η∈[-1, 1]。

使用協(xié)變量映射[29]表達(dá)式計(jì)算參考域下基函數(shù)

(7)

(8)

(9)

將微分方程(2)轉(zhuǎn)化為積分問(wèn)題, 采用伽遼金

圖1 物理單元與標(biāo)準(zhǔn)參考單元的轉(zhuǎn)換

加權(quán)殘差法, 令加權(quán)殘差為0, 即

?Ωevi·(×μ-1(×E)-iω(σ-iωε)E)ds=0,

(10)

其中,vi為加權(quán)函數(shù)。

基于格林第一定理和狄利克雷邊界條件n×E|?L=0, 把式(10)改寫成

?Ωe(μ-1×Vi)T·(×E)ds-

(11)

(12)

其中,Na是基函數(shù)的個(gè)數(shù)。對(duì)單元進(jìn)行集成, 得到最終的線性方程組

(A-B)E=0,

(13)

其中:A為剛度矩陣;B為質(zhì)量矩陣。

(14)

(15)

其中,k是單元總數(shù)。對(duì)于磁場(chǎng)H計(jì)算與上述類似的方法對(duì)Maxwell方程組中第二個(gè)方程取旋度得到, 使用了穩(wěn)定雙共軛梯度算法求解上述線性方程組, 可得到電磁場(chǎng)的分布。

2 算例分析

為了驗(yàn)證算法的正確性和有效性, 本文設(shè)計(jì)了從簡(jiǎn)單到復(fù)雜的地電模型進(jìn)行測(cè)試, 采用4階GLL多項(xiàng)式作為插值基函數(shù), 在如下的微型計(jì)算機(jī)平臺(tái)上進(jìn)行運(yùn)算:處理器 Intel(R) Core(TM)i5-4210 CPU@1.70 GHz 2.40 GHz, 內(nèi)存8 GB, Windows 10 64位操作系統(tǒng)。

2.1 算法驗(yàn)證

算例1采用國(guó)際標(biāo)準(zhǔn)模型庫(kù)中的COMMEMI 2D-0和COMMEMI 2D-1模型[30]對(duì)本文算法進(jìn)行驗(yàn)證, 模型如圖2所示。

圖2 模型COMMEMI 2D-0(a)及COMMEMI 2D-1(b)示意圖

在圖2a中，x∈[-50, -10]區(qū)域上電阻率為10 Ωm、x∈[-10, 10]區(qū)域上電阻率為1 Ωm、x∈[10, 50]上電阻率為2 Ωm。頻率取1/300 Hz, 計(jì)算地面上x=-25、 -15、 -7、 7、 15、 30 km處的視電阻率, 本文算法與積分方程法的結(jié)果對(duì)比如圖3所示, 兩種方法都能夠反映出垂直地層電阻率的變化趨勢(shì), 結(jié)果幾乎一致, 說(shuō)明了SEM算法用于2DMT正演的正確性。

圖3 IE與SEM的視電阻率

圖2b是在電阻率為100 Ωm的均勻半空間中嵌入一個(gè)電阻率為0.5 Ωm的異常體, 異常體在x方向上的位置是-0.5～0.5 km, 長(zhǎng)為2 km, 寬為1 km, 上頂面離地面 0.25 km。 計(jì)算地面上x=0.5、 1、 2、 4、 8、 16 km處的視電阻率。頻率取0.1和10 Hz時(shí)的SEM與IE的計(jì)算結(jié)果對(duì)比如圖4所示, 兩種方法在不同頻域時(shí)的模擬結(jié)果幾乎一致, 誤差小于5%, 能夠反映出地下低阻異常體的特征, 說(shuō)明了SEM算法用于2DMT正演計(jì)算的準(zhǔn)確性。

2.2 二維模型譜元法效率分析

分別采用本文算法(SEM)和雙二次插值的有限元算法(FEM)計(jì)算模型COMMEMI 2D-0, 對(duì)二者的計(jì)算結(jié)果進(jìn)行了對(duì)比。實(shí)驗(yàn)中設(shè)計(jì)了由粗到細(xì)4套不同的網(wǎng)格, 在粗網(wǎng)格1上應(yīng)用本文算法的4階譜元法, 在逐漸細(xì)化的網(wǎng)格上應(yīng)用FEM進(jìn)行計(jì)算。

圖4 IE與SEM在頻率0.1(a)和10 Hz(b)的視電阻率

表1給出了SEM和FEM的自由度和計(jì)算時(shí)間比較, 圖5是FEM在不同尺寸的網(wǎng)格(Grid2、Grid4)和SEM、 IEM的計(jì)算結(jié)果比對(duì)。結(jié)果表明, 相同網(wǎng)格尺寸下SEM比FEM(Grid2)精度高、 耗時(shí)長(zhǎng), 但在同等精度下(FEM(Grid4))SEM所用計(jì)算時(shí)間更少。因此與FEM相比, SEM在粗網(wǎng)格下能夠保證精度, 減少計(jì)算時(shí)間。

表1 譜元法與有限元法的計(jì)算時(shí)間對(duì)比

圖5 IE與SEM(Grid1)、 FEM(Grid2/Grid4)的視電阻率

2.3 水平地層中的異常體模擬

圖6 層狀地層低阻異常體(a)及高阻異常體(b)示意圖

圖7 低阻異常體(a)及高阻異常體(b)視電阻率等值線圖

2.4 復(fù)雜地形模擬

2.4.1 地壘模擬 為研究起伏地形條件下平面波電磁響應(yīng), 算例3設(shè)計(jì)了一個(gè)地壘模型, 如圖8所示。在地壘模型中, 隆起部位上端寬為10 km, 下端寬為20 km, 凸起垂直高度為5 km, 地下均勻半空間的電阻率為100 Ωm, 空氣層電阻率取1010Ωm, 頻率為f=0.1、 1和100 Hz, 計(jì)算地面上的視電阻率和相位, 結(jié)果如圖9所示。起伏地形對(duì)視電阻率和相位均發(fā)生了畸變, 視電阻率隨著地形的凸起而增大, 且頻率越大響應(yīng)越明顯, 與視電阻率相比, 相位受地形的影響更為明顯。

圖8 地壘模型示意圖

2.4.2 覆蓋層下的垂直斷層模擬 算例4設(shè)計(jì)了覆蓋層下有不同巖性的垂直斷層帶, 如圖10所示。上層覆蓋層電阻率為50 Ωm； 第二層是直立的斷層帶, 電阻率依次為30、 1、 100、 30 Ωm； 第三層為沉積層，其電阻率為100 Ωm； 第四層為基巖, 電阻率為0.1 Ωm。在頻率0.01～10 Hz范圍內(nèi)進(jìn)行正演計(jì)算, 得到如圖11所示的視電阻率分布。

可看出SEM能夠較好反映覆蓋層下的斷層, 在頻率為1 Hz的時(shí)候反映了覆蓋層的電阻率情況, 探測(cè)深度隨著頻率的減小而增大; 頻率在0.15 Hz時(shí)的視電阻率響應(yīng)接近垂直斷層, 在x方向上-60～-20 km和30～60 km顯示了30 Ωm的異常, 在-20～10 km能夠準(zhǔn)確得到低阻異常(1 Ωm), 在10～30 km清晰反映了高阻異常(100 Ωm); 再次減小頻率, 得到沉積層和基巖電阻率, 但在-15～10 km處受到低阻斷層異常的影響, 視電阻較低, 說(shuō)明在層狀地層深部受到低阻斷層的影響較大, 而在覆蓋層中受到高阻斷層的影響更加明顯, 說(shuō)明在層狀地層中淺部受到高阻異常的影響較大。

圖9 地壘視電阻率(a)及相位(b)

圖10 覆蓋層下的垂直斷層示意圖

2.4.3 簡(jiǎn)單斷層與異常體模擬 算例5給出了一個(gè)更為復(fù)雜的地電模型(圖12), 其背景電阻率為100 Ωm, 存在電阻率為10 Ωm的4個(gè)低阻體, 分別位于x∈[0,14]、z∈[8,12];x∈[8,12]、z∈[26,30];x∈[20, 40]、z∈[8, 12], 在x為20～30 km是傾斜的;x∈[30,34]、z∈[22,30]。在其深度為2～3 km處有一條高阻斷層, 電阻率為1 000 Ωm，厚度為0.5 km, 在x方向20 km處向下錯(cuò)動(dòng)0.4 km。在頻率0.01～1 Hz范圍內(nèi)計(jì)算得到頻率電阻率等值線圖(圖13)及不同頻率的相位曲線(圖14)。

圖11 頻率視電阻分布

圖12 斷層巖體下的異常體示意圖

在圖13中電阻率等值線圖能夠精確反映地電結(jié)構(gòu), 清晰地顯示出高阻斷層和低阻體的存在, 右邊低阻傾斜體沒(méi)有完全吻合的原因是由于使用階梯網(wǎng)格和網(wǎng)格尺寸較大所致。

圖14顯示了不同頻率的相位變化曲線, 相位曲線更加清晰的反映了各個(gè)深度的地電結(jié)構(gòu)。在地表1 Hz的相位清晰地反映了高阻斷層巖, 0.25 Hz的相位能夠明顯反映淺部?jī)蛇叺妥梵w和中間的背景場(chǎng)在0.063 Hz體現(xiàn)了上下低阻體中間位置的電場(chǎng)變化, 在低阻體所在位置均有變化趨勢(shì), 在右端低阻體相位變化趨勢(shì)比左端低阻體相位變化平緩, 是由于受到右端傾斜低阻體的影響。 在0.037 Hz相位(圖14b)能夠清晰地反映底部低阻體和背景場(chǎng), 能夠體現(xiàn)兩個(gè)低阻體在x方向的準(zhǔn)確位置; 在0.024 Hz仍能夠反映低阻體的存在; 在0.01 Hz相位趨于平緩, 與背景場(chǎng)趨于一致。

圖13 視電阻率等值線圖

2.4.4 深部地電結(jié)構(gòu)簡(jiǎn)化模擬 為研究深層流體穿過(guò)高阻巖石圈形成傳導(dǎo)通道時(shí)的電磁響應(yīng)特征[31-33], 設(shè)計(jì)了算例6，如圖15所示。 第一層為水平層低阻體, 層厚1 km, 電阻率為30 Ωm； 第二層為高阻體, 電阻率為2 000 Ωm, 層厚20 km, 含有2個(gè)垂直狹窄的通道ρ1(其高度為7 km), 下半部分是一個(gè)地壘結(jié)構(gòu)的低阻體, 電阻率為10 Ωm, 模型的類型取決于ρ1的選擇。觀測(cè)點(diǎn)取傳導(dǎo)通道中心點(diǎn)O和離中心點(diǎn)25 km的點(diǎn)R, 在頻率0.01～10 Hz范圍內(nèi)進(jìn)行計(jì)算，得到存在流體通道(ρ1區(qū)域)和不存在流體通道在不同觀測(cè)點(diǎn)處的視電阻率，如圖16所示。

圖14 不同頻率(1 、 0.25、 0.063、 0.037、 0.024、 0.01 Hz)相位曲線

圖15 深部地電結(jié)構(gòu)簡(jiǎn)化模型示意圖

圖中實(shí)線表示ρ1=2 000 Ωm, 即不存在垂直的流體通道情況下, 在中心點(diǎn)O和遠(yuǎn)點(diǎn)R處的視電阻曲線, 低頻時(shí)無(wú)論觀測(cè)點(diǎn)R還是O都準(zhǔn)確地反映了深部地電結(jié)構(gòu), 在高頻時(shí), 遠(yuǎn)點(diǎn)R的視電阻率高于中心點(diǎn)O, 這是因?yàn)橹行狞c(diǎn)O處的視電阻率受到底部低阻地壘的影響; 虛線表示ρ1=10 Ωm, 即存在一個(gè)雙向低阻溶液流體通道, 在中心點(diǎn)O和遠(yuǎn)點(diǎn)R處的視電阻率曲線, 視電阻率曲線和無(wú)導(dǎo)電通道的形態(tài)相似, 但在中心點(diǎn)O處, 視電阻率低于無(wú)流體通道時(shí)的視電阻率值, 在遠(yuǎn)點(diǎn)R處視電阻率與無(wú)流體通道時(shí)的值一致。 由分析可知, 大地電磁響應(yīng)對(duì)低阻流體通道有相當(dāng)高的靈敏度, 這些通道將近地表和深部導(dǎo)體連接起來(lái), 形成閉合的電流電路, 導(dǎo)致地殼和上地幔中的大地電流重新分布。

圖16 視電阻率對(duì)比

3 結(jié) 論

本文實(shí)現(xiàn)了基于譜元法的低頻平面波電磁場(chǎng)的數(shù)值模擬。在數(shù)值實(shí)驗(yàn)中, 對(duì)國(guó)際標(biāo)準(zhǔn)模型COMMEMI2D-0、 COMMEMI2D-1進(jìn)行計(jì)算, 驗(yàn)證了算法的有效性和準(zhǔn)確性, 與有限元方法對(duì)比, 表明了本文算法具有更高的計(jì)算精度和效率, 通過(guò)設(shè)計(jì)復(fù)雜地電模型, 表明了譜元法應(yīng)用于大地電磁正演的有效性, 并進(jìn)一步分析了這些典型模型的大地電磁響應(yīng)特征, 對(duì)實(shí)際工作提供了借鑒和指導(dǎo)意義。為進(jìn)一步提高算法的實(shí)用性和計(jì)算效率, 將在以后的研究中考慮將有限元和譜元法相結(jié)合。